You is perhaps acquainted with Voronoi diagrams and their uses within the geospatial analyses. If not, here is the fast TL;DR: it divides the plane into regions consisting of all points of the plane closer to a given seed than to every other. It is called after mathematician Georgy Voronoy. You may read more about it on the Wikipedia.
How does it apply to the geospatial domain? Using Voronoi diagrams you may quickly find the closest public transit stop for inhabitants of a given city at a much bigger scale, faster than calculating it individually for every constructing individually. Or you may as well use it for instance available in the market share evaluation between different brands.
On this post, I would like to indicate the differences between the everyday Voronoi diagram calculated with projected coordinates on a flat plane and the spherical one, and hopefully, show the latter’s superiority.
If we would like to see data on the map, now we have to work with projections. To point out something on the 2D plane, now we have to project the coordinates from the 3D coordinates on the globe.
The preferred projection that everyone knows and use is the Mercator projection (Web Mercator or WGS84 Mercator to be precise, because it’s utilized by many of the map providers) and the preferred coordinate system is World Geodetic System 1984 — WGS84 (or EPSG:4326). This method relies 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 suggests that angles needs to be preserved between objects on the Earth. The upper above 0° (or lower below 0°) the latitude, the larger the distortion in the realm and the distance. Since the Voronoi diagram heavily relies on the gap between the seeds, the identical distortion error is forwarded when generating the diagram.
The Earth is an irregularly shaped ellipsoid, but for our purposes, it could actually be approximated by the sphere shape. By generating the Voronoi diagram on the sphere, we are able to properly calculate the gap based on the arcs on the surface of a sphere. Later, we are able to map the generated spherical polygons to the projected 2D coordinates and we are able to make certain that the road separating two adjoining Voronoi cells will likely be perpendicular to the road connecting the 2 seeds defining these cells.
Below you may see the angles and distances problem I’ve described above. Though the lines cross at the identical point, Voronoi cells’ shapes and angles differ.
One other problem is you can’t compare the regions in several parts of the world (i.e. not laying on the identical latitude) for those who use a 2D Voronoi diagram for the reason that areas will likely be heavily distorted.
Full Jupyter notebook with code utilized in the examples below could 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 pxfrom 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 4 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",
)
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 2 on the plots.