Fixing Voronoi Polygons That Extend Forever or Miss the Study Area #

Problem statement #

You built Voronoi service areas around your facilities, and one of these is true:

>>> vor = scipy.spatial.Voronoi(xy)
>>> sum(-1 in vor.regions[r] for r in vor.point_region)
8

Eight of 29 fire stations have a region that runs off to infinity, and the polygons you can build from the rest cover 51.9% of the county's land. Or you switched to shapely.voronoi_polygons and the opposite happened: the cells cover 6.83 times the county's area. Or the polygons look right and are attached to the wrong facilities.

All of these were measured on real OpenStreetMap facilities in Chittenden County, Vermont: 29 fire stations and 106 schools inside the county, 2,903 census blocks and 168,323 people. The last failure is the dangerous one, because nothing about the map looks wrong — measured, 28 of 29 polygons did not contain the station whose row they were joined to.

Quick answer #

Build the diagram with ordered=True, extend it past the study area, clip to the land you care about, and reattach attributes by location rather than by position:

import geopandas as gpd
import shapely

stations = stations.to_crs(32145)                     # metres, never degrees
cells = shapely.voronoi_polygons(
    shapely.MultiPoint(stations.geometry.values),
    extend_to=county.union_all(),
    ordered=True,                                     # cell i belongs to point i
)
areas = gpd.GeoDataFrame(stations.drop(columns="geometry"),
                         geometry=list(cells.geoms), crs=stations.crs)
areas = areas.clip(land)                              # the county minus the lake
assert areas.contains(stations.geometry).all()

Measured: with ordered=True, 0 of 29 cells failed the containment check, and the clipped cells covered 100% of the county's land.

Step-by-step solution #

1. Project the points before you build the diagram #

Voronoi edges are perpendicular bisectors, and "perpendicular" and "halfway" only mean what you intend in a projected CRS. In degrees at 44.5° north, a degree of longitude is about 0.71 of a degree of latitude on the ground, so the bisectors lean.

Measured, building the fire-station diagram in EPSG:4326 instead of EPSG:32145 assigned 219 of 2,903 blocks — 10,080 people, 6.0% of the county — to a different station. For the 106 schools it was 13.5% of the population.

2. Recognise what scipy's -1 means #

scipy.spatial.Voronoi returns vertices and regions, and a region containing vertex index -1 is unbounded: it has edges that go to infinity. The facilities on the outside of the point set always have one.

from scipy.spatial import Voronoi

vor = Voronoi(xy)
open_regions = [i for i, r in enumerate(vor.point_region)
                if -1 in vor.regions[r] or not vor.regions[r]]

Measured: 8 of 29 fire stations and 9 of 106 schools. Dropping those regions — the usual first attempt — left the station cells covering 51.9% of the land and 78.0% of the population. The missing half is exactly the rural edge of the county, where service areas are largest and matter most.

3. Let shapely close the outer cells for you #

shapely.voronoi_polygons (GEOS) closes every cell against a bounding envelope, so there is nothing infinite to handle:

cells = shapely.voronoi_polygons(shapely.MultiPoint(xy))

The default envelope is the points' extent enlarged generously. Measured, for the fire stations it produced cells totalling 10,960 km² — 6.83× the 1,606 km² county — and for schools 12.55×. That is not a bug; it is an envelope waiting to be clipped.

4. Pass extend_to when the points do not reach the edges #

If the facilities sit well inside the study area, the default envelope might not reach its corners. extend_to guarantees it does:

cells = shapely.voronoi_polygons(mp, extend_to=county_geom, ordered=True)

extend_to only ever enlarges the envelope. Measured, passing a 20 m box left the output unchanged at 10,960 km². It cannot be used to trim the diagram; clipping does that.

5. Use ordered=True, or join by location #

This is the failure that produces a believable wrong map. Without ordered=True, GEOS returns the cells in its own internal order:

cells = shapely.voronoi_polygons(mp)                 # order is not the input order
wrong = sum(not cells.geoms[i].contains(shapely.Point(xy[i])) for i in range(len(xy)))

Measured: 28 of 29 stations and 105 of 106 schools sat outside the cell at their own index. GeoSeries.voronoi_polygons() in GeoPandas 1.1.4 has no ordered argument, and its documentation says the output order does not correspond to the input — measured, the same 28 of 29.

ordered=True needs shapely 2.1 with GEOS 3.12 or newer (this environment: shapely 2.1.2, GEOS 3.13.1). Where it is not available, attach attributes with a spatial join: each cell contains exactly one generating point.

6. Clip to the land, not to the boundary #

A county boundary often includes water. Chittenden County's TIGER polygon is 1,606 km², of which 214 km² — 13.3% — is water, nearly all of it Lake Champlain. Clipping to the boundary leaves lakeside cells with an area, and a population density, they do not have.

land = county.union_all().difference(water.union_all())   # OSM natural=water
areas = areas.clip(land)

Measured, the land left was 1,392 km², and the clipped station cells covered all of it: 100.0%. Beware of using census blocks with land area as the land mask — measured, their union still contained 142 km² of lake, because lakeside blocks carry water as well as land.

7. Remove duplicate points first #

Two facilities at the same coordinate — a station mapped as both a node and a building, say — share one cell:

>>> len(shapely.voronoi_polygons(shapely.MultiPoint(dup)).geoms)   # 32 input points
29
>>> shapely.voronoi_polygons(shapely.MultiPoint(dup), ordered=True)
shapely.errors.GEOSException: Multiple input coordinates in cell at 442497 221006

Without ordered=True the count silently drops, and any positional join is now off by three. With it, GEOS refuses and names the coordinate. Scipy kept 32 entries in point_region but only 29 distinct regions.

Code examples #

Example 1 — service areas that are bounded, ordered and clipped #

import geopandas as gpd
import numpy as np
import shapely


def voronoi_service_areas(facilities, clip_to, id_col="name", crs=None):
    """One polygon per facility, covering clip_to exactly, attributes kept."""
    if facilities.crs is None or facilities.crs.is_geographic:
        if crs is None:
            raise ValueError("project the facilities first, or pass a projected crs")
        facilities = facilities.to_crs(crs)

    pts = facilities.copy()
    pts["geometry"] = pts.geometry.representative_point()
    before = len(pts)
    pts = pts[~pts.geometry.duplicated()].reset_index(drop=True)
    if len(pts) < before:
        print(f"dropped {before - len(pts)} duplicate locations")

    boundary = clip_to.to_crs(pts.crs).union_all()
    cells = shapely.voronoi_polygons(shapely.MultiPoint(pts.geometry.values),
                                     extend_to=boundary, ordered=True)
    areas = gpd.GeoDataFrame(pts.drop(columns="geometry"),
                             geometry=list(cells.geoms), crs=pts.crs)

    misplaced = ~areas.contains(pts.geometry)
    if misplaced.any():
        raise RuntimeError(f"{misplaced.sum()} cells do not contain their facility")

    areas = areas.clip(boundary)
    coverage = areas.area.sum() / boundary.area
    print(f"{len(areas)} areas covering {coverage:.1%} of the study area")
    return areas

On the county's fire stations with the land polygon as clip_to: 29 areas covering 100.0%.

Example 2 — the fallback when ordered=True is unavailable #

import geopandas as gpd
import shapely


def voronoi_by_join(points):
    """Build unordered cells, then give each cell the attributes of the point inside it."""
    cells = gpd.GeoDataFrame(
        geometry=list(shapely.voronoi_polygons(shapely.MultiPoint(points.geometry.values)).geoms),
        crs=points.crs,
    )
    joined = gpd.sjoin(cells, points, predicate="contains", how="left")
    if joined.index.duplicated().any():
        raise ValueError("a cell contains more than one point — remove duplicates first")
    if joined["index_right"].isna().any():
        raise ValueError("a cell contains no point — check the CRS of both layers")
    return joined.drop(columns="index_right")

The two checks are what make it safe: duplicates give a cell two points, and a CRS mismatch gives cells none.

Example 3 — auditing a diagram you were given #

import numpy as np
import geopandas as gpd


def audit_service_areas(areas, facilities, study_area, population=None, pop_col="POP20"):
    """Four numbers that reveal every failure in this article.
    areas and facilities must share index labels (one row per facility)."""
    boundary = study_area.to_crs(areas.crs).union_all()
    fac = facilities.to_crs(areas.crs)

    total = areas.area.sum() / boundary.area
    covered = areas.union_all().intersection(boundary).area / boundary.area
    shared = areas.index.intersection(fac.index)
    own = areas.loc[shared].geometry.contains(fac.loc[shared].geometry, align=True).mean()
    print(f"cells: {len(areas)}  facilities: {len(fac)}  matched labels: {len(shared)}")
    print(f"total area / study area: {total:.2f}   (> 1.0 means unclipped)")
    print(f"study area covered:      {covered:.1%} (< 100% means dropped cells)")
    print(f"cells containing their own facility: {own:.1%} (< 100% means mislabelled)")

    if population is not None:
        pts = population.to_crs(areas.crs).copy()
        pts["geometry"] = pts.geometry.representative_point()
        inside = gpd.sjoin(pts, areas[["geometry"]], predicate="within")
        print(f"population assigned: {inside[pop_col].sum() / pts[pop_col].sum():.1%}")

Run on the unordered school diagram with attributes attached by position, it reported a total-area ratio of 12.55 and 0.9% of cells containing their own school — two failures from one call. On the output of Example 1 it reported 1.00 and 100.0%. The comparison is by index label, not position, because clip reorders rows.

Explanation #

Why Voronoi regions are infinite in the first place #

A Voronoi cell is the set of locations closer to one point than to any other. For a facility on the outside of the set — the convex hull — there is always a direction in which you can walk forever and stay closest to it. The region is genuinely unbounded; scipy reports that honestly with vertex -1, and it is the caller's job to close it.

That is why dropping the open regions removes precisely the peripheral facilities, and with them the largest, most rural service areas: measured, half the county's land.

Why GEOS returns the cells out of order #

GEOS builds the diagram from a Delaunay triangulation, using an incremental algorithm that inserts points in whatever order suits the triangulation, and it emits the cells in the order its internal structure yields them. Mapping back to input order is extra work, which is why it became an option (ordered=True, GEOS 3.12) rather than the default. The consequence is silent: the output is a valid, complete diagram with the wrong labels.

Why the CRS changes which facility wins #

The bisector between two points is where distances to both are equal. If distance is measured in degrees, "equal" is equal in degrees, and at this latitude an east–west degree is 29% shorter on the ground than a north–south one. Blocks near an east–west boundary between two stations switch sides. Measured, that was 6.0% of the population for 29 fire stations and 13.5% for 106 schools — more facilities means more boundaries to get wrong.

Why a Voronoi cell is still only an approximation of a service area #

Every fix here makes the diagram correct as a diagram. It remains a straight-line model: a lake, a river with one bridge, or a ridge will put people in the "wrong" cell of the right diagram. Clipping to land removes the water from the areas but does not route around it. Where travel matters, allocate blocks to their nearest facility along the road network instead.

Edge cases or notes #

• extend_to cannot shrink the diagram. Measured, a 20 m box changed nothing; clip afterwards.
• GeoDataFrame.clip reorders rows. Measured, the 29 clipped cells came back in scrambled order with their attributes intact; a positional comparison afterwards matched 1 of 29. Compare by index label.
• Polygons and lines as input contribute every vertex as a site. Convert facilities to representative_point() first, or a building footprint becomes a dozen sites.
• Collinear points produce parallel strips rather than cells; with only two facilities every cell is a half-plane.
• tolerance snaps nearby sites together, which reduces the cell count just as duplicates do.
• Clipping to a coastline or lake shore can split one cell into several parts. Keep the MultiPolygon, or you will lose islands.
• Scipy is still the right tool for Voronoi in more than two dimensions, or when you need the ridge structure itself.
• Buffering the study area before building the diagram, and including facilities just outside it, stops cells along the boundary being too large; see the edge-effects fix.

Internal links #

• How to build Voronoi service areas around facilities in Python — the how-to this fixes
• Catchment areas explained: buffers, isochrones and Voronoi compared — when a Voronoi cell is the right catchment
• Fixing accessibility scores that are wrong near the study area edge — facilities outside the boundary
• How to measure distance to the nearest facility for every home — the network alternative
• How to clip spatial data in Python with GeoPandas — the clipping step
• How to choose the right projected CRS for your study area — why the diagram must be projected
• How to perform a spatial join in GeoPandas — attaching attributes by location
• Network distance vs straight-line distance explained — the limits of any Voronoi cell

FAQ #

Why do some of my scipy Voronoi regions contain -1? #

Index -1 marks a vertex at infinity. Facilities on the outside of the point set have unbounded regions; measured, 8 of 29 fire stations did. Use shapely.voronoi_polygons, which closes them.

Why are my Voronoi polygons attached to the wrong facilities? #

GEOS returns cells in its own order. Measured, 28 of 29 cells did not contain the point at the same index. Pass ordered=True, or join attributes by location.

How do I stop Voronoi polygons extending far beyond my study area? #

Clip them. The default envelope produced cells 6.83 times the county area; extend_to only enlarges the envelope, so it cannot trim it.

Does the coordinate system matter for Voronoi polygons? #

Yes. Built in degrees, the diagram assigned 6.0% of the county's population to a different fire station than the same diagram in metres. Always project first.

Why do I get fewer polygons than points? #

Duplicate coordinates share a cell. Three duplicates turned 32 points into 29 polygons silently; with ordered=True GEOS raised an error naming the coordinate instead.

Should I clip to the county boundary or to the land? #

To the land when the boundary contains water. Chittenden County's boundary is 13.3% water, and lakeside cells clipped to it overstate their area.