II. Инструменты геоанализа (1)#
В этом занятии мы познакомимся с базовыми методами пространственного анализа данных, используя возможности библиотек GeoPandas и Shapely. Мы научимся строить буферные зоны, выполнять пространственные объединения (spatial join), и пространственную фильтрацию, а также применять эти инструменты на практике: агрегировать точечные данные по ячейкам регулярной сетки.
Эти навыки подготовят нас к анализу исторической вспышки холеры в Лондоне (1854) – второму проекту на нашем курсе.
Импортируем библиотеки
Первым делом подключим необходимые библиотеки
import geopandas as gpd
import osmnx as ox
from shapely.geometry import Point, Polygon, box
import matplotlib.pyplot as plt
GeoPandas – надстройка над pandas, расширяющая возможности DataFrame для работы с геометрией. Хранит объекты геометрии (точки, линии, полигоны) в специальном столбце `geometry` и предоставляет методы для пространственных операций (буферизация, объединение, объединение по пространственному признаку).
osmnx – библиотека для загрузки и анализа данных OpenStreetMap. Позволяет легко получать графы улиц, точки интереса и другие геоданные напрямую из OSM, а также выполнять маршрутизацию, визуализацию и базовый пространственный анализ городских сетей.
Shapely – библиотека для создания и операций над геометрическими примитивами (точки, линии, полигоны). Предоставляет методы для объединения, пересечения, измерения расстояний и вычисления буферных зон.
matplotlib.pyplot – модуль визуализации, вдохновлённый MATLAB. Позволяет строить 2D-графики (линии, гистограммы, scatter-плоты), настраивать оформление, добавлять подписи, легенды и работать с осями.
contextily – утилита для подложек карт (tile basemaps) в matplotlib. Позволяет легко добавить тайловые слои (OpenStreetMap, Stamen и др.) под ваши геоплоты, автоматически приводя проекции в соответствие.
1. Базовые инструменты пространственного анализа#
Прежде чем переходить к основным инструментам, рассмотрим простые операции с геоданными: повторим как мы работаем с разными системами координат, рассмотрим измерение расстояний, буферизацию и пространственные отношения между объектами
1.1 Создание точек и перепроецирование (повторение)#
Для дальнешей работы создадим две точки на основе координат в Санкт-Петербурге
# Задаем координаты двух точек (долгота, широта)
lon1, lat1 = 30.3359, 59.9316 # Основное здание Александринского театра (историческая сцена)
lon2, lat2 = 30.3322, 59.9385 # Государственный Русский музей (Михайловский дворец)
# Создаем GeoSeries из этих точек с указанием CRS WGS84 (EPSG:4326 - географическая система координат)
points = gpd.GeoSeries([Point(lon1, lat1), Point(lon2, lat2)], crs="EPSG:4326")
points
0 POINT (30.33590 59.93160)
1 POINT (30.33220 59.93850)
dtype: geometry
Посмотрим на точки на карте
points.explore(tiles='CartoDB positron')
Мы создали две точки в географической системе координат WGS84 – широта и долгота в градусах.
Однако, как мы помним, для многих аналитических задач (например, расчета расстояний, площадей) предпочтительно использовать систему координат в метрических единицах.
Проверим текущую систему координат наших данных
current_crs = points.crs
print("Текущая система координат:", current_crs)
Текущая система координат: EPSG:4326
Оценим подходящую UTM-зону для наших точек и перепроецируем в нее
utm_crs = points.estimate_utm_crs() # автоматически определяем UTM-зону
points_utm = points.to_crs(utm_crs) # преобразуем GeoSeries в эту проекцию
print("Новая система координат:", points_utm.crs)
Новая система координат: EPSG:32636
Теперь обе точки заданы в спроецированной системе координат (UTM Zone, в метрах). Вычислим расстояние между ними.
Сначала попробуем вычислить расстояние напрямую в географической CRS (в градусах), а затем в метрической CRS:
# Расстояние между точками в градусах (как есть, в исходной CRS)
dist_deg = points[0].distance(points[1])
# Расстояние между точками в метрах (после проекции в UTM)
dist_m = points_utm[0].distance(points_utm[1])
print(f"Расстояние в градусах: {dist_deg:.6f}°")
print(f"Расстояние в метрах: ~{dist_m:.0f} м")
Расстояние в градусах: 0.007829°
Расстояние в метрах: ~796 м
В дальнейшем, всегда перед измерением расстояний или площадей, убедимся, что данные заданы в метрической системе координат.
1.2 Буферные зоны (Buffer)#
Буферная зона – это область вокруг геометрического объекта на заданном расстоянии. Чаще всего буфер строят вокруг точек, линий или полигонов, чтобы выделить зону влияния или ближайшего окружения.
При работе с буферами важно помнить про систему координат: в географической системе координат (широта/долгота) радиус буфера указывается в градусах. Поэтому для измерения реальных расстояний буфер можно строить после в метрических системах координат (например, UTM).
buffer_500 = points_utm.buffer(500)
buffer_500.explore(tiles='CartoDB positron')
Построение буферных зон — относительно простая операция, но она регулярно применяется в пространственном анализе. При ее выполнении особенно важно внимательно отслеживать единицы измерения в используемой системе координат.
1.3 Пространственное объединение (Spatial Join)#
Пространственное объединение (spatial join) — операция, при которой к каждому объекту одного слоя присоединяются атрибуты объекта другого слоя, с которым он связан пространственно. Это может быть принадлежность к району, попадание в буфер и т.д.
Однако само по себе пространственное объединение — это не только механическое добавление атрибутов на основе пересечений геометрий. Чтобы этот процесс был действительно точным и осмысленным, важно понимать, каким именно образом соотносятся объекты в пространстве. Пересекаются ли они? Один внутри другого? Или, быть может, лишь касаются границ? Для того чтобы корректно интерпретировать такие случаи, используются пространственные предикаты — ключевой инструмент в пространственном анализе, позволяющий задать и проверить топологические отношения между геометриями.
Пространственные предикаты – это функции, позволяющие определить топологическое отношение между двумя геометрическими объектами (точками, линиями, полигонами и т.д.). Предикаты позволяют отвечать на вопросы вроде: «Пересекаются ли эти объекты?», «Находится ли один объект внутри другого?», «Касаются ли они границ друг друга?» и т.п. Знание пространственных отношений помогает выполнять пространственные запросы и объединения данных (spatial join), фильтрацию на основе взаимного расположения и др.
Ниже приведена таблица основных бинарных предикатов (отношений) с кратким описанием:
Предикат |
Описание (условие истинности) |
---|---|
intersects (пересекается) |
Имеют общую точку – Геометрии пересекаются, если они имеют хотя бы одну общую точку в пространстве (любое наложение границ или внутренних областей). Эквивалентно логическому отрицанию disjoint (см. ниже). |
disjoint (не пересекается) |
Не имеют общих точек – Геометрии не пересекаются, если у них нет ни одной общей точки. То есть полностью разделены в пространстве. |
within (внутри) |
A внутри B – Геометрия A находится целиком внутри геометрии B, не выходя за ее границы. Все точки A (и граница, и внутренняя часть) лежат во внутренней области B. (Эквивалентно B.contains(A)). |
contains (содержит) |
A содержит B – Геометрия A полностью содержит геометрию B. То есть ни одна часть B не лежит вне A, при этом B не просто граничит с A, а имеет хотя бы одну точку строго внутри A. (Эквивалентно B.within(A)). |
touches (касается) |
Касание границ – Геометрии имеют хотя бы одну общую точку на границе, но не пересекаются по внутренней области. Иными словами, их границы соприкасаются, но пересечения интерьеров нет. |
overlaps (перекрывает) |
Частичное перекрытие – Геометрии частично перекрываются, если они имеют общую часть, но ни одна не содержит другую целиком. При этом геометрии одного типа (полигон-полигон или линия-линия) – например, два полигона, которые пересекаются по площади, но каждый выходит за пределы другого. |
crosses (перекрещивается) |
Пересечение с выходом – Геометрии пересекаются внутренними точками таким образом, что пересечение имеет меньшую размерность. Например, линия пересекает полигон, входя и выходя из него (часть линии внутри, часть снаружи), или две линии пересекаются в одной точке. |
Какой предикат нам понадобится для пространствнного объединения точек метро с полигонами районов города?
(ответ будет ниже по ходу работы)
Перейдем к задаче – определим, в каких районах города находятся станции метро
#Прочитаем данные о станциях метро
metro = gpd.read_file('data/spb_metro.geojson')
#Прочитаем данные о городских районах
admin_district = gpd.read_file('data/spb_admin.gpkg', layer="district")
C помощью пространственного объединения опрееделим, какая станция находится в каком районе
metro_with_district = gpd.sjoin(metro, admin_district, predicate='within')
metro_with_district[['name', 'geometry', 'NAME']].head()
name | geometry | NAME | |
---|---|---|---|
0 | пр.Ветеранов-1 | POINT (30.25319 59.84135) | Кировский район |
1 | пр.Ветеранов-2 | POINT (30.25056 59.84150) | Кировский район |
2 | Ленинский пр.-1 | POINT (30.26675 59.85135) | Кировский район |
3 | Ленинский пр.-2 | POINT (30.26951 59.84956) | Кировский район |
4 | Автово | POINT (30.26136 59.86731) | Кировский район |
Понимание пространственных предикатов особенно важно не только для объединения данных, но и для выполнения фильтрации объектов на основе их взаимного расположения. Иногда нам вовсе не нужно присоединять атрибуты одного слоя к другому — достаточно отобрать те объекты, которые, например, находятся в пределах заданной области, пересекаются с буферной зоной или касаются границ интересующей территории. Такая фильтрация позволяет гибко и точно выделять пространственно значимые подмножества данных, что делает её неотъемлемым инструментом пространственного анализа.
1.4 Фильтрация объектов по пространственному признаку#
Одно из распространённых действий в пространственном анализе – отбор объектов, удовлетворяющих определённому пространственному условию. Например, можно захотеть выбрать все точки (объекты) внутри заданного полигона. Это называется пространственной фильтрацией (фильтрация по местоположению).
Давайте отфильтруем все кафе Санкт-Петербурга, которые находятся в зоне 500-метровой доступности от станций метро
Подготовим данные
#Прочитаем данные о станциях метро
metro = gpd.read_file('data/spb_metro.geojson')
Построим буферы 500-метровой доступности вокруг станций (предварительно проверим систему координат и перепроецируем, если нужно)
metro.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Данные в географической системе координат, нам это не подходит, поэтому перепроецируем их перед тем как строить буферы
utm_crs = metro.estimate_utm_crs() # автоматически определяем UTM-зону
metro_utm = metro.to_crs(utm_crs)
Построим буферы
metro_buff_500 = metro_utm.buffer(500)
metro_buff_500.explore(tiles='CartoDB positron')
Выгрузим информацию о кафе из OSM, отфильтруем только точки и перепроецируем
tags = {"amenity": "cafe"}
cafes = ox.features_from_place("Санкт-Петербург, Россия", tags)
#фильтрация точек
cafes = cafes[cafes.geometry.type == 'Point']
#перепроецирование
cafes_utm = cafes.to_crs(utm_crs)
cafes_utm.explore(tiles='CartoDB positron')
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:717: DeprecationWarning: Passing a BlockManager to GeoDataFrame is deprecated and will raise in a future version. Use public APIs instead.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:717: DeprecationWarning: Passing a BlockManager to GeoDataFrame is deprecated and will raise in a future version. Use public APIs instead.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:717: DeprecationWarning: Passing a BlockManager to GeoDataFrame is deprecated and will raise in a future version. Use public APIs instead.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:717: DeprecationWarning: Passing a BlockManager to GeoDataFrame is deprecated and will raise in a future version. Use public APIs instead.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:717: DeprecationWarning: Passing a BlockManager to GeoDataFrame is deprecated and will raise in a future version. Use public APIs instead.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:717: DeprecationWarning: Passing a BlockManager to GeoDataFrame is deprecated and will raise in a future version. Use public APIs instead.
warnings.warn(