How to Build Voronoi Service Areas Around Facilities in Python

On this page

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

The diagram itself is one line; the steps around it decide whether the areas mean anything.

What you get — from The data used on this page

41 features: {'Polygon': 25, 'Point': 16}

One point per facility, neighbours included, clipped to land

1. Reduce every facility to one point

OpenStreetMap mixes building outlines and points: there were 41 supermarket features, 25 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, and outlines mapped twice over, which land on the same point. 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 6 ms. Areal weighting shares each block's population by the fraction of the block's area in each service area, and took 352 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. 565 of the 2,239 blocks with a road route, 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.48 minutes — but the 90th percentile was 3.60 minutes and the worst 18.5. 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,504

The same homes, assigned by straight line and by road: Shaw's gains more than a quarter, Jake's ONE Market loses more than a quarter.

Supermarket service areas and their populations

The data used on this page

Chittenden County, Vermont. The county, its water and its people are the TIGER/Line 2020 blocks: the county is the union of its blocks, and each populated block is used both as a polygon and as a representative point with its POP20. Supermarkets are OpenStreetMap shop=supermarket features within 10 km of the county. Drive times come from OpenStreetMap's drive network for the county plus 10.5 km, with OSMnx's imputed free-flow speeds, road-only between nearest nodes; two blocks (553 people) with no road route to any supermarket are left out of the road comparison. 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 pandas as pd
import scipy.sparse as sp
import shapely
from scipy.sparse.csgraph import dijkstra

# 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-service-areas-python-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
county = blocks_all.union_all()
populated_blocks = blocks_all[blocks_all.POP20 > 0].reset_index(drop=True)
homes = populated_blocks.copy()
homes["geometry"] = homes.geometry.representative_point()

reach = gpd.GeoSeries([county.buffer(10_500)], crs=CRS).to_crs(4326).iloc[0]
features = ox.features_from_polygon(reach, tags={"shop": "supermarket"}).to_crs(CRS)
features = features[shapely.distance(features.geometry.representative_point().values, county) <= 10_000]
features = features.sort_index().reset_index()[["element", "id", "name", "geometry"]]
print(f"{len(features)} features: {features.geom_type.value_counts().to_dict()}")

G = ox.add_edge_travel_times(ox.add_edge_speeds(ox.graph_from_polygon(reach, network_type="drive")))
G = ox.project_graph(G, to_crs=CRS)
index = {n: i for i, n in enumerate(G.nodes)}
secs = {}
for u, v, d in G.edges(data=True):                      # shortest of any parallel edges
    secs[(index[u], index[v])] = min(secs.get((index[u], index[v]), np.inf), d["travel_time"])
rows, cols = zip(*secs)
R = sp.csr_array((list(secs.values()), (rows, cols)), shape=(len(index), len(index))).T.tocsr()
h_idx = [index[n] for n in ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y)]

41 features: {'Polygon': 25, 'Point': 16}

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)


stores = facility_points(features)

41 features -> 39 points (dropped: 1 inside an outline, 1 duplicates)

s_idx = [index[n] for n in ox.distance.nearest_nodes(G, stores.geometry.x, stores.geometry.y)]
minutes = dijkstra(R, directed=True, indices=s_idx)[:, h_idx].T / 60      # homes x stores, in stores' order
routable = np.isfinite(minutes).all(axis=1)

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[routable], areas, minutes[routable]).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.48, 90th percentile 3.60, longest 18.5
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.

Example 4 — the other measurements on this page

import time

start = time.perf_counter()
service_areas(stores, land)
print(f"build and clip: {(time.perf_counter() - start) * 1000:.0f} ms")
inside = shapely.distance(stores.geometry.values, county) == 0
outside_areas = areas.index[~inside[areas.index]]
w = homes.POP20.to_numpy(float)
start = time.perf_counter()
joined = gpd.sjoin(homes, areas[["geometry"]], predicate="within")
t_join = time.perf_counter() - start
by_point = joined.groupby("index_right").POP20.sum().reindex(areas.index, fill_value=0)
start = time.perf_counter()
by_area = people_by_area(areas, populated_blocks)
t_area = time.perf_counter() - start
print(f"{len(areas)} areas, {len(outside_areas)} of them for supermarkets outside the county, holding "
      f"{by_point[outside_areas].sum():,} residents; {(~inside).sum()} supermarkets outside in all")
print(f"sjoin {t_join * 1000:.0f} ms, areal weighting {t_area * 1000:.0f} ms; residents moved by point counting "
      f"{(by_point - by_area).abs().sum() / 2:,.0f} ({(by_point - by_area).abs().sum() / 2 / w.sum():.1%})")
split = gpd.overlay(populated_blocks[["geometry"]].reset_index(names="b"), areas[["geometry"]], how="intersection",
                    keep_geom_type=True).groupby("b").size()
print(f"blocks split between areas: {(split > 1).sum()}; areas with no block point: "
      f"{', '.join(stores.name[i] + f' ({by_area[i]:.0f} by area)' for i in areas.index[by_point == 0])}")
print(f"by point: median {by_point.median():,.0f}, largest {stores.name[by_point.idxmax()]} {by_point.max():,}")
print(f"areas: median {areas.km2.median():.1f} km2, from {areas.km2.min():.2f} to {areas.km2.max():.1f}; "
      f"county {county.area / 1e6:,.0f} km2, land {land.area / 1e6:,.0f} km2")

# without the supermarkets outside the county
inner = service_areas(stores[inside], land)
j_in = gpd.sjoin(homes, inner[["geometry"]], predicate="within")
j_all = joined.loc[j_in.index]
changed = j_in.index_right.to_numpy() != j_all.index_right.to_numpy()
print(f"leaving out the outside supermarkets moves {changed.sum()} blocks, "
      f"{homes.POP20[j_in.index[changed]].sum() / w.sum():.1%} of residents")

# straight-line distance to the area's own supermarket
d = homes.loc[joined.index].distance(stores.geometry.reindex(joined.index_right).set_axis(joined.index))
wj = homes.loc[joined.index, "POP20"].to_numpy(float)
print(f"distance to own supermarket: median {weighted_quantile(d.to_numpy(), wj, .5):,.0f} m, "
      f"90th {weighted_quantile(d.to_numpy(), wj, .9):,.0f} m, longest {d.max():,.0f} m")

by_road = minutes[routable].argmin(axis=1)
area_of = joined.index_right.reindex(homes.index).to_numpy()[routable]
print(f"blocks nearer by road to another store: {(by_road != area_of).sum()} of {routable.sum()}; "
      f"residents changing store {(people.by_road - people.voronoi).abs().sum() / 2:,.0f}")
table = pd.DataFrame({"supermarket": stores.name[areas.index].to_numpy(), "km2": areas.km2.round(1).to_numpy(),
                      "points": by_point.to_numpy(), "area": by_area.round().to_numpy(),
                      "by road": people.by_road.round().to_numpy()}).sort_values("points", ascending=False)
print(table.head(6).to_string(index=False))

build and clip: 75 ms
38 areas, 10 of them for supermarkets outside the county, holding 4,149 residents; 11 supermarkets outside in all
sjoin 6 ms, areal weighting 352 ms; residents moved by point counting 4,751 (2.8%)
blocks split between areas: 452; areas with no block point: Mac's Market (21 by area)
by point: median 3,788, largest Jake's ONE Market 14,012
areas: median 18.0 km2, from 0.21 to 201.2; county 1,606 km2, land 1,534 km2
leaving out the outside supermarkets moves 128 blocks, 2.5% of residents
distance to own supermarket: median 1,537 m, 90th 4,857 m, longest 12,062 m
blocks nearer by road to another store: 565 of 2239; residents changing store 27,646
          supermarket   km2  points    area  by road
    Jake's ONE Market   7.5   14012 13933.0  10110.0
          City Market   3.7   13658 14334.0  15025.0
               Shaw's  33.7   12279 11797.0  15654.0
            Hannaford  31.1   10020 10027.0  10246.0
            Hannaford 126.8    9884 10056.0   9472.0
Shelburne Supermarket 167.1    8310  8341.0   9504.0

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.

Boundary and counting choices moved a few per cent of residents; road distance moved a quarter.

Duplicates, ordering and degree coordinates

Exact duplicate points with ordered=True raise GEOSException: Multiple input coordinates in cell; remove them first.
GeoSeries.voronoi_polygons does 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
How to build a morphological tessellation — Voronoi cells around buildings

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