Earth Isn’t Flat, and Neither Should Your Voronoi Diagrams Be

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

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.

Dimensions & projections — why does it matter?

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.

Prerequisites

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

First example

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.

The first thing that can be seen is that the 2D Voronoi diagram doesn’t loop back around the globe, since it works on a flat Cartesian plane. The Spherical Voronoi diagram properly covers the Earth and doesn’t break at the antimeridian line (where the longitude switches from 180° to -180°).

To numerically quantify the difference we can calculate the IoU (Intersection over Union) metric (or Jaccard Index) to measure the difference between the shapes of the polygons. The value of this metric falls between 0 and 1, where 0 means no overlap and 1 means full overlap.

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
)

The calculated value is 0.423, which is pretty low and on the big scale, those polygons are different from each other, which can be easily seen in the plots above.

Real data example: splitting the globe using AEDs (Automated External Defibrillators) positions

The data used in this example comes from the OpenAEDMap and is based on OpenStreetMap data. The prepared file has filtered positions (80694 to be exact) without duplicated nodes defined on top of each other.

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

Generate Voronoi diagrams for the AEDs

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()

Let’s compare these Voronoi diagrams.

The difference is quite obvious when looking at the plots. All borders in the 2D version are straight while spherical ones look quite bendy in the WGS84 coordinates. You can also clearly see that on the flat version, a lot of regions converge on the poles (orthogonal projection focuses on the south pole), while the spherical one doesn’t. Another visible difference is continuity around antimeridian, which was mentioned in the first example. The regions emerging from New Zealand are abruptly cut on the flat version.

Let’s see the IoU value:

calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)

The calculated value is 0.511, which is slightly better than the first example, but still, the polygons match roughly 50%.

Zooming into the city scale

Let’s see the difference on a smaller scale. We can select all the AEDs that are located in London and plot it.

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

The value is 0.675. It’s getting better, but it still is a noticeable difference. Since the AEDs are placed denser, the shapes and distances are getting smaller, so the differences between Voronoi diagrams calculated in the projected 2D plane and on a sphere diminish.

Let’s look at some individual examples overlayed on top of each other.

The areas of those polygons mostly match, but you can see the differences in angles and shapes. Those discrepancies could be important in the spatial analysis and might change the results of them. The bigger the area of interest, the bigger becomes the difference.

Summary

I hope that now you can see why the spherical Voronoi diagram is better suited for use in the geospatial domain than the flat one.

Most of the analyses in the domain are currently made using Voronoi diagrams in a projected 2D flat plane, which could lead to wrong results.

For a long time, there was no simple solution for spherical Voronoi diagrams available for geospatial data scientists and analysts working in Python. Now it’s as easy as installing one library.
Sure, it calculates a little bit longer than the flat solution since it has to project points to and from spherical coordinates, while properly clipping polygons intersecting antimeridian, but it shouldn’t matter if you want to preserve precision in your analyses.
For JavaScript users, there is an already available spherical Voronoi D3.js implementation.

Disclaimer

I’m one of the maintainers of the srai library.

Earth Isn’t Flat, and Neither Should Your Voronoi Diagrams Be was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.