Земля не плоская, и ваша диаграмма Вороного не должна быть такой же

Мир не плоский, а ваша диаграмма Вороного не должна быть такой же

История о точности, раскрывающая силу сферических геопространственных диаграмм Вороного с использованием Python

Земля со сферической диаграммой Вороного, движущейся между 2 проекциями: ортогональной и эквиректангулярной. Сгенерировано автором с использованием библиотеки D3.js.

Возможно, вы уже знакомы с диаграммами Вороного и их использованием в геопространственном анализе. Если нет, вот краткое объяснение: они делят плоскость на регионы, состоящие из всех точек плоскости, ближайших к данному зерну, чем к любому другому. Они названы в честь математика Георгия Вороного. Больше информации можно найти на Википедии.

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

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

Размеры и проекции – почему это важно?

Если мы хотим увидеть данные на карте, нам необходимо работать с проекциями. Чтобы показать что-то на 2D плоскости, мы должны проецировать координаты из 3D координат на глобусе. Самая популярная проекция, которую мы все знаем и используем – проекция Меркатора (Веб-Меркатор или точнее WGS84 Меркатор, так как им пользуются большинство поставщиков карт), а самая популярная система координат – Мировая геодезическая система 1984 – WGS84 (или EPSG:4326). Эта система основана на градусах, и долгота в ней варьируется от -180° до 180° (от запада к востоку), а широта от -90° до 90° (с юга на север).

Каждая проекция на 2D плоскость имеет некоторые искажения. Проекция Меркатора является конформной картографической проекцией, что означает, что углы между объектами на Земле должны сохраняться. Чем выше выше 0° (или ниже 0°) широта, тем больше искажение в области и расстоянии. Поскольку диаграмма Вороного сильно зависит от расстояния между зернами, то же искажение передается при ее генерации.

Земля имеет неправильную эллипсоидную форму, но для наших целей она может быть приближена сферой. Генерируя диаграмму Вороного на сфере, мы можем правильно вычислить расстояние на основе дуг на поверхности сферы. Позже мы можем сопоставить сгенерированные сферические полигоны с проецированными 2D координатами, и мы можем быть уверены, что линия, разделяющая две смежные ячейки Вороного, будет перпендикулярна линии, соединяющей два зерна, определяющие эти ячейки.

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

Пример разницы между углами и расстояниями в обоих версиях диаграммы Вороного. Изображение автора.

Еще одна проблема состоит в том, что вы не можете сравнивать регионы в разных частях мира (т.е. не находящихся на одной широте), если вы используете 2D диаграмму Вороного, поскольку области будут сильно искажены.

Полный Jupyter-блокнот с использованным кодом в примерах ниже можно найти на GitHub. Здесь некоторые функции пропущены для краткости.

Необходимые условия

Установите необходимые библиотеки

pip install -q srai[voronoi,osm,plotting] geodatasets

Импортировать необходимые модули и функции

import geodatasetsimport geopandas as gpdimport matplotlib.pyplot as pltimport plotly.express as pxfrom shapely.geometry import MultiPoint, Pointfrom shapely.ops import voronoi_diagramfrom srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf

Первый пример

Определим шесть точек на земном шаре: Северный и Южный полюса, и четыре точки на экваторе.

earth_points_gdf = gpd.GeoDataFrame(    geometry=[        Point(0, 0),        Point(90, 0),        Point(180, 0),        Point(-90, 0),        Point(0, 90),        Point(0, -90),    ],    index=[1, 2, 3, 4, 5, 6],    crs="EPSG:4326",)
Изображение автора.

Создание диаграммы Вороного с использованием функции voronoi_diagram из библиотеки Shapely

def generate_flat_voronoi_diagram_regions(    seeds_gdf: gpd.GeoDataFrame,) -> gpd.GeoDataFrame:    points = MultiPoint(seeds_gdf.geometry.values)    # Генерация 2D диаграммы    regions = voronoi_diagram(points)    # Преобразование геометрий в GeoDataFrame    flat_voronoi_regions = gpd.GeoDataFrame(        geometry=list(regions.geoms),        crs="EPSG:4326",    )    # Применение индексов из исходного GeoDataFrame    flat_voronoi_regions.index = gpd.pd.Index(        flat_voronoi_regions.sjoin(seeds_gdf)["index_right"],        name="region_id",    )    # Обрезка границ земли    flat_voronoi_regions.geometry = flat_voronoi_regions.geometry.clip_by_rect(        xmin=-180, ymin=-90, xmax=180, ymax=90    )    return flat_voronoi_regions

earth_poles_flat_voronoi_regions = generate_flat_voronoi_diagram_regions(    earth_points_gdf)

Создание диаграмм Вороного с использованием VoronoiRegionalizer из библиотеки srai. Внутри используется реализация SphericalVoronoi из библиотеки scipy с правильным преобразованием координат WGS84 в сферическую систему координат и обратно.

earth_points_spherical_voronoi_regions = VoronoiRegionalizer(    seeds=earth_points_gdf).transform()

Давайте посмотрим на разницу между этими двумя диаграммами на графиках.

Разница между диаграммами Вороного в плоской (слева) и сферической (справа) версиях в координатах WGS84. Создано автором с использованием библиотеки GeoPandas.

Разница между диаграммами Вороного в плоском (слева) и сферическом (справа) вариантах в ортогональной проекции. Создано автором с использованием Plotly.

Первое, что можно увидеть, это то, что двумерная диаграмма Вороного не проходит по всей поверхности Земли, поскольку она работает на плоской картезианской плоскости. Сферическая диаграмма Вороного правильно охватывает Землю и не обрывается на антимеридиане (где долгота переключается с 180° на -180°).

Чтобы численно измерить разницу, мы можем вычислить метрику IoU (Intersection over Union) (или индекс Жаккара), чтобы измерить разницу между формами полигонов. Значение этой метрики находится между 0 и 1, где 0 означает отсутствие перекрытия, а 1 означает полное перекрытие.

def calculate_iou(    flat_regions: gpd.GeoDataFrame, spherical_regions: gpd.GeoDataFrame) -> float:    total_intersections_area = 0    total_unions_area = 0    # Перебор всех регионов    for index in spherical_regions.index:        # Нахождение совпадающего сферического и плоского региона Вороного        spherical_region_geometry = spherical_regions.loc[index].geometry        flat_region_geometry = flat_regions.loc[index].geometry        # Вычисление площади их пересечения        intersections_area = spherical_region_geometry.intersection(            flat_region_geometry        ).area        # Вычисление площади их объединения        # Альтернативный код:        # spherical_region_geometry.union(flat_region_geometry).area        unions_area = (            spherical_region_geometry.area            + flat_region_geometry.area            - intersections_area        )        # Добавление к общим суммам        total_intersections_area += intersections_area        total_unions_area += unions_area    # Деление площади пересечения на площадь объединения    return round(total_intersections_area / total_unions_area, 3)

calculate_iou(    earth_points_flat_voronoi_regions, earth_points_spherical_voronoi_regions)

Вычисленное значение составляет 0,423, что достаточно низко и на большой шкале эти полигоны отличаются друг от друга, что легко видно на графиках выше.

Пример реальных данных: разделение земного шара с помощью позиций AED (автоматических внешних дефибрилляторов)

Используемые данные в этом примере взяты с OpenAEDMap и основаны на данных OpenStreetMap. Подготовленный файл содержит отфильтрованные позиции (в точности 80694), не имеющие дублирующих узлов, определенных друг на друге.

# Загрузка позиций AED в GeoDataFrameaed_world_gdf = gpd.read_file(    "https://raw.githubusercontent.com/RaczeQ/medium-articles/main/articles/spherical-geovoronoi/aed_world.geojson")

Генерация диаграммы Вороного для позиций AED

aed_flat_voronoi_regions = generate_flat_voronoi_diagram_regions(aed_world_gdf)aed_spherical_voronoi_regions = VoronoiRegionalizer(    seeds=aed_world_gdf, max_meters_between_points=1_000).transform()

Давайте сравним эти диаграммы Вороного.

Разница между диаграммами Вороного в плоском (слева) и сферическом (справа) вариантах в WGS84 координатах. Создано автором с использованием GeoPandas.

Разница между диаграммами Вороного в плоской (слева) и сферической (справа) версиями в ортогональной проекции. Создано автором с использованием Plotly.

Разница довольно очевидна при просмотре графиков. Все границы в 2D-версии прямые, тогда как сферические границы выглядят довольно изогнутыми в координатах WGS84. Вы также можете ясно видеть, что в плоской версии много регионов сходится на полюсах (ортогональная проекция фокусируется на южном полюсе), тогда как сферическая версия этого не делает. Еще одна заметная разница – это непрерывность вокруг антимеридиана, которая была упомянута в первом примере. Регионы, появляющиеся из Новой Зеландии, резко обрываются в плоской версии.

Давайте посмотрим на значение IoU:

calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)

Вычисленное значение составляет 0.511, что немного лучше, чем в первом примере, но все же полигоны соответствуют примерно 50%.

Увеличение до уровня города

Давайте посмотрим разницу на меньшем масштабе. Мы можем выбрать все ПНИ, находящиеся в Лондоне, и построить их.

greater_london_area = geocode_to_region_gdf("Greater London")aeds_in_london = aed_world_gdf.sjoin(greater_london_area)
2D- и сферические диаграммы Вороного, наложенные друг на друга красным и синим цветами. Изображение автором.
calculate_iou(    aed_flat_voronoi_regions.loc[aeds_in_london.index],    aed_spherical_voronoi_regions.loc[aeds_in_london.index])

Значение равно 0.675. Оно становится лучше, но все еще есть заметная разница. Поскольку ПНИ располагаются плотнее, формы и расстояния становятся меньше, поэтому различия между диаграммами Вороного, рассчитанными на проекционной 2D-плоскости и на сфере, уменьшаются.

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

Изображение автором.

Площади этих полигонов в основном совпадают, но можно заметить различия в углах и формах. Эти расхождения могут быть важными в пространственном анализе и могут изменить результаты анализа. Чем больше интересующая область, тем больше становится разница.

Итог

Я надеюсь, теперь вы видите, почему сферическая диаграмма Вороного более подходит для использования в геопространственной области, чем плоская.

Большая часть анализов в этой области в настоящее время выполняется с использованием диаграмм Вороного на проекционной 2D-плоскости, что может привести к неверным результатам.

Долгое время не было простого решения для сферических диаграмм Вороного, доступных для геопространственных ученых и аналитиков данных, работающих в Python. Теперь это так же просто, как установка одной библиотеки. Конечно, расчет занимает немного больше времени, чем плоское решение, поскольку он должен проецировать точки на сферу и обратно из сферических координат, при этом правильно обрезая полигоны, пересекающие антимеридиан, но это не должно иметь значения, если вы хотите сохранить точность ваших анализов. Для пользователей JavaScript уже доступна реализация сферической диаграммы Вороного в D3.js.

Отказ от ответственности

Я являюсь одним из сотрудников библиотеки srai.