--- title: "Voronoi (Thiessen) polygons in Python: catchment rainfall from EA rain gauges and fire station areas in Cumbria" description: "Build Voronoi (Thiessen) polygons in Python with shapely and GeoPandas: Thiessen catchment rainfall from EA gauges and Voronoi vs drive-time service areas." url: https://www.technicalgis.com/tutorials/voronoi-thiessen-polygons-python/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-11-03 updated: 2026-11-03 language: en-GB --- # Voronoi (Thiessen) polygons in Python: catchment rainfall from EA rain gauges and fire station areas in Cumbria > Build Voronoi (Thiessen) polygons in Python with shapely and GeoPandas: Thiessen catchment rainfall from EA gauges and Voronoi vs drive-time service areas. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/voronoi-thiessen-polygons-python/) ## Key takeaways 1. shapely.voronoi\_polygons(MultiPoint, extend\_to=boundary, ordered=True) returns one cell per point in input order (shapely 2.1+), so the cells can go straight back onto your GeoDataFrame. Intersect them with the boundary to get Thiessen polygons. GeoSeries.voronoi\_polygons does not keep the order. 2. Use a projected CRS. Voronoi cells are built from straight-line distances, so latitude and longitude give the wrong shapes; British National Grid (EPSG:27700) is right for Great Britain. 3. Include points outside the boundary. For the River Eden to Sheepmount (2,276 km²), 31 of 47 gauges have a cell inside the catchment and 13 of those 31 stand outside it. Leaving them out would give the edge gauges far too much weight. 4. Check the gauges first. Two EA gauges recorded 726 and 524 mm in 28 days against about 115 mm at their neighbours. Left in, they raise the Thiessen catchment rainfall from 141 to 191 mm. 5. The method matters: 141 mm by Thiessen, 160 mm by IDW and 167 mm for the plain average of the gauges inside, because several gauges cluster in the wet south-west fells. Dropping any one gauge moves the Thiessen figure by at most 3.6 mm. 6. A Voronoi cell is the nearest point in a straight line, not by road. For Cumbria's 38 fire stations, 4.8% of residents (24,090 people) live in a Voronoi cell whose station is not their quickest by road; for 4,680 of them it is more than 5 minutes slower. ## Before you start - Comfortable with GeoPandas basics: GeoDataFrames, CRSs and overlays (see the spatial join tutorial) - Python 3.10+ with geopandas, shapely 2.1+, scipy, osmnx, osmium, networkx, matplotlib and requests - An internet connection: the first run downloads about 60 MB from the Environment Agency, the NRFA, the ONS, Nomis and Geofabrik, and takes a few minutes **Short answer:** to make Voronoi (Thiessen) polygons in Python, put your points in a projected CRS, pass them to **`shapely.voronoi_polygons`** as one `MultiPoint` with `ordered=True`, and intersect the cells with your boundary: ```python import geopandas as gpd, shapely pts = gpd.read_file("gauges.gpkg").to_crs(27700) # metres, not degrees boundary = gpd.read_file("catchment.gpkg").to_crs(27700).union_all() cells = shapely.voronoi_polygons(shapely.MultiPoint(list(pts.geometry)), extend_to=boundary, ordered=True) # one cell per point, same order pts["cell"] = gpd.GeoSeries(list(cells.geoms), index=pts.index, crs=pts.crs).intersection(boundary) pts["weight"] = pts["cell"].area / boundary.area # Thiessen weights, summing to 1 ``` `ordered=True` (shapely 2.1 or later) is the important part: it returns the cells in the order of the input points, so each one goes back on the right row. GeoPandas’ own `GeoSeries.voronoi_polygons()` builds the same diagram but, as its documentation says, does not keep the input order, so you would have to match the cells back with a spatial join. The rest of this tutorial uses that function on two real UK problems: 1. **Catchment rainfall.** The classic Thiessen use: weight 47 Environment Agency rain gauges by area to get the rainfall over the **River Eden** catchment above Carlisle for four weeks of autumn 2026, and see what the gauge average and IDW give instead. 2. **Nearest fire station.** Voronoi cells for Cumbria’s **38 fire stations**, compared with the station that is actually quickest to reach by road, to show where straight-line “nearest” goes wrong. ![Four panels. a) 47 blue gauge points scattered inside and around the green outline of the River Eden catchment. b) The same points joined by grey lines into 83 Delaunay triangles. c) Red points at the triangle circumcentres, joined into red Voronoi edges that run halfway between neighbouring gauges and out past the edges of the frame. d) The Voronoi cells clipped to the catchment outline and shaded blue by each gauge's rainfall total, darkest in the south-west; the unclipped cells outside are drawn in thin grey.](https://www.technicalgis.com/_astro/02-voronoi-construction.DKeEAszZ_Z1EOMO1.webp) _Building Thiessen polygons step by step for the 47 EA rain gauges in and within 15 km of the Eden catchment: a) the points and the boundary, b) the Delaunay triangulation, c) the triangle circumcentres, which join up into the Voronoi edges, d) the cells clipped to the catchment, shaded by each gauge's 28-day total._ You will learn how to: 1. Get a catchment boundary from the NRFA and rain gauges and readings from the EA flood-monitoring API. 2. Check gauge completeness and flag faulty gauges by comparing each one with its Delaunay neighbours. 3. Build Voronoi cells from the Delaunay triangulation and clip them to a boundary. 4. Compute Thiessen weights and the catchment rainfall, and compare them with the gauge average and IDW. 5. Test how sensitive the result is to each gauge. 6. Build Voronoi “nearest station” areas and compare them with drive-time service areas from OSMnx and SciPy’s Dijkstra. 7. Count the people a straight-line Voronoi assigns to the wrong station, using Census 2021 output areas. ## Set-up ```bash py -3.12 -m venv .venv-voronoi .venv-voronoi\Scripts\activate pip install -r requirements.txt python voronoi_thiessen.py ``` On macOS or Linux use `python3.12 -m venv .venv-voronoi` and `source .venv-voronoi/bin/activate`. Everything is in [`voronoi_thiessen.py`](https://www.technicalgis.com/code/voronoi-thiessen-polygons-python/), which caches each download in `data/`, so a second run takes under a minute. Everything is done in **British National Grid (EPSG:27700)**. A Voronoi cell is defined by straight-line distance, and the cell areas become weights, so both need a CRS in metres. In longitude and latitude a degree east is only about 0.57 of a degree north at Carlisle, and the cells come out the wrong shape. See [coordinate reference systems in Python](https://www.technicalgis.com/tutorials/coordinate-reference-systems-python/) if you need a refresher. ## Part 1: Thiessen-weighted catchment rainfall ### Step 1: the catchment and the gauges The catchment is the **River Eden to the Sheepmount gauging station** in Carlisle (NRFA station 76007), which drains the Eden valley, the northern Pennines and the eastern Lake District fells. The National River Flow Archive publishes its boundary through a WFS. The EA’s flood-monitoring API lists every rain gauge within a radius of a point: `voronoi_thiessen.py` ```python cat = get_json(NRFA_WFS, DATA / f"nrfa_{NRFA_ID}.geojson") catchment = gpd.GeoDataFrame.from_features(cat["features"], crs=CRS) basin = shapely.make_valid(catchment.geometry.iloc[0]) lon, lat = Transformer.from_crs(CRS, 4326, always_xy=True).transform(*basin.centroid.coords[0]) stations = get_json(f"{EA}/id/stations?parameter=rainfall&lat={lat:.4f}&long={lon:.4f}&dist=80&_limit=1000", DATA / "ea_rain_stations.json")["items"] # ... convert each station's lat/long to British National Grid ... gauges = gauges[gauges.distance(basin) <= GAUGE_BUFFER_M] # 15 km ``` `Output — Step 1` ```text Eden at Sheepmount (NRFA 76007): polygon area 2,276.0 km2 (NRFA catchment area 2,286.5 km2) EA rainfall stations within 80 km of the catchment centroid: 125; inside or within 15 km of the catchment: 61 ``` Keep the gauges **outside** the catchment. A Voronoi cell depends on all of a point’s neighbours, so a gauge 5 km beyond the watershed can own a strip of the catchment and stops the nearest inside gauge from claiming land it does not represent. A 15 km buffer was enough here; check that cells from the outermost ring do not reach the boundary. ### Step 2: readings, completeness and a neighbour check Each EA gauge is a tipping-bucket recorder that reports rain every 15 minutes. We take 28 whole days, **6 October to 2 November 2026**, which is 2,688 readings per gauge: `voronoi_thiessen.py` ```python for sid in gauges.index: items = get_json(f"{EA}/id/stations/{sid}/readings?startdate={START}&enddate={END}&_sorted&_limit=10000")["items"] for it in items: if "15_min" in it["measure"] and isinstance(it.get("value"), (int, float)): recs.append((sid, it["dateTime"], it["value"])) # ... totals and completeness per gauge ... gauges = gauges[gauges.complete >= 0.95].copy() ``` A gauge with gaps under-reports the total, so we drop any gauge with less than 95% of the expected readings. Then comes a check that is easy to skip and costly to miss. Rain totals are spatially smooth over a few kilometres, so a gauge that records several times what its neighbours do is suspect. The **Delaunay triangulation** gives each gauge a natural set of neighbours (the gauges whose Voronoi cells touch its cell): `voronoi_thiessen.py` ```python edges0 = shapely.delaunay_triangles(shapely.MultiPoint(list(gauges.geometry)), only_edges=True) # ... map each edge's two ends back to station ids in nbrs ... gauges["nbr_median"] = [gauges.total[list(nbrs[s])].median() for s in gauges.index] gauges["ratio"] = gauges.total / gauges.nbr_median suspect = (gauges.ratio > 2) & (gauges.max15 > 25) gauges = gauges[~suspect].copy() ``` `Output — Step 2` ```text 136,062 readings from 52 gauges (cached in ea_rain_2026-10-06_2026-11-02.csv); max single 15-minute value 57.2 mm, negative values 0 dropped 008850 (NY657958): 83% of the 2,688 expected readings dropped 051729 (NZ008020): 62% of the 2,688 expected readings dropped 047281 (SD829866): 26% of the 2,688 expected readings gauges kept (>=95% complete): 49 (19 inside the catchment) 28-day totals: min 76.6 mm, median 148.4 mm, max 725.5 mm highest totals relative to their Delaunay neighbours: gridref total nbr_median ratio max15 station 602253 NY609428 725.5 115.1 6.3 57.0 613905 NY558774 523.9 112.7 4.6 57.2 013045 NY757386 163.9 90.6 1.8 5.6 592464 NY225135 534.8 323.4 1.7 11.8 588006 SD284986 416.2 256.4 1.6 5.5 flagged 602253 (NY609428): 726 mm against a neighbour median of 115 mm, with 15-minute readings up to 57.0 mm flagged 613905 (NY558774): 524 mm against a neighbour median of 113 mm, with 15-minute readings up to 57.2 mm Thiessen mean with the flagged gauges left in: 191.4 mm; 47 gauges go forward ``` Two gauges stand out: **six and nearly five times** their neighbours’ totals, each with single 15-minute readings of about 57 mm. That is a cloudburst intensity, near the most extreme ever measured in the UK, and none of the neighbours saw anything like it, so they are almost certainly faulty or contain spurious tips. The test needs both conditions. A high ratio alone would also catch real orographic rain: gauge 592464 in the Lake District fells recorded 535 mm, 1.7 times its neighbours, with nothing above 12 mm in any 15 minutes, which is what a wet mountain gauge looks like. ![A map of the Eden catchment as a pale green area outlined in green, with a dashed grey 15 km buffer around it. Blue circles mark 47 rain gauges, labelled with their 28-day totals and shaded from pale blue (about 80 mm) in the east and north to dark blue (over 400 mm) in the south-west near the Lake District fells. Two red crosses labelled 524? and 726? mark the gauges that failed the neighbour check, one near the northern tip and one in the middle of the eastern edge. A red triangle marks the Sheepmount gauging station at the catchment outlet in the north-west.](https://www.technicalgis.com/_astro/01-gauges-and-catchment.IySO9780_NSxPh.webp) _The 47 EA rain gauges kept and their totals for 6 October to 2 November 2026 (mm). The two red crosses failed the neighbour check; the Sheepmount gauging station marks the catchment outlet._ ### Step 3: from Delaunay triangles to Voronoi cells Figure 2 at the top shows the construction. The Delaunay triangulation joins the gauges into triangles whose circumcircles contain no other gauge. The centre of each circumcircle is equidistant from the triangle’s three gauges, so it is a corner of the Voronoi diagram. Joining the circumcentres of triangles that share an edge gives the Voronoi edges, each the perpendicular bisector of a Delaunay edge: `voronoi_thiessen.py` ```python mp = shapely.MultiPoint(list(gauges.geometry)) tri = shapely.delaunay_triangles(mp) # the dual of the Voronoi diagram def circumcentre(poly): (ax_, ay), (bx, by), (cx, cy) = list(poly.exterior.coords)[:3] d = 2 * (ax_ * (by - cy) + bx * (cy - ay) + cx * (ay - by)) ux = ((ax_**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)) / d uy = ((ax_**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax_ - cx) + (cx**2 + cy**2) * (bx - ax_)) / d return ux, uy centres = np.array([circumcentre(p) for p in tri.geoms]) # the Voronoi vertices ``` You never need to do that in practice: GEOS does the whole thing inside `shapely.voronoi_polygons`. The helper used for both parts of this tutorial is four lines: `voronoi_thiessen.py` ```python def thiessen(points, boundary): """Voronoi cells for `points` (GeoSeries), clipped to `boundary`, returned in the same order.""" cells = shapely.voronoi_polygons(shapely.MultiPoint(list(points)), extend_to=boundary, ordered=True) out = gpd.GeoSeries(list(cells.geoms), index=points.index, crs=points.crs) return out.intersection(boundary) gauges["cell"] = thiessen(gauges.geometry, basin) gauges["weight"] = gauges.cell.area / basin.area ``` `extend_to` makes the outer cells reach at least the boundary’s bounding box (by default they stop a little beyond the points), and the intersection then cuts them to the catchment. `Output — Step 3` ```text Delaunay triangulation: 83 triangles from 47 points Voronoi vertices are the triangle circumcentres: 83 of them cells that reach the catchment: 31 of 47 (13 from gauges outside it) clipped cells cover 100.0000% of the catchment ``` The last line is the check to always run: the clipped cells should cover the boundary exactly, with no gaps or overlaps, so the weights sum to 1. ### Step 4: Thiessen weights and catchment rainfall The Thiessen estimate is the area-weighted mean of the gauge totals: each gauge’s total times its share of the catchment. `voronoi_thiessen.py` ```python w = gauges[gauges.weight > 0].sort_values("weight", ascending=False) thiessen_mm = float((w.weight * w.total).sum()) ``` `Output — Step 4 (first 12 of 31 rows)` ```text gridref total area_km2 weight station 604742 NY492308 120.9 203.7 0.090 604142 NY546551 112.7 192.0 0.084 598691 NY721194 88.0 188.3 0.083 605936 NY435360 148.4 163.5 0.072 598928 NY673144 85.4 157.6 0.069 603315 NY642703 130.1 136.4 0.060 601304 NY507159 228.0 135.0 0.059 013336 NY668423 108.0 135.0 0.059 606299 NY389564 94.0 132.8 0.058 598145 NY844121 91.9 129.4 0.057 605543 NY345358 167.3 126.6 0.056 604040 NY575537 114.4 91.0 0.040 ... Thiessen catchment rainfall: 141.1 mm (= sum of weight x gauge total) largest single weight 9.0%, smallest 0.00% ``` ![Two panels. Left: the Eden catchment divided into Thiessen cells shaded by gauge total, darker blue in the south-west, with each larger cell labelled with its share of the area, from 2% to 9%. Black dots mark the gauges, several of them outside the catchment. Right: a horizontal bar chart of the 15 largest weights, from 9.0% for gauge NY492308 (121 mm) down to about 3% for NY426265 (182 mm), with one bar, NY668423, in green to show that the gauge is outside the catchment.](https://www.technicalgis.com/_astro/03-thiessen-weights.6hOVJURq_Z1tcmA4.webp) _a) The Thiessen cells labelled with their area weights; b) the 15 largest weights with each gauge's total. Green marks a gauge that stands outside the catchment but still represents 5.9% of it._ No gauge carries more than 9% of the catchment, which is a healthy spread. The large lowland cells in the middle and east of the catchment, around Penrith and Appleby, carry the most weight, and their totals of 85 to 130 mm pull the estimate well below the median gauge total. ### Step 5: Thiessen against the gauge average and IDW, and sensitivity Two common alternatives are the plain average of the gauges inside the catchment, and inverse distance weighting (IDW) to a grid that is then averaged: `voronoi_thiessen.py` ```python arith = gauges.loc[gauges.inside, "total"].mean() gx, gy = np.meshgrid(xs, ys) # 500 m cell centres inside = shapely.contains_xy(basin, gx, gy) d = np.hypot(gx[inside][:, None] - px[None], gy[inside][:, None] - py[None]) wts = 1 / np.maximum(d, 1) ** 2 idw_mm = ((wts * pv).sum(1) / wts.sum(1)).mean() loo = {} # leave one gauge out for sid in w.index: rest = gauges.drop(sid) loo[sid] = float((thiessen(rest.geometry, basin).area / basin.area * rest.total).sum()) ``` `Output — Step 5` ```text arithmetic mean of the 18 gauges inside: 167.4 mm IDW (power 2, 500 m grid, 9,072 cells, all 47 gauges): 160.3 mm Thiessen: 141.1 mm leave-one-out Thiessen (drop each contributing gauge in turn): 139.7 to 144.7 mm biggest change from dropping 604742 (+3.6 mm) half the gauges (200 random draws of 23): 5th-95th percentile 133.9 to 158.8 mm rain volume over the catchment: 321 million m3 in 28 days ``` ![Two panels. Left: the catchment shaded by an IDW rainfall surface, mostly pale blue with dark blue patches in the south-west where the wettest gauges are, and the Thiessen cell edges drawn over it in thin black lines. Right: a bar chart of catchment rainfall by method: 191 mm for Thiessen before quality control, 167 mm for the arithmetic mean of gauges inside, 160 mm for IDW and 141 mm for Thiessen. Beside the Thiessen bar, a short red error bar shows the leave-one-out range and a longer grey one the 5 to 95% range from using half the gauges.](https://www.technicalgis.com/_astro/04-method-comparison.DX8ZlZbo_ZlPmv.webp) _a) The IDW surface (power 2) with the Thiessen cells over it; b) catchment rainfall for 6 October to 2 November 2026 from the same gauges by each method. The error bars show the leave-one-out range (red) and the 5–95% range when only half the gauges are used (grey)._ The three methods differ by 26 mm, or 19%, and the reason is the network. Six of the 18 gauges inside the catchment are in the wet fells of the south-west, around Ullswater, Haweswater and Shap, with 182 to 346 mm. They are a third of the gauges, so a plain average gives them a third of the say, but their Thiessen cells cover under 18% of the catchment. IDW is pulled the same way, because power 2 lets each wet gauge dominate the cells around it and the wettest gauges outside the catchment, to the south-west, still count. The two flagged gauges matter more than the method: they would have added 50 mm. Dropping any single good gauge moves the Thiessen result by at most 3.6 mm, and even using a random half of the network gives 134 to 159 mm. For a fuller treatment of interpolation, including the uncertainty that Thiessen and IDW cannot give you, see [kriging in Python](https://www.technicalgis.com/tutorials/kriging-spatial-interpolation-python/). ## Part 2: Voronoi vs network service areas for fire stations Voronoi cells are also the quickest way to answer “which facility is nearest to here?”. That is fine for a first look, but people travel by road. Cumbria is a good test: a large rural county with lakes, fells and estuaries in the way. ### Step 6: fire stations from OpenStreetMap, checked `voronoi_thiessen.py` ```python for o in osmium.FileProcessor(str(pbf), osmium.osm.NODE | osmium.osm.WAY).with_locations().with_filter( osmium.filter.TagFilter(("amenity", "fire_station"))): ... # node position, or the mean of a building's nodes fs = fs[fs.within(county)] EXCLUDE = {"n986915522": "former Ulverston station (The Ellers), replaced by the Blue Light Hub in 2021"} airport = fs.fstype == "airport" fs = fs[~fs.osm.isin(EXCLUDE) & ~airport].copy() # a station is often mapped twice (a node and the building outline): merge features within 150 m ``` `Output — Step 6` ```text Cumbria (Cumberland + Westmorland and Furness): 6,825 km2 amenity=fire_station objects in the Geofabrik extract: 41 removed n986915522: former Ulverston station (The Ellers), replaced by the Blue Light Hub in 2021 removed 1 airport fire service (fire_station:type=airport, Walney airfield) names added for 13 features mapped without one inside Cumbria after merging duplicates within 150 m: 38 stations ``` Do not trust a point layer just because it comes from a good source. Before the cleaning, OSM had one closed station still mapped, an airfield fire service that does not answer 999 calls, several stations mapped twice, and 13 without a name. After the cleaning the count matches the fire and rescue service’s own list of **38 stations**: 8 crewed full-time and 30 on call. ### Step 7: Voronoi “nearest station” areas The same helper builds the cells, this time clipped to the county: `voronoi_thiessen.py` ```python stations["cell"] = thiessen(stations.geometry, county) ``` `Output — Step 7` ```text Voronoi cell area: median 155 km2, smallest 13 km2 (Walney), largest 486 km2 (Brampton) ``` ### Step 8: the quickest station by road For the comparison we need drive times from every station to every road junction. We build a drive network for the county from the Geofabrik extract, load it into OSMnx with its default speeds by road type, and run Dijkstra from all 38 stations at once with SciPy, which is much faster than looping in NetworkX: `voronoi_thiessen.py` ```python G = ox.graph_from_xml(osm_file, simplify=True) G = ox.routing.add_edge_travel_times(ox.routing.add_edge_speeds(G)) G = ox.project_graph(ox.truncate.largest_component(G, strongly=True), to_crs=CRS) nodes, edges = ox.graph_to_gdfs(G) # ... build a sparse matrix M of travel times between node indices ... snode = ox.distance.nearest_nodes(G, stations.geometry.x, stations.geometry.y) T = dijkstra(M, directed=True, indices=[idx[n] for n in snode]) + s_snap[:, None] # seconds ``` `Output — Step 8` ```text drive graph: 28,724 nodes, 66,211 edges travel times from 38 stations to 28,724 nodes: 1,091,512 values ``` The routing is directed, from the station outwards, so one-way streets count in the direction an appliance drives. The short hop between each station or home and its nearest node is added at 20 km/h. For the network tools themselves, see [street networks with OSMnx](https://www.technicalgis.com/tutorials/street-networks-osmnx/). ### Step 9: who does the Voronoi send to the wrong station? Census 2021 gives the population of every **output area** (about 300 people each), and the ONS publishes population-weighted centroids for them. For each centroid we compare the station whose Voronoi cell it falls in with the station that is quickest by road: `voronoi_thiessen.py` ```python Toa = T[:, [idx[n] for n in onode]] + o_snap[None] # station x output area, seconds pwc["net_station"] = Toa.argmin(0) pwc["net_min"] = Toa.min(0) / 60 cells = gpd.GeoDataFrame({"vor_station": stations.index}, geometry=stations.cell, crs=CRS) pwc = gpd.sjoin(pwc, cells, predicate="within", how="left") pwc["extra_min"] = Toa[pwc.vor_station.values, np.arange(len(pwc))] / 60 - pwc.net_min mis = pwc[pwc.net_station != pwc.vor_station] ``` `Output — Step 9` ```text 1,742 output areas in Cumbria, population 499,363 output areas whose Voronoi station is not their quickest station by road: 82 (4.7%), 24,090 people (4.8%) extra drive time if sent to the Voronoi station instead: median 1.1 min, 95th percentile 17.1 min, max 18.6 min people with more than 5 extra minutes: 4,680 drive time to the quickest station, population-weighted: median 3.5 min; share within 10 min 97.8%, within 20 min 100.0% Voronoi cell of Carlisle East: 4,521 people are quicker to reach from another station Voronoi cell of Millom: 4,075 people are quicker to reach from another station Voronoi cell of Barrow: 3,529 people are quicker to reach from another station area where the Voronoi cell and the network service area agree: 77.8% of the county ``` ![Two maps of Cumbria. Left: 38 Voronoi cells coloured in four colours so no neighbours match, each with a red triangle and the station name, from Longtown and Brampton in the north to Barrow and Walney in the south and Alston and Kirkby Stephen in the east. Cells are small along the west coast and the south and large in the fells and the north-east. Right: the same county coloured by the quickest station by road for each output area, with the Voronoi edges drawn in black over it. The coloured areas mostly follow the black lines but bulge across them in many places, along valleys and main roads.](https://www.technicalgis.com/_astro/05-voronoi-vs-network.D_clJFDV_Z1dhFy5.webp) _a) Voronoi cells for the 38 fire stations, with four colours so that no two neighbours match; b) the quickest station by road for each 2021 output area (same colours), with the Voronoi edges in black. Where colours cross a black line, the straight-line answer and the road answer differ._ ![A map of Cumbria in pale grey with the Voronoi cell edges and the fire stations. Coloured patches mark the output areas assigned to the wrong station, coloured by the extra drive time: most are small yellow patches of under 2 minutes along cell edges, a few orange ones of 2 to 10 minutes, and one large dark red area of over 10 minutes on the east side of the Duddon estuary in the south-west, opposite Millom.](https://www.technicalgis.com/_astro/06-mismatch.DYLEs_3l_Z1f9dQ2.webp) _Output areas whose Voronoi station is not their quickest station by road, coloured by how much longer the drive from the Voronoi station would be. The dark red area east of the Duddon estuary is in Millom's Voronoi cell, but the estuary lies between them._ For 95% of residents the straight-line answer is right, and most of the mistakes cost a minute or two at the edges of cells where two stations are about equally far. The serious errors are where geography breaks the straight line. The worst is the **Duddon estuary**: Millom’s Voronoi cell reaches across the sands to the Furness side, but a fire engine from Millom has to drive round the head of the estuary, more than 10 minutes further than the trip from Barrow or Ulverston. The same happens on a smaller scale where Arnside’s cell crosses the Kent estuary towards Grange, and where Windermere’s cell reaches across the lake towards Coniston. So use Voronoi cells as a quick nearest-facility map, and for station catchments in open, flat country with a dense road grid. Use network service areas when the answer affects people, especially near water, mountains or motorways with few junctions. If what you want is not “which is nearest?” but “how much service can each place reach, given the demand on it?”, that is a different question: see [spatial accessibility with 2SFCA](https://www.technicalgis.com/tutorials/spatial-accessibility-2sfca-python/), which uses the same OSMnx drive times, and [location-allocation in Python](https://www.technicalgis.com/tutorials/location-allocation-python/) for choosing where facilities should go. ## Limitations and honest caveats - **Short, unchecked rain records.** The EA real-time data are provisional and only about four weeks are online. Use the quality-checked archive in the Hydrology Data Explorer for any formal water-balance work, and longer periods than 28 days. - **The neighbour check is a screen, not a verdict.** The thresholds (twice the neighbour median and a 15-minute reading over 25 mm) were set for this autumn period; a summer thunderstorm can legitimately give one gauge far more than its neighbours. Look at the hyetograph of any flagged gauge before you drop it. - **Thiessen ignores terrain.** See the note in Step 5: in mountain catchments, an elevation-aware gridded product is usually better than any point weighting. - **Drive time is not response time.** On-call crews have to reach the station first, typically several minutes, and the OSMnx speeds are free-flow defaults by road type, not observed speeds. - **Edge effects.** Stations in Dumfries and Galloway, Northumberland, Durham, North Yorkshire and Lancashire are left out, so cells and service areas at the county edge are bigger than in reality; cross-border cover is common. - **OSM needs checking.** The station list here was checked by hand against the service’s published list in November 2026. Repeat the check if you rerun it later. ## Frequently asked questions ### How do I make Voronoi polygons in Python with GeoPandas? Project the points to a CRS in metres, then call `shapely.voronoi_polygons(shapely.MultiPoint(list(gdf.geometry)), extend_to=boundary, ordered=True)` and put `list(cells.geoms)` back into the GeoDataFrame as a new geometry column; with `ordered=True` (shapely 2.1+) the cells come back in the same order as the points. `gdf.voronoi_polygons()` in GeoPandas 1.0+ also works, but its output order does not match the input, so join the cells back to the points with `gpd.sjoin(cells, points, predicate="contains")`. ### What is the difference between Voronoi and Thiessen polygons? There is none in the geometry: they are the same cells, where every location is closer to its own point than to any other. “Voronoi” is the mathematical name, after Georgy Voronoi; “Thiessen polygons” is the name used in hydrology and geography, after Alfred Thiessen, who used them in 1911 to weight rain gauges for areal rainfall. “Dirichlet tessellation” is a third name for the same thing. ### How do I clip Voronoi polygons to a boundary? Pass the boundary as `extend_to` so the outer cells are long enough to reach it, then intersect: `GeoSeries(cells).intersection(boundary)`. Check that the clipped areas add up to the boundary’s area. Include the points just outside the boundary when you build the diagram, otherwise the edge cells inside are too big. ### How do I calculate Thiessen-weighted rainfall for a catchment? Build the Voronoi cells for all gauges in and around the catchment, clip them to it, and divide each cell’s area by the catchment area to get the weights, which add up to 1. The catchment rainfall is the sum of weight × gauge total. Drop gauges with missing readings and check for faulty ones first; here two bad gauges raised the result from 141 to 191 mm. ### Why do my Voronoi polygons look wrong or stretched? Usually because the points are in longitude and latitude. Voronoi cells use straight-line distance in the coordinate units, and a degree of longitude is much shorter than a degree of latitude in the UK, so the cells are squashed. Reproject to a projected CRS such as EPSG:27700 first. Duplicate points can also produce missing cells; drop them with `drop_duplicates` on the geometry. ### Is a Voronoi diagram the same as a service area? No. A Voronoi cell is the area closest in a straight line; a service area is the area quickest to reach along a network. They agree where roads are dense and direct, and disagree around estuaries, lakes, hills and one-way systems. In Cumbria they matched for 77.8% of the county’s area and 95.2% of its residents. ## Next steps - Measure access, not just nearness, with [spatial accessibility with 2SFCA](https://www.technicalgis.com/tutorials/spatial-accessibility-2sfca-python/), and choose new facility sites with [location-allocation](https://www.technicalgis.com/tutorials/location-allocation-python/). - Draw a catchment from a DEM instead of downloading it with [watershed delineation in Python](https://www.technicalgis.com/tutorials/watershed-delineation-python/), then weight the gauges inside it. - Replace Thiessen with a statistical surface and its uncertainty in [kriging in Python](https://www.technicalgis.com/tutorials/kriging-spatial-interpolation-python/), or average a gridded rainfall product over the catchment with [zonal statistics](https://www.technicalgis.com/tutorials/zonal-statistics-python/). - Test whether the gauges or stations are clustered or regular with [point pattern analysis](https://www.technicalgis.com/tutorials/point-pattern-analysis-python/), or swap Voronoi cells for a regular grid with [H3 hexagons](https://www.technicalgis.com/tutorials/h3-hexagon-grid-python/). ## Data and licences - **Rainfall and gauges**: Environment Agency real-time flood-monitoring API, readings for 6 October to 2 November 2026, downloaded on 3 November 2026. This uses Environment Agency flood and river level data from the real-time data API (Beta). Licensed under the [Open Government Licence v3.0](https://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/). - **Catchment boundary**: NRFA station 76007, Eden at Sheepmount, from the [National River Flow Archive](https://nrfa.ceh.ac.uk/data/station/info/76007) WFS. © UK Centre for Ecology & Hydrology; used under the NRFA’s terms, which require this attribution. - **Fire stations and roads**: © OpenStreetMap contributors, available under the [Open Database Licence](https://www.openstreetmap.org/copyright). Cumbria extract from [Geofabrik](https://download.geofabrik.de/europe/united-kingdom/england/cumbria.html), downloaded on 3 November 2026. Station list checked against [Cumbria Fire and Rescue Service](https://www.cumbriafire.gov.uk/our-fire-stations). - **Boundaries and centroids**: Office for National Statistics, Local Authority Districts (May 2026) and Output Areas (2021) boundaries and population-weighted centroids. Contains OS data © Crown copyright and database right 2026. Source: Office for National Statistics licensed under the Open Government Licence v3.0. - **Population**: Census 2021 table TS001 (usual residents) by output area, from [Nomis](https://www.nomisweb.co.uk/). Office for National Statistics, Open Government Licence v3.0. - Code: [`voronoi_thiessen.py`](https://www.technicalgis.com/code/voronoi-thiessen-polygons-python/), tested with Python 3.12.0, geopandas 1.2.0, shapely 2.2.0, pyogrio 0.13.0, pyproj 3.8.0, numpy 2.5.4, pandas 3.1.0, scipy 1.18.1, networkx 3.7, osmnx 2.1.1, osmium 4.3.2, matplotlib 3.11.2 and requests 2.34.2. ## Key terms - **Voronoi diagram**: A division of the plane into cells, one per input point, where every location in a cell is closer to that cell's point than to any other point. Each cell edge lies halfway between two neighbouring points. - **Thiessen polygons**: The name used in hydrology and geography for Voronoi cells, after the meteorologist Alfred Thiessen, who used them in 1911 to weight rain gauges by the area each one represents. - **Delaunay triangulation**: The triangulation of a set of points in which no point lies inside any triangle's circumcircle. It is the dual of the Voronoi diagram: two points share a Voronoi edge exactly when they are joined by a Delaunay edge. - **Circumcentre**: The centre of the circle through a triangle's three corners. The circumcentres of the Delaunay triangles are the corners (vertices) of the Voronoi cells. - **Areal rainfall**: The average depth of rain over an area such as a catchment, in mm. It is estimated from point gauges by a weighting scheme such as Thiessen polygons, inverse distance weighting or kriging. - **Service area**: The area that is quickest to reach from (or to) a facility along a network, such as roads. Unlike a Voronoi cell, it follows the network, so rivers, estuaries and hills bend its edges. ## Data for this tutorial - [National River Flow Archive (NRFA)](https://www.technicalgis.com/data-sources/united-kingdom/#src-nrfa) · Catchment boundary and area for station 76007, Eden at Sheepmount, from the NRFA WFS - [Environment Agency real-time flood-monitoring API](https://www.technicalgis.com/data-sources/united-kingdom/#src-ea-flood-monitoring-api) · 15-minute tipping-bucket rainfall from 52 gauges, 6 October to 2 November 2026 - [Geofabrik OpenStreetMap extract (England)](https://www.technicalgis.com/data-sources/united-kingdom/#src-geofabrik-openstreetmap-extract-england) · Cumbria .osm.pbf: amenity=fire\_station and the drivable road network - [ONS Open Geography Portal](https://www.technicalgis.com/data-sources/united-kingdom/#src-open-geography-portal) · Local authority boundaries (May 2026), 2021 output area boundaries and population-weighted centroids - [Nomis](https://www.technicalgis.com/data-sources/united-kingdom/#src-nomis) · Census 2021 TS001, usual residents by output area Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/voronoi-thiessen-polygons-python/