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.

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

)

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.