How to Build Voronoi Service Areas Around Facilities in Python
Problem statement
A Voronoi service area gives each facility the territory that is closer to it, in a straight line, than to any other facility. It is the quickest way to split a region between facilities and count who each one serves. The geometry takes milliseconds; the answer depends on three decisions around it — which facilities go in, what the areas are clipped to, and how people are counted — and on whether straight lines are a fair stand-in for travel.
Measured with the 39 supermarkets in and within 10 km of Chittenden County, Vermont, and the county's 168,323 residents:
- Building and clipping the 38 areas that reach the county's land took 42 milliseconds.
- 10 of those areas belonged to supermarkets outside the county, and 4,149 residents lived in them.
- Counting people by block point rather than by block area moved 4,751 residents (2.8%) between areas.
- For 26.5% of residents another supermarket was nearer by road than the one whose area they lived in.
Quick answer
Build the cells in input order, attach the facilities' attributes, clip to land, and join people by location:
import geopandas as gpd
import shapely
cells = shapely.voronoi_polygons(shapely.MultiPoint(stores.geometry.to_list()),
extend_to=land.buffer(10_000), ordered=True)
areas = gpd.GeoDataFrame(stores.drop(columns="geometry"), geometry=list(cells.geoms), crs=stores.crs).clip(land)
people = gpd.sjoin(homes, areas[["geometry"]], predicate="within").groupby("index_right").POP20.sum()stores holds one projected point per facility (including those just outside the study area), land is the study area without open water, and homes holds one point per populated census block with its population in POP20.
Step-by-step solution
1. Reduce every facility to one point
OpenStreetMap mixes building outlines and points: there were 40 supermarket features, 24 outlines and 16 points. Take a representative point of each outline — always inside the shape, unlike a centroid — and drop points that fall inside an outline, which are the same shop mapped twice. That left 39 supermarkets. Check for facilities very close together too: 2 of them had another supermarket within 50 m, and exact duplicates make ordered=True fail (see fixing Voronoi polygons that extend forever or miss the study area).
2. Include facilities beyond the study area
People near a boundary shop across it. 11 of the 39 supermarkets were outside the county but within 10 km, and 10 of them claimed part of the county's land. Leaving them out would have put 128 blocks, holding 2.5% of residents, in the area of a more distant supermarket inside the county.
3. Build the diagram in a projected CRS, in input order
Use a projected coordinate system so that "closer" means closer on the ground — the county data here are in EPSG:32145, Vermont State Plane in metres. shapely.voronoi_polygons with ordered=True returns the cells in the order of the input points, so cell i belongs to facility i; every cell contained its own supermarket. extend_to stretches the outer cells past the study area instead of stopping at the envelope of the points.
4. Clip to land
The county polygon covers 1,606 km², or 1,534 km² once census blocks that are entirely water are removed. Clip to the land version, or the lakeside areas count open water as territory. After clipping, 38 areas remained, with a median of 18.0 km² and a range from 0.21 km² to 201.2 km².
5. Count the people in each area
A block point inside an area is the fast count: the spatial join took 4 ms. Areal weighting shares each block's population by the fraction of the block's area in each service area, and took 323 ms. The two agree for most areas, but 452 blocks were split between areas and the point count moved 4,751 residents — 2.8% of the county — relative to area weighting. Small areas feel it most: no block point fell in the area of Mac's Market, which held 21 people by area. By point count, the median area held 3,788 residents and the largest, Jake's ONE Market's, 14,012.
6. Report what the areas hide
A service area says nothing about how far its residents are from its facility. The population-weighted median straight-line distance to the area's supermarket was 1,537 m, the 90th percentile 4,857 m and the longest 12,062 m.
7. Check the areas against travel time
Voronoi areas assume straight lines. Assign each block to its nearest supermarket by drive time instead and compare. 566 of the 2,241 blocks, holding 26.5% of residents, were nearer by road to a different supermarket. For most of them the difference was small — a median of 0.53 minutes — but the 90th percentile was 3.60 minutes and the worst 18.6. Store totals moved more than that suggests: 27,646 residents changed store, and 20 of the 38 areas gained or lost more than a quarter of their population.
supermarket area km2 people (points) people (area) people (by road)
Jake's ONE Market 7.5 14,012 13,933 10,110
City Market 3.7 13,658 14,334 15,025
Shaw's 33.7 12,279 11,797 15,654
Hannaford 31.1 10,020 10,027 10,246
Hannaford 126.8 9,884 10,056 9,472
Shelburne Supermarket 167.1 8,310 8,341 9,828
Code examples
Example 1 — one point per facility
import geopandas as gpd
import numpy as np
import pandas as pd
def facility_points(features, tolerance=1.0):
"""Outlines to representative points; drop points inside an outline and near-exact duplicates."""
outlines = features[features.geom_type.isin(["Polygon", "MultiPolygon"])].copy()
points = features[features.geom_type == "Point"]
inside = gpd.sjoin(points, outlines[["geometry"]], predicate="within").index.unique()
outlines["geometry"] = outlines.geometry.representative_point()
out = pd.concat([outlines, points.drop(index=inside)], ignore_index=True)
snapped = pd.DataFrame(np.round(np.column_stack([out.geometry.x, out.geometry.y]) / tolerance))
repeated = snapped.duplicated().to_numpy()
print(f"{len(features)} features -> {len(out) - repeated.sum()} points "
f"(dropped: {len(inside)} inside an outline, {repeated.sum()} duplicates)")
return gpd.GeoDataFrame(out[~repeated], geometry="geometry", crs=features.crs).reset_index(drop=True)On the supermarket layer, in EPSG:32145, it printed 40 features -> 39 points (dropped: 1 inside an outline, 0 duplicates).
Example 2 — areas and their populations
import shapely
def service_areas(facilities, land, reach=10_000):
cells = shapely.voronoi_polygons(shapely.MultiPoint(facilities.geometry.to_list()),
extend_to=land.buffer(reach), ordered=True)
areas = gpd.GeoDataFrame(facilities.drop(columns="geometry"), geometry=list(cells.geoms), crs=facilities.crs)
areas = areas.clip(land)
areas["km2"] = areas.area / 1e6
return areas
def people_by_area(areas, blocks, population="POP20"):
"""Share each block's population by the fraction of its area inside each service area."""
b = blocks[[population, "geometry"]].assign(block_m2=blocks.area)
pieces = gpd.overlay(b, areas[["geometry"]].reset_index(names="area_id"), how="intersection", keep_geom_type=True)
pieces["people"] = pieces[population] * pieces.area / pieces.block_m2
return pieces.groupby("area_id").people.sum().reindex(areas.index, fill_value=0)
land = county.difference(blocks_all[(blocks_all.ALAND20 == 0) & (blocks_all.AWATER20 > 0)].union_all())
areas = service_areas(stores, land)
areas["people"] = people_by_area(areas, populated_blocks)
print(f"{len(areas)} areas, {areas.people.sum():,.0f} people; median area {areas.km2.median():.1f} km2, "
f"median population {areas.people.median():,.0f}")38 areas, 168,323 people; median area 18.0 km2, median population 3,711
By area weighting the median area held 3,711 people, against 3,788 by block point. areas keeps the facilities' index, so its labels match rows of stores, whatever order clip returns them in.
Example 3 — the same homes by road
def weighted_quantile(values, weights, q):
order = np.argsort(values)
cum = np.cumsum(weights[order]) / weights.sum()
return float(values[order][np.searchsorted(cum, q)])
def compare_with_drive_time(homes, areas, minutes, population="POP20"):
"""minutes: travel time from each home (row) to each facility (column, in the facilities' order)."""
joined = gpd.sjoin(homes, areas[["geometry"]], predicate="within", how="left")
by_area = joined[~joined.index.duplicated()].index_right.to_numpy(dtype=int)
by_road = minutes.argmin(axis=1)
w = homes[population].to_numpy(dtype=float)
moved = by_area != by_road
extra = minutes[np.arange(len(minutes)), by_area] - minutes.min(axis=1)
print(f"nearer by road to another facility: {w[moved].sum() / w.sum():.1%} of people; extra minutes for them: "
f"median {weighted_quantile(extra[moved], w[moved], 0.5):.2f}, "
f"90th percentile {weighted_quantile(extra[moved], w[moved], 0.9):.2f}, longest {extra.max():.1f}")
return pd.DataFrame({"voronoi": np.bincount(by_area, weights=w, minlength=minutes.shape[1]),
"by_road": np.bincount(by_road, weights=w, minlength=minutes.shape[1])})
people = compare_with_drive_time(homes, areas, minutes).loc[areas.index]
change = (people.by_road - people.voronoi).abs() / people.voronoi.clip(lower=1)
print(f"areas whose population changes by more than 25%: {(change > 0.25).sum()} of {len(people)}")nearer by road to another facility: 26.5% of people; extra minutes for them: median 0.53, 90th percentile 3.60, longest 18.6
areas whose population changes by more than 25%: 20 of 38
minutes came from one shortest-path search per supermarket on the drive network, as in measuring distance to the nearest facility.
Explanation
What a Voronoi area assumes
That people use the facility nearest in a straight line, that every facility pulls equally, and that there is no capacity. Those assumptions make the areas a clean partition — every point belongs to exactly one — which is why they are useful for maps and territories, and why they are a weak model of behaviour.
Why block points and block areas disagree
A block point stands for everyone in the block. Where blocks are small, the error is small and cancels between neighbours; where blocks are large and rural, a single point can fall in one area while most of the houses are in another. The difference is largest for the smallest areas, where a single block is a large share of the total.
Why road distance disagrees with straight lines
Roads bend around water and hills, bridges are few, and the road network is denser in some directions than others. Road distance to the nearest pharmacy in the same county was a median 1.46 times the straight-line distance. A lake between a home and the nearest store in a straight line can make a store on the same shore the faster choice.
When straight lines are good enough
For a first division of a region, a sales-territory map, or anywhere facilities are dense and the road network is a regular grid, Voronoi areas are quick and defensible. For questions about access or travel burden — who is poorly served, how long trips are — assign by network travel time instead.
Edge cases or notes
- Exact duplicate points with
ordered=TrueraiseGEOSException: Multiple input coordinates in cell; remove them first. GeoSeries.voronoi_polygonsdoes not preserve order in geopandas 1.1; only 2.6% of its cells contained the supermarket in the same row, so attach attributes with a spatial join.- Latitude–longitude coordinates move boundaries. Build in a projected CRS.
- Size and brand are ignored. Use a Huff model when attractiveness matters.
- Empty areas are real. A facility can win no populated block at all.
- A network Voronoi is just nearest by travel time. Assign homes to their fastest facility; there is no polygon step.
- Closed or planned facilities change every neighbour. Rebuild the whole diagram, not one cell.
Internal links
- Fixing Voronoi polygons that extend forever or miss the study area — unbounded, unordered and unclipped cells
- Catchment areas explained: buffers, isochrones and Voronoi compared — when an assignment area is the right catchment
- How to measure distance to the nearest facility for every home — the travel-time assignment
- How to count the population within reach of each facility — overlapping catchments instead of a partition
- How to estimate market share with a Huff model in Python — assignment by probability
- Fixing accessibility scores that are wrong near the study area edge — facilities beyond the boundary
- Location analysis explained: catchments, accessibility and site selection — where service areas fit
- Census geographies explained: blocks, tracts, output areas and why they nest — the population units being counted
FAQ
What is a Voronoi service area?
The part of a region closer, in a straight line, to one facility than to any other. Around 39 supermarkets, the 38 areas reaching Chittenden County's land had a median size of 18.0 km².
How do I keep Voronoi cells matched to the right facility?
Pass ordered=True to shapely's voronoi_polygons, which returns cells in input order, or attach attributes afterwards with a spatial join between cells and facility points.
Should I include facilities outside my study area?
Yes. Ten supermarkets outside the county claimed part of its land, and leaving them out put 2.5% of residents in the area of a more distant store.
How should I count the population of each service area?
Areal weighting of block polygons is the safer count; block points are faster. On the county's areas the point count moved 2.8% of residents between areas.
Are Voronoi areas the same as nearest by road?
No. For 26.5% of residents a different supermarket was nearer by road, usually by under a minute but by up to 18.6 minutes.
When should I use drive time instead of Voronoi areas?
Whenever the question is about travel or access, or the road network is irregular. Twenty of the 38 supermarket areas changed population by more than a quarter when homes were assigned by road.