(Пере)проецирование пространственных

В предыдущем разделе мы познакомились с основами картографических проекций и систем координат.

В этом разделе мы рассмотрим, как:

выбирать подходящую проекцию для пространственного анализа;
перепроецировать данные в другую систему координат;
назначать CRS, если она отсутствует;

0. Импорт библиотек и подготовка данных¶

0.1 Импорт библиотек¶

import geopandas as gpd
import pandas as pd
import osmnx as ox

0.2 Подготовка данных¶

Загрузим границу города из OpenStreetMap

area_name = "Центральный район, Санкт-Петербург"

admin_border = ox.geocode_to_gdf(area_name)
admin_border.explore(tiles="cartodbpositron")

1. Проверка системы координат¶

При работе с пространственными данными одним из первых шагов является проверка системы координат, в которой они находятся. Это позволяет понять, как интерпретировать координаты и подходят ли данные для дальнейшего пространственного анализа.

Перед началом работы важно проверить:

указана ли CRS в данных;
какой тип CRS используется;
подходят ли единицы измерения для анализа;
совпадает ли CRS у всех используемых слоёв.

Проверим систему координат данных:

admin_border.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

Наши данные находятся в географической системе координат WGS 84 (EPSG:4326).
Координаты объектов записаны в градусах широты и долготы.

2. Определение подходящей UTM-зоны¶

При работе с городами или локальными территориями одной из наиболее часто используемых систем координат является UTM.

Как мы уже упоминали в предыдущем разделе, данные в географической системе координат измеряются в градусах, а расстояние, соответствующее одному градусу долготы, меняется в зависимости от широты. Поэтому такие координаты не подходят для точных измерений расстояний и площадей.

Система координат UTM (Universal Transverse Mercator) использует метры. Это позволяет значительно точнее рассчитывать расстояния, площади и направления при анализе данных на локальном уровне.

Теперь разберёмся, как определить подходящую UTM-зону для данных.

2.1. С помощью формулы¶

UTM-зону можно определить самостоятельно на основе координат точки.

Чтобы найти номер зоны, достаточно знать долготу.
Широта используется для определения того, находится ли точка в северном или южном полушарии.

Номер UTM-зоны можно вычислить по следующей формуле:

\text{zone} = \left\lfloor \frac{\text{longitude} + 180}{6} \right\rfloor + 1

(1)

Сначала долгота сдвигается из диапазона [-180, +180] в [0, 360], , затем результат делится на 6, поскольку каждая UTM-зона имеет ширину 6° долготы. После этого значение округляется вниз, чтобы определить номер зоны, а затем прибавляется +1, так как нумерация UTM-зон начинается с 1, а не с 0.

После определения номера зоны можно получить соответствующий EPSG-код:

EPSG:326xx — для зон северного полушария
EPSG:327xx — для зон южного полушария

где xx — номер UTM-зоны.

Напишем небольшую функцию, которая по долготе и широте точки определяет номер UTM-зоны и возвращает соответствующий EPSG-код.

def utm_epsg(longitude, latitude):
    zone = int((longitude + 180) / 6) + 1
    
    if latitude >= 0:
        epsg = 32600 + zone   # северное полушарие
    else:
        epsg = 32700 + zone   # южное полушарие
        
    return epsg

Рассмотрим пример для координат Санкт-Петербурга:

lon = 30.3
lat = 59.9

utm_epsg(lon, lat)

32636

Этот EPSG код можно сразу использовать для перепроецирования данных на следующих шагах.

2.2. Автоматическое определение¶

На практике вычислять UTM-зону вручную обычно не требуется. GeoPandas может автоматически определить подходящую систему координат UTM с помощью метода .estimate_utm_crs().

Этот метод анализирует географическое положение данных и определяет UTM-зону на основе центра геометрии слоя.

data_crs = admin_border.estimate_utm_crs()

data_crs

<Projected CRS: EPSG:32636>
Name: WGS 84 / UTM zone 36N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 30°E and 36°E, northern hemisphere between equator and 84°N, onshore and offshore. Belarus. Cyprus. Egypt. Ethiopia. Finland. Israel. Jordan. Kenya. Lebanon. Moldova. Norway. Russian Federation. Saudi Arabia. Sudan. Syria. Türkiye (Turkey). Uganda. Ukraine.
- bounds: (30.0, 0.0, 36.0, 84.0)
Coordinate Operation:
- name: UTM zone 36N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Как видно, результат совпадает с EPSG-кодом, полученным ранее с помощью формулы. Его также можно использовать для перепроецирования данных.

3. Перепроецирование¶

Перепроецирование — это преобразование пространственных данных из одной системы координат в другую.

На практике это означает пересчёт координат каждой точки, чтобы они соответствовали другой картографической проекции.

Давайте посмотрим, как записаны координаты геометрий в поле geometry.

admin_border.geometry

0    POLYGON ((30.30827 59.94106, 30.30908 59.94051...
Name: geometry, dtype: geometry

Мы видим, что каждая вершина полигона представлена парой координат. В географической системе координат EPSG:4326 эти значения соответствуют долготе и широте, выраженным в градусах.

Для перепроецирования данных в GeoPandas используется метод .to_crs().

Этот метод пересчитывает координаты всех геометрий, чтобы данные корректно отображались в новой системе координат.

Теперь перепроецируем данные в подходящую UTM-зону, которую мы ранее определили с помощью метода estimate_utm_crs() и записали в переменную data_crs.

admin_border_utm = admin_border.to_crs(data_crs)

Сначала убедимся, что система координат действительно изменилась.

admin_border_utm.crs

<Projected CRS: EPSG:32636>
Name: WGS 84 / UTM zone 36N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 30°E and 36°E, northern hemisphere between equator and 84°N, onshore and offshore. Belarus. Cyprus. Egypt. Ethiopia. Finland. Israel. Jordan. Kenya. Lebanon. Moldova. Norway. Russian Federation. Saudi Arabia. Sudan. Syria. Türkiye (Turkey). Uganda. Ukraine.
- bounds: (30.0, 0.0, 36.0, 84.0)
Coordinate Operation:
- name: UTM zone 36N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Теперь данные находятся в проецированной системе координат UTM, где координаты выражены в метрах.

Посмотрим, как теперь записаны координаты геометрии.

admin_border_utm.geometry

0    POLYGON ((349622.326 6647904.934, 349664.793 6...
Name: geometry, dtype: geometry

В отличие от предыдущего вывода, координаты теперь представлены в метрах, а не в градусах. Это можно заметить по порядку величины чисел: их значения становятся значительно больше — обычно сотни тысяч или миллионы.

Мы успешно перепроецировали слой в подходящую UTM-зону. Это означает, что данные теперь находятся в системе координат, которая подходит для точных пространственных вычислений — например, измерения расстояний, расчёта площадей и выполнения других аналитических операций.

4. Назначение системы координат (CRS)¶

Иногда пространственные данные не содержат информации о системе координат (CRS). В таком случае GeoPandas не может корректно интерпретировать координаты геометрий, поскольку неизвестно, в какой системе координат они заданы

4.1 Для данных без указанной системы координат¶

Иногда встречаются пространственные данные, у которых не указана система координат (CRS). В этом случае необходимо явно задать CRS, чтобы корректно интерпретировать координаты.

Рассмотрим пример с CSV-файлом, содержащим координаты театров в Санкт-Петербурге.

В этом разделе будут использованы файлы из нашего репозитория:

spb_theaters.csv — данные о театрах Санкт-Петербурга.
Источник: Портал открытых данных Санкт-Петербурга

Сначала прочитаем данные и создадим из них GeoDataFrame:

df_csv = pd.read_csv('./data/spb_theaters.csv')

gdf_csv = gpd.GeoDataFrame(
    df_csv,
    geometry=gpd.points_from_xy(df_csv['longitude'], df_csv['latitude'])
)

Проверим, какая система координат указана у gdf_csv:

print(gdf_csv.crs)

None

Результат — None, потому что при создании GeoDataFrame мы не указали систему координат.

На практике это означает, что GeoPandas не знает, как интерпретировать координаты геометрий.

В нашем примере координаты представлены как долгота и широта, поэтому можно предположить, что они заданы в географической системе координат WGS 84 (EPSG:4326).

Если CRS не задана, её можно добавить двумя способами:

Способ 1: использовать метод .set_crs()
Способ 2: напрямую присвоить значение через атрибут .crs

(Однако важно помнить, что стоит CRS сразу при создании GeoDataFrame. Это помогает избежать ошибок при дальнейшей работе с пространственными данными.)

poi_gdf = gdf_csv.set_crs(epsg=4326) # option 1

poi_gdf.crs = "EPSG:4326" # option 2

print(poi_gdf.crs)

EPSG:4326

Назначение CRS не изменяет координаты, а только сообщает GeoPandas, в какой системе координат они записаны

Самое важное различие при работе с CRS
Перепроецирование (to_crs)
Преобразование координат из одной системы координат (CRS) в другую.
При этом изменяется сама CRS и значения координат всех геометрий.
Назначение CRS (set_crs)
Указание системы координат без изменения координат.
Мы просто сообщаем GeoPandas, как интерпретировать существующие координаты.

Итог¶

В этом разделе мы рассмотрели, как перепроецировать пространственные данные или задавать систему координат, если она отсутствует.

Мы научились:

определять подходящую UTM-зону для набора данных с помощью метода;
использовать .to_crs() для перепроецирования данных;
различать назначение CRS (указание системы координат без изменения координат) и перепроецирование (фактическое преобразование координат).

Корректная работа с системами координат является важным этапом подготовки пространственных данных