Create Voronoi Polygons and Delaunay Triangles in PyQGIS
Given a set of points, two structures describe how they share space. Voronoi polygons — also called Thiessen polygons — give every point the region of the plane closer to it than to any other point. The Delaunay triangulation connects points whose Voronoi cells touch, forming triangles that are as close to equilateral as the points allow. They are two views of the same geometry, and between them they answer a surprising range of practical questions: which weather station's reading applies to this field, which school is nearest each address, which points are neighbours of which.
This recipe belongs to Spatial Statistics & Pattern Analysis. It builds both structures with Processing, clips Voronoi cells to a meaningful area, uses them for proximity allocation, and turns the triangulation into a list of neighbours.
Prerequisites
- QGIS 3.34 LTR or newer, or the QGIS 4 series.
- A point layer in a projected CRS. Voronoi cells drawn in degrees are distorted, especially at high latitudes.
- No duplicate points. Two points at the same location have no meaningful boundary between them; remove duplicates first.
Build Voronoi polygons
The Voronoi algorithm takes a point layer and a buffer percentage that controls how far the outer cells extend beyond the points' extent. Attributes of each point are copied to its cell.
import processing
from qgis.core import QgsProject
stations = QgsProject.instance().mapLayersByName("rain_gauges")[0]
cells = processing.run("native:voronoipolygons", {
"INPUT": stations,
"BUFFER": 10, # percent of extent added around the points
"TOLERANCE": 0,
"COPY_ATTRIBUTES": True,
"OUTPUT": "memory:gauge_cells",
})["OUTPUT"]
print(cells.featureCount(), "cells for", stations.featureCount(), "gauges")
QgsProject.instance().addMapLayer(cells)
Breakdown: Each output polygon carries the attributes of the point it surrounds, so a cell can be styled or joined by the gauge's id or reading directly. Cells along the outside of the point set are unbounded in theory; the buffer clips them to a rectangle slightly larger than the points' extent. A count mismatch between cells and points usually means duplicate points. TOLERANCE snaps nearly coincident points before building and can be left at zero for clean data; on QGIS releases before the algorithm moved to the native provider, the id is qgis:voronoipolygons and the tolerance parameter is absent.
Clip cells to a meaningful area
The rectangle the buffer produces is arbitrary. For analysis, cells should stop at the boundary of the area they describe — a catchment, a country, a coastline.
basin = QgsProject.instance().mapLayersByName("river_basin")[0]
clipped = processing.run("native:clip", {
"INPUT": cells, "OVERLAY": basin, "OUTPUT": "memory:gauge_cells_basin"})["OUTPUT"]
areas = processing.run("native:fieldcalculator", {
"INPUT": clipped, "FIELD_NAME": "cell_km2", "FIELD_TYPE": 0,
"FIELD_LENGTH": 12, "FIELD_PRECISION": 3,
"FORMULA": "$area / 1e6", "OUTPUT": "memory:gauge_cells_final"})["OUTPUT"]
total = sum(f["cell_km2"] for f in areas.getFeatures())
print(f"basin area covered: {total:.1f} km²")
QgsProject.instance().addMapLayer(areas)
Breakdown: Clipping keeps attributes, so the cells remain linked to their gauges. Area must be recalculated after clipping, because a cell that extended into the sea is now smaller; $area uses the project's ellipsoid settings, while area($geometry) gives planar area in layer units. The sum of cell areas should equal the basin area — a quick check that no cell was lost. A gauge outside the basin may still own a cell that reaches inside it; that is correct, as the nearest gauge to those locations is outside.
Weight values by Thiessen area
The classic use of Voronoi cells is the Thiessen method for areal averages: each gauge's rainfall is weighted by the area of its cell within the basin, giving a basin-average rainfall that does not over-weight clusters of gauges.
weighted = sum(f["rain_mm"] * f["cell_km2"] for f in areas.getFeatures())
simple = sum(f["rain_mm"] for f in areas.getFeatures()) / areas.featureCount()
print(f"Thiessen mean: {weighted / total:.1f} mm | simple mean: {simple:.1f} mm")
Breakdown: Where gauges are clustered — several in a city, one on a mountain — a simple mean lets the cluster dominate. The Thiessen weights reflect how much of the basin each gauge represents. The difference between the two means is itself informative: a large gap shows that the network is unevenly spread and that the simple mean would mislead. Interpolation methods such as IDW give a continuous surface instead of constant cells.
Allocate locations to their nearest site
Voronoi cells are a fast way to allocate many locations to the nearest of a few sites by straight-line distance: build cells around the sites once, then a point-in-polygon join does the rest.
schools = QgsProject.instance().mapLayersByName("primary_schools")[0]
addresses = QgsProject.instance().mapLayersByName("addresses")[0]
school_cells = processing.run("native:voronoipolygons", {
"INPUT": schools, "BUFFER": 20, "COPY_ATTRIBUTES": True,
"OUTPUT": "memory:"})["OUTPUT"]
allocated = processing.run("native:joinattributesbylocation", {
"INPUT": addresses, "JOIN": school_cells, "PREDICATE": [5], # within
"JOIN_FIELDS": ["school_name"], "METHOD": 1,
"OUTPUT": "memory:addresses_nearest_school"})["OUTPUT"]
QgsProject.instance().addMapLayer(allocated)
Breakdown: For thousands of addresses and dozens of sites, building cells once and joining is much faster than a nearest-neighbour search per address, and the cells themselves are a useful map of the allocation. The limitation is the distance measure: cells use straight-line distance, so a river without a bridge or a motorway without a crossing is ignored. For allocation by travel time, use service areas or an origin–destination matrix instead.
Measure how catchments cross administrative lines
Once sites have cells, overlaying them with administrative units shows how a site's catchment is split — useful for funding formulas, reporting by district, or spotting a school whose nearest-pupils area straddles a council boundary.
districts = QgsProject.instance().mapLayersByName("districts")[0]
pieces = processing.run("native:intersection", {
"INPUT": school_cells, "OVERLAY": districts,
"INPUT_FIELDS": ["school_name"], "OVERLAY_FIELDS": ["district"],
"OUTPUT": "memory:"})["OUTPUT"]
shares = defaultdict(dict)
for f in pieces.getFeatures():
shares[f["school_name"]][f["district"]] = f.geometry().area()
for school, parts in shares.items():
total = sum(parts.values())
if len(parts) > 1:
split = ", ".join(f"{d} {a / total:.0%}" for d, a in sorted(parts.items(), key=lambda kv: -kv[1]))
print(f"{school}: {split}")
Breakdown: The intersection splits each cell along district boundaries and keeps one attribute from each side, so every piece knows its school and its district. Summing piece areas per school gives the share of the catchment in each district; printing only schools with more than one district focuses on the cross-boundary cases. Area shares are a proxy — weighting by population, by intersecting with an address layer instead of measuring area, gives the share of people rather than land, which is usually what matters. Add from collections import defaultdict if you run this section on its own.
Build a Delaunay triangulation
The triangulation algorithm connects the points into triangles. Its output is a polygon layer of triangles; the edges of those triangles are the neighbour relationships.
triangles = processing.run("native:delaunaytriangulation", {
"INPUT": stations, "TOLERANCE": 0, "ADD_ATTRIBUTES": True,
"OUTPUT": "memory:gauge_tin"})["OUTPUT"]
edges = processing.run("native:polygonstolines", {
"INPUT": triangles, "OUTPUT": "memory:"})["OUTPUT"]
exploded = processing.run("native:explodelines", {
"INPUT": edges, "OUTPUT": "memory:"})["OUTPUT"]
unique = processing.run("native:deleteduplicategeometries", {
"INPUT": exploded, "OUTPUT": "memory:tin_edges"})["OUTPUT"]
lengths = sorted(f.geometry().length() for f in unique.getFeatures())
print(len(lengths), "edges; median length", round(lengths[len(lengths) // 2]))
Breakdown: Converting triangles to lines and exploding them gives one segment per triangle side; each interior edge appears twice (once per adjacent triangle), so deleting duplicate geometries leaves one per neighbour pair. With ADD_ATTRIBUTES, each triangle records the ids of its three corner points on recent releases, which makes building a neighbour list by id straightforward. The edge-length distribution is a quick diagnostic: long edges at the outside of the point set span empty space and are usually excluded when the triangulation is used for neighbour relations.
Turn edges into a neighbour list
Spatial statistics such as the Getis-Ord hot spot test need to know which features are neighbours. Delaunay edges, with overly long ones removed, are a natural definition.
from collections import defaultdict
from qgis.core import QgsSpatialIndex
MAX_EDGE = lengths[int(len(lengths) * 0.95)] # drop the longest 5 %
index = QgsSpatialIndex(stations.getFeatures())
pts = {f.id(): f.geometry() for f in stations.getFeatures()}
neighbours = defaultdict(set)
for e in unique.getFeatures():
g = e.geometry()
if g.length() > MAX_EDGE:
continue
line = g.asPolyline()
a = index.nearestNeighbor(line[0], 1)[0]
b = index.nearestNeighbor(line[-1], 1)[0]
neighbours[a].add(b)
neighbours[b].add(a)
print(sum(len(v) for v in neighbours.values()) / max(len(neighbours), 1), "neighbours on average")
Breakdown: Each edge's end points coincide with two input points, so a nearest-neighbour lookup recovers their ids. Dropping the longest edges stops points on opposite sides of a gap — across a lake, between two towns — from being treated as neighbours. Delaunay neighbours average about six per point, a well-behaved number for spatial weights.
QGIS version compatibility
native:voronoipolygons and native:delaunaytriangulation are the GEOS-based implementations used on QGIS 3.34 LTR, 3.40 LTR and QGIS 4; earlier 3.x releases used qgis:voronoipolygons and qgis:delaunaytriangulation with slightly different parameters. Check with processing.algorithmHelp(...) on your version. native:clip, native:fieldcalculator and the line algorithms are stable across all.
Troubleshooting
- Fewer cells than points. Duplicate or near-coincident points; remove duplicates or set a tolerance.
- Cells look stretched. The layer is in a geographic CRS; reproject first.
- Allocation ignores barriers. Voronoi uses straight-line distance; use network analysis instead.
- Neighbour lists link distant points. Long outer edges were kept; filter by length.
Conclusion
Build Voronoi cells with native:voronoipolygons, clip them to the real study area and recompute areas, use them for Thiessen weighting and fast nearest-site allocation, build Delaunay triangles for neighbour relationships and surfaces, and filter long edges before treating them as neighbours.
Frequently Asked Questions
Are Thiessen and Voronoi polygons the same? Yes. "Thiessen polygons" is the name used in hydrology and climatology for the same construction.
Can I weight Voronoi cells by site capacity?
Not with this algorithm; weighted (power) diagrams need a library such as scipy or a custom implementation.
How do I get a TIN surface from the triangles? Use TIN interpolation, which builds and rasterises the triangulation in one step.
Do the cells change if I add a point? Only the cells around the new point change, but the algorithm rebuilds everything; rerun it.