Una historia sobre precisión, que revela el poder de los diagramas esféricos geoespaciales de Voronoi con Python
Quizás esté familiarizado con los diagramas de Voronoi y sus usos en los análisis geoespaciales. Si no, aquí está el TL;DR rápido: divide el plano en regiones que consisten en todos los puntos del plano más cercanos a una semilla determinada que a cualquier otra. Lleva el nombre del matemático Georgy Voronoy. Puedes leer más sobre esto en el Wikipedia.
¿Cómo se aplica al dominio geoespacial? Usando los diagramas de Voronoi puedes encontrar rápidamente la parada de transporte público más cercana para los habitantes de una ciudad determinada a mayor escala, más rápido que calcularla individualmente para cada edificio por separado. O también puedes utilizarlo por ejemplo en el análisis de cuota de mercado entre diferentes marcas.
En este post quiero mostrar las diferencias entre el típico diagrama de Voronoi calculado con coordenadas proyectadas en un plano y el esférico, y con suerte, mostrar la superioridad de este último.
Dimensiones y proyecciones: ¿por qué es importante?
Si queremos ver datos en el mapa, tenemos que trabajar con proyecciones. Para mostrar algo en el plano 2D, tenemos que proyectar las coordenadas a partir de las coordenadas 3D en el globo.
La proyección más popular que todos conocemos y utilizamos es la proyección de Mercator (Mercator Web o WGS84 Mercator para ser precisos, ya que es utilizado 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 oscila en longitud de -180° a 180° (de oeste a este) y en latitud de -90° a 90° (de sur a norte).
Cada proyección al plano 2D tiene algunas distorsiones. El Mercator es un conforme Proyección cartográfica, lo que significa que se deben preservar los ángulos entre los objetos de la Tierra. Cuanto mayor sea la latitud por encima de 0° (o menor por debajo de 0°), mayor será la distorsión en el área y el distancia. Debido a que el diagrama de Voronoi depende en gran medida de la distancia entre las semillas, se transmite el mismo error de distorsión al generar el diagrama.
La Tierra es un elipsoide de forma irregular, pero para nuestros propósitos, se puede aproximar a ella por su forma esférica. Al generar el diagrama de Voronoi en la esfera, podemos calcular correctamente la distancia en función de los arcos en la superficie de una esfera. Posteriormente, podemos asignar los polígonos esféricos generados a las coordenadas 2D proyectadas y podemos estar seguros de que la línea que separa dos celdas de Voronoi adyacentes será perpendicular a la línea que conecta las dos semillas que definen estas celdas.
A continuación puedes ver el problema de ángulos y distancias que describí anteriormente. Aunque las líneas se cruzan en el mismo punto, las formas y ángulos de las células de Voronoi difieren.
Otro problema es que no se pueden comparar regiones en diferentes partes del mundo (es decir, que no se encuentran 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 siguientes se puede encontrar en GitHub. Aquí se omiten algunas funciones por motivos de brevedad.
Requisitos previos
Instalar las bibliotecas necesarias
pip install -q srai[voronoi,osm,plotting] geodatasets
Importar módulos y funciones necesarios
import geodatasets
import geopandas as gpd
import matplotlib.pyplot as plt
import plotly.express as px
from shapely.geometry import MultiPoint, Point
from shapely.ops import voronoi_diagram
from srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf
Primer ejemplo
Definamos seis puntos del 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 usando diagrama_voronoi de la biblioteca Shapely
def generate_flat_voronoi_diagram_regions(
seeds_gdf: gpd.GeoDataFrame,
) -> gpd.GeoDataFrame:
points = MultiPoint(seeds_gdf.geometry.values)
# Generate 2D diagram
regions = voronoi_diagram(points)
# Map geometries to GeoDataFrame
flat_voronoi_regions = gpd.GeoDataFrame(
geometry=list(regions.geoms),
crs="EPSG:4326",
)
# Apply indexes from the seeds dataframe
flat_voronoi_regions.index = gpd.pd.Index(
flat_voronoi_regions.sjoin(seeds_gdf)["index_right"],
name="region_id",
)
# Clip to Earth boundaries
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
)
Genere diagramas de Voronoi usando Voronoi Regionalizador de la biblioteca srai.
Debajo del capó, utiliza el EsféricoVoronoi implementación desde 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 las tramas.
Lo primero que se puede ver es que el diagrama de Voronoi 2D no da la vuelta al mundo, ya que funciona en un plano. plano cartesiano. El diagrama esférico de Voronoi cubre adecuadamente la Tierra y no se rompe en el anti meridiano línea (donde la longitud cambia de 180° a -180°).
Para cuantificar numéricamente la diferencia podemos calcular el pagaré (Intersección sobre Unión) métrica (o Índice Jaccard) para medir la diferencia entre las formas de los polígonos. El valor de esta métrica está entre 0 y 1, donde 0 significa que no hay superposición y 1 significa superposición total.
def calculate_iou(
flat_regions: gpd.GeoDataFrame, spherical_regions: gpd.GeoDataFrame
) -> float:
total_intersections_area = 0
total_unions_area = 0
# Iterate all regions
for index in spherical_regions.index:
# Find matching spherical and flat Voronoi region
spherical_region_geometry = spherical_regions.loc[index].geometry
flat_region_geometry = flat_regions.loc[index].geometry
# Calculate their intersection area
intersections_area = spherical_region_geometry.intersection(
flat_region_geometry
).area
# Calculate their union area
# Alternative code:
# spherical_region_geometry.union(flat_region_geometry).area
unions_area = (
spherical_region_geometry.area
+ flat_region_geometry.area
- intersections_area
)
# Add to the total sums
total_intersections_area += intersections_area
total_unions_area += unions_area
# Divide the intersection area by the union area
return round(total_intersections_area / total_unions_area, 3)
calculate_iou(
earth_points_flat_voronoi_regions, earth_points_spherical_voronoi_regions
)
El valor calculado es 0.423que es bastante bajo y, a gran escala, esos polígonos son diferentes entre sí, lo que se puede ver fácilmente en los gráficos de arriba.
Ejemplo de datos reales: dividir el globo usando posiciones de DEA (Desfibriladores Externos Automáticos)
Los datos utilizados en este ejemplo provienen de la AbrirAEDMapa y se basa en Abrir mapa de calles datos. El archivo preparado tiene posiciones filtradas (80694 para ser exactos) sin nodos duplicados definidos uno encima del otro.
# Load AEDs positions to GeoDataFrame
aed_world_gdf = gpd.read_file(
"https://raw.githubusercontent.com/RaczeQ/medium-articles/main/articles/spherical-geovoronoi/aed_world.geojson"
)
Generar diagramas de Voronoi para los DEA
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()
Comparemos estos diagramas de Voronoi.
La diferencia es bastante obvia al observar las tramas. 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 centra en el polo sur), mientras que la esférica no. Otra diferencia visible es la continuidad alrededor del antimeridiano, que se mencionó en el primer ejemplo. En la versión plana, las regiones que emergen de Nueva Zelanda están abruptamente recortadas.
Veamos el valor de IoU:
calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)
El valor calculado es 0.511que es ligeramente mejor que el primer ejemplo, pero aún así, los polígonos coinciden aproximadamente en un 50%.
Acercándonos a la escala de la ciudad
Veamos la diferencia a menor escala. Podemos seleccionar todos los DEA que se encuentran en Londres y trazarlo.
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 todavía hay una diferencia notable. Dado que los DEA se colocan más densos, las formas y distancias se hacen 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 uno encima del otro.
Las áreas de esos polígonos coinciden en su mayoría, pero puedes ver las diferencias en ángulos y formas. Esas discrepancias podrían ser importantes en el análisis espacial y podrían cambiar los resultados de los mismos. Cuanto mayor sea el área de interés, mayor será la diferencia.
Resumen
Espero que ahora puedan ver por qué el diagrama esférico de Voronoi 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 erróneos.
Durante mucho tiempo, no hubo una solución sencilla para los diagramas esféricos de Voronoi disponibles para los científicos y analistas de datos geoespaciales que trabajaban en Python. Ahora es tan fácil como instalar una biblioteca.
Claro, calcula un poco más que la solución plana, ya que tiene que proyectar puntos hacia y desde coordenadas esféricas, mientras recorta adecuadamente los polígonos que se cruzan con el antimeridiano, pero no debería importar si desea preservar la precisión en sus análisis.
Para los usuarios de JavaScript, ya existe un Voronoi esférico disponible Implementación de D3.js.
Descargo de responsabilidad
Soy uno de los mantenedores de la biblioteca srai.
La Tierra no es plana y sus diagramas de Voronoi tampoco deberían serlo fue publicado originalmente en Hacia la ciencia de datos en Medium, donde las personas continúan la conversación resaltando y respondiendo a esta historia.