La Tierra no es plana, y tus diagramas de Voronoi tampoco deberían serlo

La Tierra no es plana, y tus diagramas de Voronoi tampoco deberían serlo.

Una historia sobre precisión, revelando el poder de los diagramas Voronoi geoespaciales esféricos con Python

Tal vez estés familiarizado con los diagramas Voronoi y su uso en los análisis geoespaciales. Si no es así, aquí tienes un resumen rápido: divide el plano en regiones que consisten en todos los puntos del plano que están más cerca de una semilla dada que de cualquier otra. Recibe su nombre del matemático Georgy Voronoy. Puedes leer más al respecto en la página de Wikipedia.

¿Cómo se aplica al ámbito geoespacial? Usando diagramas Voronoi, puedes encontrar rápidamente la parada de transporte público más cercana para los habitantes de una ciudad dada a una escala mayor, más rápido que calcularlo individualmente para cada edificio por separado. O también puedes usarlo, por ejemplo, en el análisis de participación de mercado entre diferentes marcas.

En esta publicación, quiero mostrar las diferencias entre el típico diagrama Voronoi calculado con coordenadas proyectadas en un plano plano y el diagrama esférico, y espero mostrar la superioridad de este último.

Dimensiones y proyecciones: ¿por qué importa?

Si queremos ver datos en el mapa, tenemos que trabajar con proyecciones. Para mostrar algo en el plano 2D, tenemos que proyectar las coordenadas desde las coordenadas 3D en el globo. La proyección más popular que todos conocemos y usamos es la proyección de Mercator (Mercator Web o Mercator WGS84 para ser precisos, ya que es utilizada por la mayoría de los proveedores de mapas) y el sistema de coordenadas más popular es el Sistema Geodésico Mundial 1984 – WGS84 (o EPSG:4326). Este sistema se basa en grados y varía en longitud de -180° a 180° (de oeste a este) y en latitud de -90° a 90° (de sur a norte).

• Segmentar cualquier cosa en 3D para nubes de puntos Guía completa (SAM 3D)
• Crea una interfaz de usuario web para interactuar con LLMs utilizando Amazon SageMaker JumpStart
• Conozca a NANA, el avatar de recepcionista con inteligencia artificial de Moonshine Studio

Cada proyección a un plano 2D tiene algunas distorsiones. La proyección de Mercator es una proyección cartográfica conforme, lo que significa que los ángulos entre objetos en la Tierra deben conservarse. Cuanto más arriba de 0° (o más abajo de 0°) esté la latitud, mayor será la distorsión en el área y la distancia. Debido a que el diagrama Voronoi depende en gran medida de la distancia entre las semillas, el mismo error de distorsión se transmite al generar el diagrama.

La Tierra es una elipsoide con forma irregular, pero para nuestros propósitos, se puede aproximar con la forma de una esfera. Al generar el diagrama de Voronoi en la esfera, podemos calcular correctamente la distancia basada en los arcos en la superficie de la esfera. Luego, podemos mapear los polígonos esféricos generados a las coordenadas proyectadas en 2D y asegurarnos de que la línea que separa dos celdas de Voronoi adyacentes sea perpendicular a la línea que conecta las dos semillas que definen estas celdas.

A continuación, puedes ver el problema de los ángulos y distancias que he descrito anteriormente. Aunque las líneas se cruzan en el mismo punto, las formas y los ángulos de las celdas de Voronoi difieren.

Otro problema es que no se pueden comparar las regiones en diferentes partes del mundo (es decir, que no se encuentren en la misma latitud) si se utiliza un diagrama de Voronoi 2D, ya que las áreas estarán muy distorsionadas.

El cuaderno completo de Jupyter con el código utilizado en los ejemplos se puede encontrar en GitHub. Aquí se omiten algunas funciones por brevedad.

Prerrequisitos

Instalar las bibliotecas requeridas

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

Importar los módulos y funciones requeridos

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

Primer ejemplo

Definamos seis puntos en el globo: los polos norte y sur, y cuatro puntos en el ecuador.

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",)

Generar diagrama de Voronoi utilizando voronoi_diagram de la biblioteca Shapely

def generate_flat_voronoi_diagram_regions(    seeds_gdf: gpd.GeoDataFrame,) -> gpd.GeoDataFrame:    points = MultiPoint(seeds_gdf.geometry.values)    # Generar diagrama 2D    regions = voronoi_diagram(points)    # Mapear geometrías a GeoDataFrame    flat_voronoi_regions = gpd.GeoDataFrame(        geometry=list(regions.geoms),        crs="EPSG:4326",    )    # Aplicar índices del GeoDataFrame inicial    flat_voronoi_regions.index = gpd.pd.Index(        flat_voronoi_regions.sjoin(seeds_gdf)["index_right"],        name="region_id",    )    # Recortar a los límites terrestres    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)

Generar diagramas de Voronoi utilizando VoronoiRegionalizer de la biblioteca srai. Bajo la capa, utiliza la implementación SphericalVoronoi de la biblioteca scipy y transforma adecuadamente las coordenadas WGS84 hacia y desde el sistema de coordenadas esféricas.

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

Veamos la diferencia entre los dos en los gráficos.

Lo primero que se puede observar es que el diagrama de Voronoi 2D no se extiende alrededor del globo, ya que funciona en un plano cartesiano plano. El diagrama de Voronoi esférico cubre correctamente la Tierra y no se ve afectado por la línea del antimeridiano (donde la longitud cambia de 180° a -180°).

Para cuantificar numéricamente la diferencia, podemos calcular la métrica de IoU (Intersección sobre Unión) (o Índice de Jaccard) para medir la diferencia entre las formas de los polígonos. El valor de esta métrica varía entre 0 y 1, donde 0 significa que no hay superposición y 1 significa una superposición completa.

def calcular_iou(    regiones_planas: gpd.GeoDataFrame, regiones_esféricas: gpd.GeoDataFrame) -> float:    area_total_intersecciones = 0    area_total_union = 0    # Iterar todas las regiones    for indice in regiones_esféricas.index:        # Encontrar la correspondencia entre la región de Voronoi esférica y plana        geometría_region_esférica = regiones_esféricas.loc[indice].geometría        geometría_region_plana = regiones_planas.loc[indice].geometría        # Calcular el área de intersección        area_intersecciones = geometría_region_esférica.intersection(            geometría_region_plana        ).área        # Calcular el área de unión        # Código alternativo:        # geometría_region_esférica.union(geometría_region_plana).área        area_union = (            geometría_region_esférica.área            + geometría_region_plana.área            - area_intersecciones        )        # Agregar a las sumas totales        area_total_intersecciones += area_intersecciones        area_total_union += area_union    # Dividir el área de intersección por el área de unión    return round(area_total_intersecciones / area_total_union, 3)

calcular_iou(    regiones_voronoi_planas_terrestres, regiones_voronoi_esféricas_terrestres)

El valor calculado es 0.423, lo cual es bastante bajo, lo que indica que estos polígonos son diferentes entre sí, como se puede ver fácilmente en los gráficos anteriores.

Ejemplo de datos reales: dividir el globo utilizando posiciones de AED (Desfibriladores Externos Automatizados)

Los datos utilizados en este ejemplo provienen de OpenAEDMap y se basan en datos de OpenStreetMap. El archivo preparado tiene posiciones filtradas (80694 para ser exactos) sin nodos duplicados definidos uno encima del otro.

# Cargar posiciones de AED en GeoDataFrameaed_mundial_gdf = gpd.read_file(    "https://raw.githubusercontent.com/RaczeQ/medium-articles/main/articles/spherical-geovoronoi/aed_world.geojson")

Generar diagramas de Voronoi para los AED

regiones_voronoi_planas_aed = generar_regiones_diagrama_voronoi_plano(aed_mundial_gdf)regiones_voronoi_esféricas_aed = VoronoiRegionalizer(    semillas=aed_mundial_gdf, metros_máximos_entre_puntos=1_000).transform()

Comparemos estos diagramas de Voronoi.

La diferencia es bastante obvia al observar las representaciones gráficas. Todos los bordes en la versión 2D son rectos, mientras que los esféricos parecen bastante curvados en las coordenadas WGS84. También se puede ver claramente que en la versión plana, muchas regiones convergen en los polos (la proyección ortogonal se enfoca en el polo sur), mientras que en la versión esférica no sucede esto. Otra diferencia visible es la continuidad alrededor del anti-meridiano, que se mencionó en el primer ejemplo. Las regiones que emergen de Nueva Zelanda se cortan abruptamente en la versión plana.

Veamos el valor de IoU:

calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)

El valor calculado es 0.511, que es ligeramente mejor que el primer ejemplo, pero aún así, los polígonos coinciden aproximadamente en un 50%.

Acercamiento a la escala de la ciudad

Veamos la diferencia en una escala más pequeña. Podemos seleccionar todos los AED que se encuentren en Londres y representarlos.

greater_london_area = geocode_to_region_gdf("Greater London")aeds_in_london = aed_world_gdf.sjoin(greater_london_area)

calculate_iou(    aed_flat_voronoi_regions.loc[aeds_in_london.index],    aed_spherical_voronoi_regions.loc[aeds_in_london.index],)

El valor es 0.675. Está mejorando, pero aún es una diferencia notable. Dado que los AED están ubicados más densamente, las formas y distancias se vuelven más pequeñas, por lo que las diferencias entre los diagramas de Voronoi calculados en el plano 2D proyectado y en una esfera disminuyen.

Veamos algunos ejemplos individuales superpuestos.

Las áreas de esos polígonos coinciden en su mayoría, pero se pueden ver las diferencias en ángulos y formas. Esas discrepancias podrían ser importantes en el análisis espacial y podrían cambiar los resultados. A medida que aumenta el área de interés, la diferencia se vuelve más grande.

Resumen

Espero que ahora puedas ver por qué el diagrama de Voronoi esférico es más adecuado para su uso en el dominio geoespacial que el plano.

La mayoría de los análisis en el dominio se realizan actualmente utilizando diagramas de Voronoi en un plano 2D proyectado, lo que podría conducir a resultados incorrectos.

Durante mucho tiempo, no hubo una solución sencilla para los diagramas de Voronoi esféricos disponibles para científicos de datos y analistas geoespaciales que trabajan en Python. Ahora es tan fácil como instalar una biblioteca. Por supuesto, se calcula un poco más lento que la solución plana, ya que debe proyectar puntos hacia y desde coordenadas esféricas, mientras recorta correctamente los polígonos que intersecan el anti-meridiano, pero no importa si quieres preservar la precisión en tus análisis. Para los usuarios de JavaScript, ya está disponible una implementación de diagrama de Voronoi esférico con D3.js.

Aviso legal

Soy uno de los responsables de la biblioteca srai.

We will continue to update Zepes; if you have any questions or suggestions, please contact us!