Why Spherical GeoVoronoi diagrams are better?


A story about precision, unveiling the power of spherical geospatial Voronoi diagrams with Python

An animation showing Earth changing state between two projections.
Earth with Spherical Voronoi diagram moving between 2 projections: Orthogonal and Equirectangular. Generated by the author using the D3.js library.

You might be familiar with Voronoi diagrams and their uses in the geospatial analyses. If not, here is the quick TL;DR: it divides the plane into regions consisting of all points of the plane closer to a given seed than to any other. It is named after mathematician Georgy Voronoy. You can read more about it on the Wikipedia.

How does it apply to the geospatial domain? Using Voronoi diagrams you can quickly find the closest public transit stop for inhabitants of a given city at a bigger scale, faster than calculating it individually for each building separately. Or you can also use it for example in the market share analysis between different brands.

In this post, I want to show the differences between the typical Voronoi diagram calculated with projected coordinates on a flat plane and the spherical one, and hopefully, show the latter’s superiority.

If we want to see data on the map, we have to work with projections. To show something on the 2D plane, we have to project the coordinates from the 3D coordinates on the globe.
The most popular projection that we all know and use is the Mercator projection (Web Mercator or WGS84 Mercator to be precise, since it’s used by most of the map providers) and the most popular coordinate system is World Geodetic System 1984 — WGS84 (or EPSG:4326). This system is based on degrees and it ranges in longitude from -180° to 180° (West to East) and in latitude from -90° to 90° (South to North).

Each projection to the 2D plane has some distortions. The Mercator is a Conformal map projection which means that angles should be preserved between objects on the Earth. The higher above 0° (or lower below 0°) the latitude, the bigger the distortion in the area and the distance. Because the Voronoi diagram heavily relies on the distance between the seeds, the same distortion error is forwarded when generating the diagram.

The Earth is an irregularly shaped ellipsoid, but for our purposes, it can be approximated by the sphere shape. By generating the Voronoi diagram on the sphere, we can properly calculate the distance based on the arcs on the surface of a sphere. Later, we can map the generated spherical polygons to the projected 2D coordinates and we can be sure that the line separating two adjacent Voronoi cells will be perpendicular to the line connecting the two seeds defining these cells.

Below you can see the angles and distances problem I’ve described above. Even though the lines cross at the same point, Voronoi cells’ shapes and angles differ.

Chart showing the difference in angles and distances between 2D and Spherical Voronoi diagrams.
Example of angles and distances difference in both versions of Voronoi diagram. Image by the author.

Another problem is that you can’t compare the regions in different parts of the world (i.e. not laying on the same latitude) if you use a 2D Voronoi diagram since the areas will be heavily distorted.

Full Jupyter notebook with code used in the examples below can be found on GitHub. Here some functions are skipped for brevity.

Install required libraries

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

Import required modules and functions

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

Let’s define six points on the globe: the North and South Poles, and four points on the equator.

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",
)
Image by the author.

Generate Voronoi diagram using voronoi_diagram from the Shapely library

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
)

Generate Voronoi diagrams using VoronoiRegionalizer from the srai library.
Under the hood, it uses the SphericalVoronoi implementation from the scipy library and properly transforms the WGS84 coordinates to and from the spherical coordinate system.

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

Let’s see the difference between the two on the plots.



Source link

Leave a Comment