La Terra non è piatta, e nemmeno dovrebbero esserlo i tuoi diagrammi di Voronoi

Il Mondo non è piatto, e nemmeno dovrebbero esserlo i tuoi diagrammi di Voronoi

Un racconto sulla precisione, svelando il potere dei diagrammi di Voronoi geospaziali sferici con Python

Potresti essere familiare con i diagrammi di Voronoi e il loro utilizzo nelle analisi geospaziali. In caso contrario, ecco il rapido riassunto: divide il piano in regioni composte da tutti i punti del piano più vicini a un dato seme rispetto a qualsiasi altro. Prende il nome dal matematico Georgy Voronoy. Puoi leggere di più su di esso su Wikipedia.

Come si applica al campo geospaziale? Utilizzando i diagrammi di Voronoi, puoi trovare rapidamente la fermata dei mezzi pubblici più vicina per gli abitanti di una data città su una scala più ampia, più velocemente che calcolarla singolarmente per ogni edificio separatamente. O puoi anche utilizzarlo ad esempio nell’analisi della quota di mercato tra diverse marche.

In questo articolo, voglio mostrare le differenze tra il tipico diagramma di Voronoi calcolato con coordinate proiettate su un piano piatto e quello sferico e spero di dimostrare la superiorità di quest’ultimo.

Dimensioni e proiezioni – perché è importante?

Se vogliamo visualizzare i dati sulla mappa, dobbiamo lavorare con le proiezioni. Per mostrare qualcosa sul piano 2D, dobbiamo proiettare le coordinate dalle coordinate 3D sul globo. La proiezione più popolare che tutti conosciamo e usiamo è la proiezione di Mercatore (Web Mercator o Mercatore WGS84 per essere precisi, poiché è utilizzata dalla maggior parte dei fornitori di mappe) e il sistema di coordinate più popolare è il sistema geodetico mondiale 1984 – WGS84 (o EPSG:4326). Questo sistema si basa sui gradi e varia in longitudine da -180° a 180° (Ovest a Est) e in latitudine da -90° a 90° (Sud a Nord).

• Segmentare qualsiasi cosa in 3D per le nuvole di punti Guida completa (SAM 3D)
• 6 Lavori remoti di intelligenza artificiale da cercare nel 2024
• Guidare nell’era dell’IA La paura dell’IA cambierà i costi delle vite?

Ogni proiezione sul piano 2D presenta alcune distorsioni. Il Mercatore è una proiezione di mappa conforme, il che significa che gli angoli dovrebbero essere preservati tra gli oggetti sulla Terra. Più alto sopra 0° (o più basso al di sotto di 0°) la latitudine, maggiore è la distorsione nell’area e nella distanza. Poiché il diagramma di Voronoi si basa pesantemente sulla distanza tra i semi, lo stesso errore di distorsione viene trasmesso durante la generazione del diagramma.

La Terra è un ellissoide dalla forma irregolare, ma per i nostri scopi può essere approssimata alla forma di una sfera. Generando il diagramma di Voronoi sulla sfera, possiamo calcolare correttamente la distanza in base agli archi sulla superficie di una sfera. Successivamente, possiamo mappare i poligoni sferici generati sulle coordinate 2D proiettate e possiamo essere sicuri che la linea che separa due celle adiacenti di Voronoi sarà perpendicolare alla linea che connette i due semi che definiscono queste celle.

Di seguito puoi vedere il problema degli angoli e delle distanze che ho descritto in precedenza. Anche se le linee si intersecano nello stesso punto, le forme e gli angoli delle celle di Voronoi differiscono.

Un altro problema è che non è possibile confrontare le regioni in diverse parti del mondo (cioè non posizionate sulla stessa latitudine) se si utilizza un diagramma di Voronoi 2D, poiché le aree saranno fortemente distorte.

L’intero notebook Jupyter con il codice utilizzato negli esempi sottostanti può essere trovato su GitHub. Qui alcune funzioni vengono saltate per brevità.

Prerequisiti

Installa le librerie richieste

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

Importa i moduli e le funzioni richieste

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

Primo esempio

Definiamo sei punti sul globo: i poli nord e sud e quattro punti sull’equatore.

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

Genera il diagramma di Voronoi utilizzando voronoi_diagram dalla libreria Shapely

def generate_flat_voronoi_diagram_regions(    seeds_gdf: gpd.GeoDataFrame,) -> gpd.GeoDataFrame:    points = MultiPoint(seeds_gdf.geometry.values)    # Genera il diagramma 2D    regions = voronoi_diagram(points)    # Mappa le geometrie a GeoDataFrame    flat_voronoi_regions = gpd.GeoDataFrame(        geometry=list(regions.geoms),        crs="EPSG:4326",    )    # Applica gli indici dal dataframe dei seeds    flat_voronoi_regions.index = gpd.pd.Index(        flat_voronoi_regions.sjoin(seeds_gdf)["index_right"],        name="region_id",    )    # Taglia ai confini terrestri    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)

Genera i diagrammi di Voronoi utilizzando VoronoiRegionalizer dalla libreria srai. Sotto il cofano, utilizza l’implementazione SphericalVoronoi dalla libreria scipy e trasforma correttamente le coordinate WGS84 da e verso il sistema di coordinate sferiche.

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

Vediamo la differenza tra i due grafici.

La prima cosa che si può notare è che il diagramma di Voronoi 2D non si ripiega attorno al globo, poiché funziona su un piano cartesiano piatto Cartesiano. Il diagramma di Voronoi sferico copre correttamente la Terra e non si spezza nella linea dell’antimeridiano (dove la longitudine passa da 180 ° a -180 °).

Per quantificare numericamente la differenza possiamo calcolare la metrica IoU (Intersection over Union) (o Jaccard Index) per misurare la differenza tra le forme dei poligoni. Il valore di questa metrica varia tra 0 e 1, dove 0 significa nessuna sovrapposizione e 1 significa sovrapposizione completa.

def calcola_iou(    regioni_piane: gpd.GeoDataFrame, regioni_sferiche: gpd.GeoDataFrame) -> float:    area_intersezioni_totali = 0    area_totali_unite = 0    # Itera tutte le regioni    for indice in regioni_sferiche.indice :        # Trova regione sferica corrispondente e regione Voronoi piatta        geometria_regionale_sferica = regioni_sferiche.loc[indice].geometria        geometria_regionale_piana = regioni_piane.loc[indice].geometria        # Calcola l'area di intersezione        area_intersezioni = geometria_regionale_sferica.intersection(            geometria_regionale_piana        ).area        # Calcola l'area di unione        # Codice alternativo:        # geometria_regionale_sferica.unione(geometria_regionale_piana).area        area_unioni = (           geometria_regionale_sferica.area            + geometria_regionale_piana.area            - area_intersezioni        )        # Aggiungi alle somme totali        area_intersezioni_totali += area_intersezioni        area_totali_unite += area_unioni    # Dividi l'area di intersezione per l'area di unione    return round(area_intersezioni_totali / area_totali_unite, 3)

calcola_iou(regioni_voronoi_piane_terrestri, regioni_voronoi_sferiche_terrestri)

Il valore calcolato è 0,423, che è piuttosto basso e su scala grande, quei poligoni sono diversi l’uno dall’altro, come si può facilmente vedere nei grafici sopra.

Esempio di dati reali: suddivisione del globo usando le posizioni dei defibrillatori esterni automatici (AED)

I dati utilizzati in questo esempio provengono da OpenAEDMap e si basano su dati di OpenStreetMap. Il file preparato ha le posizioni filtrate (esattamente 80694) senza nodi duplicati definiti uno sopra l’altro.

# Carica le posizioni AED in GeoDataFrameaed_world_gdf = gpd.read_file(    "https://raw.githubusercontent.com/RaczeQ/medium-articles/main/articles/spherical-geovoronoi/aed_world.geojson")

Genera i diagrammi di Voronoi per gli AED

aed_flat_voronoi_regions = genera_regioni_diagramma_voronoi_piane(aed_world_gdf)aed_spherical_voronoi_regions = VoronoiRegionalizer(    seed=aed_world_gdf, max_metri_tra_punti=1_000).trasforma()

Confrontiamo questi diagrammi di Voronoi.

La differenza è abbastanza evidente osservando i grafici. Tutti i confini nella versione 2D sono retti mentre quelli sferici sembrano piuttosto curvi nelle coordinate WGS84. Puoi anche vedere chiaramente che nella versione piana, molte regioni convergono ai poli (la proiezione ortogonale si concentra sul polo sud), mentre nella versione sferica non avviene. Un’altra differenza visibile è la continuità intorno all’antimeridiano, che è stato menzionato nel primo esempio. Le regioni che emergono dalla Nuova Zelanda vengono tagliate bruscamente nella versione piana.

Vediamo il valore di IoU:

calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)

Il valore calcolato è 0,511, leggermente migliore rispetto al primo esempio, ma comunque i poligoni corrispondono approssimativamente al 50%.

Avvicinandosi alla scala della città

Vediamo la differenza su una scala più piccola. Possiamo selezionare tutti gli AED che si trovano a Londra e tracciarli.

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

Il valore è 0,675. Sta migliorando, ma è comunque una differenza evidente. Poiché gli AED sono posizionati più densamente, le forme e le distanze si riducono, quindi le differenze tra i diagrammi di Voronoi calcolati nel piano 2D proiettato e su una sfera diminuiscono.

Osserviamo alcuni esempi individuali sovrapposti l’uno sull’altro.

Le aree di quei poligoni corrispondono principalmente, ma puoi vedere le differenze negli angoli e nelle forme. Queste discrepanze potrebbero essere importanti nell’analisi spaziale e potrebbero modificare i risultati. Più grande è l’area di interesse, maggiore diventa la differenza.

Sommario

Spero che ora tu possa capire perché il diagramma di Voronoi sferico è più adatto all’uso nel dominio geospaziale rispetto a quello piatto.

La maggior parte delle analisi nel settore viene attualmente effettuata utilizzando diagrammi di Voronoi in un piano piatto 2D proiettato, il che potrebbe portare a risultati errati.

Per molto tempo non c’è stata una soluzione semplice per i diagrammi di Voronoi sferici disponibili per gli scienziati dei dati geospaziali e gli analisti che lavorano in Python. Ora è facile come installare una libreria. Certo, calcola un po’ più a lungo rispetto alla soluzione piatta poiché deve proiettare punti verso e da coordinate sferiche, mentre ritaglia correttamente i poligoni che si intersecano con l’antimeridiano, ma non dovrebbe importare se desideri preservare la precisione nelle tue analisi. Per gli utenti di JavaScript, è disponibile un’implementazione di diagrammi di Voronoi sferici tramite D3.js.

Disclaimer

Sono uno dei responsabili della libreria srai.