🔧 Fix Location & accessibility Geometry & topology Location & Accessibility · step 16 Notebook By · Updated

Fixing Voronoi Polygons That Extend Forever or Miss the Study Area

On this page

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 of the 29 fire stations

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.

Three treatments of the same Voronoi diagram: open regions dropped, unclipped envelope, clipped to land. same 29 fire stations, three outcomesscipy, open regions droppedgaps at the edge51.9% of land coveredshapely default envelope6.83× the countyordered, clipped to land100% of land, 0 mislabelledThe outer stations are the ones with open regions — and the largest rural service areas spatialworkflow.io
The first two are the same diagram; only the treatment of the outer cells differs.

Project, close the outer cells, extend and clip

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 77.9% 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. A 20 m box among the stations left it at 10,960 km²; one at the CRS origin stretched it to 137,276 km². Clipping trims; extend_to cannot.

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 145 km² of OpenStreetMap water, 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:

xy = np.column_stack([stations.geometry.x, stations.geometry.y])
dup = np.vstack([xy, xy[:3]])                        # three stations mapped twice
print(len(shapely.voronoi_polygons(shapely.MultiPoint(dup)).geoms), "cells from", len(dup), "points")
try:
    shapely.voronoi_polygons(shapely.MultiPoint(dup), ordered=True)
except shapely.errors.GEOSException as error:
    print(f"GEOSException: {error}")
vor_dup = Voronoi(dup)
print(f"scipy: {len(vor_dup.point_region)} entries in point_region, {len(set(vor_dup.point_region))} distinct regions")
29 cells from 32 points
GEOSException: Multiple input coordinates in cell at 442497 221006
scipy: 32 entries in point_region, 29 distinct regions

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.

How many of 29 Voronoi cells fail to contain the fire station whose row they were joined to. cells not containing their own station, of 29shapely, default28 of 29GeoSeries method28 of 29clip, then by position28 of 29ordered=True0 of 29compared by label0 of 29029Nothing on the map looks wrong — only the attributes are spatialworkflow.io
A positional join after an unordered diagram gives every station a neighbour’s area and population.

Bounded, ordered and clipped service areas

The data used on this page

Chittenden County, Vermont. The county is the union of its TIGER/Line 2020 blocks, all of which are used for population (POP20). Facilities are OpenStreetMap amenity=fire_station and amenity=school features inside the county, by representative point, with stations mapped twice at the same point merged; water is OpenStreetMap natural=water. Everything is in NAD83 / Vermont (EPSG:32145). Overpass via the maps.mail.ru mirror while overpass-api.de was refusing connections. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import requests
import urllib.request
import zipfile

import geopandas as gpd
import numpy as np
import osmnx as ox
import scipy.spatial
import shapely
from scipy.spatial import Voronoi

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-voronoi-polygons-unbounded-or-clipped-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.overpass_rate_limit = False
ox.settings.http_user_agent = "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"
CRS = 32145                                               # NAD83 / Vermont, metres

if not os.path.exists("tl_2020_50_tabblock20.zip"):
    urllib.request.urlretrieve("https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/"
                               "tl_2020_50_tabblock20.zip", "tl_2020_50_tabblock20.zip")
vt_blocks = gpd.read_file("tl_2020_50_tabblock20.zip").to_crs(CRS)
blocks_all = vt_blocks[vt_blocks.COUNTYFP20 == "007"]                   # Chittenden
blocks = blocks_all.reset_index(drop=True)
county = blocks.union_all()
county_gdf = gpd.GeoDataFrame(geometry=[county], crs=CRS)
county_ll = county_gdf.to_crs(4326).iloc[0].geometry


def osm_points(tags):
    f = ox.features_from_polygon(county_ll, tags=tags).to_crs(CRS)
    f["geometry"] = f.geometry.representative_point()
    f = f[shapely.contains(county, f.geometry.values)].sort_index().reset_index()
    return f[~f.geometry.duplicated()].reset_index(drop=True)[["element", "id", "name", "geometry"]]


stations = osm_points({"amenity": "fire_station"})
schools = osm_points({"amenity": "school"})
water = ox.features_from_polygon(county_ll, tags={"natural": "water"}).to_crs(CRS)
water = water[water.geom_type.isin(["Polygon", "MultiPolygon"])]
land = county.difference(water.union_all())
land_gdf = gpd.GeoDataFrame(geometry=[land], crs=CRS)
print(f"{len(stations)} fire stations, {len(schools)} schools, {len(blocks):,} blocks, {blocks.POP20.sum():,} people")
29 fire stations, 106 schools, 2,903 blocks, 168,323 people

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


station_areas = voronoi_service_areas(stations, land_gdf)
29 areas covering 100.0% of the study area

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%}")


unordered = gpd.GeoDataFrame(schools.drop(columns="geometry"), crs=schools.crs,
                             geometry=list(shapely.voronoi_polygons(shapely.MultiPoint(schools.geometry.values)).geoms))
audit_service_areas(unordered, schools, county_gdf)
audit_service_areas(station_areas, stations, land_gdf, blocks)
cells: 106  facilities: 106  matched labels: 106
total area / study area: 12.55   (> 1.0 means unclipped)
study area covered:      100.0% (< 100% means dropped cells)
cells containing their own facility: 0.9% (< 100% means mislabelled)
cells: 29  facilities: 29  matched labels: 29
total area / study area: 1.00   (> 1.0 means unclipped)
study area covered:      100.0% (< 100% means dropped cells)
cells containing their own facility: 100.0% (< 100% means mislabelled)
population assigned: 99.9%

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.

Example 4 — the measurements quoted above

xy = np.column_stack([stations.geometry.x, stations.geometry.y])
sxy = np.column_stack([schools.geometry.x, schools.geometry.y])
land_area, county_area = land.area, county.area
pts = blocks.representative_point()
w = blocks.POP20.to_numpy(float)


def owner(cells_geoms, points):
    """Index of the cell containing each point, by location."""
    out = np.full(len(points), -1)
    for i, c in enumerate(cells_geoms):
        out[shapely.contains(c, points.values)] = i
    return out


# degrees against metres
for label, layer in (("stations", stations), ("schools", schools)):
    ll = layer.to_crs(4326)
    metres = shapely.voronoi_polygons(shapely.MultiPoint(layer.geometry.values), ordered=True, extend_to=county.buffer(20_000))
    degrees = shapely.voronoi_polygons(shapely.MultiPoint(ll.geometry.values), ordered=True,
                                       extend_to=county_gdf.to_crs(4326).union_all().buffer(0.3))
    a = owner(metres.geoms, pts)
    b = owner(degrees.geoms, pts.to_crs(4326))
    moved = a != b
    print(f"{label}: built in degrees, {moved.sum()} of {len(blocks):,} blocks move "
          f"({w[moved].sum():,.0f} people, {w[moved].sum() / w.sum():.1%})")

# scipy's open regions
for label, arr in (("stations", xy), ("schools", sxy)):
    v = Voronoi(arr)
    opened = [i for i, r in enumerate(v.point_region) if -1 in v.regions[r] or not v.regions[r]]
    print(f"{label}: {len(opened)} of {len(arr)} regions unbounded")
v = Voronoi(xy)
closed = [shapely.Polygon(v.vertices[v.regions[r]]) for r in v.point_region if -1 not in v.regions[r] and v.regions[r]]
kept = shapely.union_all(closed).intersection(land)
inside = shapely.contains(kept, pts.values)
print(f"closed scipy cells cover {kept.area / land_area:.1%} of the land and {w[inside].sum() / w.sum():.1%} of the people")

# shapely's default envelope, and extend_to
for label, arr in (("stations", xy), ("schools", sxy)):
    cells = shapely.voronoi_polygons(shapely.MultiPoint(arr))
    print(f"{label}: unclipped cells total {cells.area / 1e6:,.0f} km², {cells.area / county_area:.2f}x the county")
cx, cy = xy.mean(axis=0)
tiny = shapely.voronoi_polygons(shapely.MultiPoint(xy), extend_to=shapely.box(cx - 10, cy - 10, cx + 10, cy + 10))
far = shapely.voronoi_polygons(shapely.MultiPoint(xy), extend_to=shapely.box(0, 0, 20, 20))
print(f"extend_to a 20 m box among the stations: {tiny.area / 1e6:,.0f} km²; "
      f"a 20 m box at the CRS origin (0, 0): {far.area / 1e6:,.0f} km²")

# order
for label, arr in (("stations", xy), ("schools", sxy)):
    cells = shapely.voronoi_polygons(shapely.MultiPoint(arr))
    wrong = sum(not cells.geoms[i].contains(shapely.Point(arr[i])) for i in range(len(arr)))
    ordered = shapely.voronoi_polygons(shapely.MultiPoint(arr), ordered=True)
    right = sum(not ordered.geoms[i].contains(shapely.Point(arr[i])) for i in range(len(arr)))
    print(f"{label}: {wrong} of {len(arr)} outside their own cell unordered, {right} with ordered=True")
gs = stations.geometry.voronoi_polygons()
print(f"GeoSeries.voronoi_polygons: {sum(not gs.iloc[i].contains(stations.geometry.iloc[i]) for i in range(len(stations)))} "
      f"of {len(stations)} (geopandas {gpd.__version__}, shapely {shapely.__version__}, GEOS {shapely.geos_version_string})")

# water
water_in = county.intersection(water.union_all())
print(f"county {county_area / 1e6:,.0f} km², water {water_in.area / 1e6:,.0f} km² ({water_in.area / county_area:.1%}), "
      f"land {land_area / 1e6:,.0f} km²")
land_blocks = blocks_all[blocks_all.ALAND20 > 0].union_all()
print(f"union of blocks with land area still holds {land_blocks.intersection(water.union_all()).area / 1e6:,.0f} km² of OSM water")
stations: built in degrees, 219 of 2,903 blocks move (10,080 people, 6.0%)
schools: built in degrees, 324 of 2,903 blocks move (22,804 people, 13.5%)
stations: 8 of 29 regions unbounded
schools: 9 of 106 regions unbounded
closed scipy cells cover 51.9% of the land and 77.9% of the people
stations: unclipped cells total 10,960 km², 6.83x the county
schools: unclipped cells total 20,157 km², 12.55x the county
extend_to a 20 m box among the stations: 10,960 km²; a 20 m box at the CRS origin (0, 0): 137,276 km²
stations: 28 of 29 outside their own cell unordered, 0 with ordered=True
schools: 105 of 106 outside their own cell unordered, 0 with ordered=True
GeoSeries.voronoi_polygons: 28 of 29 (geopandas 1.1.4, shapely 2.1.2, GEOS 3.13.1)
county 1,606 km², water 214 km² (13.3%), land 1,392 km²
union of blocks with land area still holds 145 km² of OSM water

Why Voronoi cells are infinite

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.

Four checks for a Voronoi service-area layer, and the positional join that breaks it. four checks before the areas are usedbuilt in a projected CRSdegrees moved 6.0% of peopleordered=True, or joined by location28 of 29 wrong otherwiseduplicate locations removed32 points gave 29 cellsclipped to land, not the boundarythe county is 13.3% waterrows compared by position after clipclip reorders rowsEach check is one line; together they take about a second spatialworkflow.io
The four checks together take a second on thousands of cells.

Extents, reordering clips and collinear points

  • 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.

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.