Polygon tessellation is the stage that turns a set of synthetic seed points into a complete, non-overlapping partition of geographic space — the zones, catchments, and administrative cells that downstream models aggregate over. This page is part of Spatial Distribution & Pattern Generation: where point process simulation models emit the discrete seeds and density mapping smooths them into intensity surfaces, tessellation is the step that converts seeds into bounded vector regions. Done in a notebook it is a one-line scipy.spatial.Voronoi call; done inside a generation pipeline it has to be deterministic, topologically valid, area-true, and reproducible across machines — or every feature engineered from the resulting zones quietly drifts.
The engineering failure this page solves is the leaky partition — a tessellation that looks like a clean set of zones on screen but is not actually a gap-free, non-overlapping cover of the extent. The symptoms are familiar to anyone who has aggregated synthetic events into generated zones and watched the totals not add up:
Sliver polygons. Floating-point imprecision and aggressive boundary snapping leave hair-thin polygons between intended neighbours. They carry a handful of events each, inflate the zone count, and break any join keyed on “one event, one zone.”
Gaps and overlaps at constraint edges. When a tessellation is clipped to a coastline or administrative boundary with naive coordinate comparison rather than robust predicates, edges cross, cells overlap, and the union of all zones no longer equals the input extent.
Unbounded outer cells. Raw Voronoi diagrams produce infinite cells around the convex hull. Forget to clip them and the outermost zones have no finite area, so per-area normalization silently divides by infinity or NaN.
Non-reproducible topology. Let the seed RNG draw from wall-clock time, or feed an unstable row order into the Delaunay step, and two runs on the same configuration produce different adjacency graphs. Anything downstream that assumed a stable zone-to-neighbour map is now wrong.
Each defect is individually small and invisible in a one-off plot. Composed inside a continuous pipeline they produce zones that “look about right” but fail area-conservation checks, break statistical-parity gates, and corrupt every aggregation. The remedy is the same separation-of-concerns discipline that governs CRS contract enforcement elsewhere in the stack: make the projection, the seed RNG, the tolerance, and the topology gate explicit, pinned, and tested.
Tessellation sits on the standard synthetic-spatial toolchain. Pin majors explicitly so geometric predicates and snapping tolerances do not shift under you between environments:
Two environment decisions matter more than the versions. First, never tessellate in geographic degrees. Computing areas, distances, or snapping tolerances in EPSG:4326 makes every threshold latitude-dependent — a 1e-7 tolerance is centimetres at the equator and metres near the poles. Reproject to a metric, equal-area or local CRS (a UTM zone regionally, EPSG:6933 globally) before any geometry runs, exactly as the realism metrics evaluation stage expects metric inputs. Second, pin GEOS through Shapely 2.x: make_valid and the overlay predicates changed behaviour across GEOS minor versions, and an unpinned environment will repair geometries differently on CI than on a developer laptop.
Every tessellation in this family is a function of three things — a set of seeds, a partition rule, and a set of constraints — and getting the partition right means controlling all three deliberately.
Voronoi diagrams assign every location to its nearest seed, producing convex cells whose boundaries are the perpendicular bisectors between adjacent seeds. They are the natural model for catchment areas, service regions, and synthetic administrative zoning, because “nearest facility” is exactly the decision rule a Voronoi cell encodes. The drill-down on optimizing Voronoi tessellation for synthetic zoning maps covers the half-edge data structures and incremental insertion that keep this step memory-bounded at scale.
Delaunay triangulation is the dual graph: connect two seeds whenever their Voronoi cells share an edge and you get the triangulation that maximizes the minimum angle, avoiding the sliver triangles that wreck interpolation. It is rarely the final output, but it is the workhorse intermediate — for mesh refinement, for barycentric interpolation surfaces, and as the substrate for constrained tessellation, where you force specific edges (a river, a road, a jurisdictional line) to appear in the mesh so cells never straddle a real-world barrier.
Grid tessellations — hexagonal, quadtree, or adaptive irregular grids — trade the data-driven shape of Voronoi cells for uniform area and computational regularity. Choose them when the consumer needs comparable cell areas (binning for a feature store) rather than meaningful zone boundaries.
The decisive lever is where the seeds come from. Drawing seeds from a uniform distribution yields artificially regular cells that no real geography exhibits. Drawing them from a Poisson, Thomas, or Matérn cluster process — the same generators described under point process simulation models — makes the resulting partition inherit realistic spatial clustering, so the zones preserve the second-order statistics that a graph neural network or spatial-aggregation feature later depends on. Tessellation does not add spatial structure; it inherits it from the seed process, which is why the seed step is part of this pipeline and not an afterthought.
For privacy-sensitive zoning, the seed-weighting step is also where you apply protection. Perturbing seed positions or per-seed weights with calibrated noise under a differential privacy budget ensures the generated zones cannot be inverted back to a real cadastral or demographic boundary, while the partition as a whole still reproduces the prescribed intensity.
The pipeline runs as four auditable stages. Each block below is copy-pasteable and assumes the pinned toolchain above.
Step 1 — Ingest the boundary and normalize the CRS. Load the constraint geometry, reproject to a metric CRS, and confirm it is a valid planar polygon before anything else touches it.
python
import geopandas as gpd
from shapely.validation import make_valid
METRIC_CRS ="EPSG:6933"# equal-area; swap for a local UTM zone regionally
SNAP_TOL_M =1e-6# snapping tolerance in metres, NOT degreesdefload_boundary(path:str)-> gpd.GeoDataFrame:
boundary = gpd.read_file(path).to_crs(METRIC_CRS)
boundary["geometry"]= boundary.geometry.apply(make_valid)assert boundary.geometry.is_valid.all(),"constraint geometry is not planar-valid"return boundary
Step 2 — Generate and weight seeds deterministically. Seed the RNG from an explicit, version-controlled integer so the same configuration always yields the same point set. Here a Poisson process fills the boundary’s bounding box and is clipped to the boundary.
python
import numpy as np
from shapely.geometry import MultiPoint
defpoisson_seeds(boundary: gpd.GeoDataFrame, intensity:float, rng_seed:int)-> np.ndarray:
rng = np.random.default_rng(rng_seed)# PCG64; document the algorithm
minx, miny, maxx, maxy = boundary.total_bounds
area =(maxx - minx)*(maxy - miny)
n = rng.poisson(intensity * area)# count is itself stochastic
pts = np.column_stack([
rng.uniform(minx, maxx, n),
rng.uniform(miny, maxy, n),])
inside = gpd.GeoSeries(gpd.points_from_xy(pts[:,0], pts[:,1]), crs=METRIC_CRS)
mask = inside.within(boundary.union_all())return pts[mask.values]
Step 3 — Tessellate with spatial-index acceleration. Shapely 2.x exposes voronoi_polygons, which delegates to GEOS and returns finite, clipped cells when an envelope is supplied — sidestepping the unbounded-outer-cell trap of raw scipy.spatial.Voronoi.
python
from shapely import voronoi_polygons
from shapely.geometry import MultiPoint
deftessellate(seeds: np.ndarray, boundary: gpd.GeoDataFrame)-> gpd.GeoDataFrame:
envelope = boundary.union_all()
cells = voronoi_polygons(
MultiPoint(seeds),
extend_to=envelope,# clip outer cells to the constraint, no infinities
only_edges=False,)
clipped = gpd.GeoSeries(list(cells.geoms), crs=METRIC_CRS).intersection(envelope)
out = gpd.GeoDataFrame(geometry=clipped[~clipped.is_empty], crs=METRIC_CRS)return out.reset_index(drop=True)
For very large seed counts, build an R-tree over the constraint boundary first (gpd.GeoDataFrame.sindex) so the per-cell clip is a bounded-box query rather than an O(n) scan against the full boundary.
Step 4 — Repair topology and snap to a gap-free cover. Dissolve slivers below an area floor into their largest neighbour, then assert the union equals the input extent within tolerance.
python
defrepair(cells: gpd.GeoDataFrame, min_area_m2:float)-> gpd.GeoDataFrame:
cells["geometry"]= cells.geometry.apply(make_valid)
cells = cells.set_precision(SNAP_TOL_M)# GEOS precision-grid snap
slivers = cells[cells.area < min_area_m2]iflen(slivers):# merge each sliver into the neighbour it shares the longest edge with
cells = _dissolve_slivers(cells, slivers, min_area_m2)return cells.reset_index(drop=True)
The set_precision snap is what makes neighbouring edges share identical coordinates rather than merely-close ones — the difference between a true partition and a pile of nearly-adjacent polygons.
A tessellation is only trustworthy if a CI gate proves it is a complete, non-overlapping, reproducible cover. Four assertions catch the failure modes above.
python
import numpy as np
defgate_partition(cells: gpd.GeoDataFrame, boundary: gpd.GeoDataFrame,
area_tol:float=1e-3)->None:# 1. Every cell is geometrically valid.assert cells.geometry.is_valid.all(),"invalid geometry survived repair"# 2. Area is conserved: union of cells equals the input extent.
covered = cells.union_all().area
extent = boundary.union_all().area
assertabs(covered - extent)/ extent < area_tol,"partition leaks area"# 3. No overlaps: summed cell area equals the covered area.assertabs(cells.area.sum()- covered)/ covered < area_tol,"cells overlap"# 4. No empty or zero-area cells slipped through.assert(cells.area >0).all(),"zero-area cell present"
Reproducibility gets its own test: run Step 2 through Step 4 twice with the same rng_seed and assert the geometries are byte-identical via well-known binary.
python
deftest_tessellation_is_deterministic(boundary):
a = repair(tessellate(poisson_seeds(boundary,1e-6,42), boundary),50.0)
b = repair(tessellate(poisson_seeds(boundary,1e-6,42), boundary),50.0)assertlist(a.geometry.apply(lambda g: g.wkb))== \
list(b.geometry.apply(lambda g: g.wkb)),"topology not reproducible"
Beyond geometric correctness, gate on statistical fidelity: the distribution of cell areas and the seed nearest-neighbour function should match the prescribed seed process. Comparing those summaries against the baseline prior — the kind of distance-based check the realism metrics evaluation stage formalizes — is what separates “the polygons are valid” from “the polygons are right.” Wire all of these as build-failing checks through the CI/CD integration pipeline so an invalid partition can never be promoted as an artifact.
Continental and global extents push tessellation into memory and latency territory that a single-pass, in-memory approach cannot survive.
Spatial chunking with buffered halos. Partition the extent into tiles, generate seeds and tessellate each tile independently, but extend every tile by a halo of at least one expected cell diameter before tessellating. Cells near a seam depend on seeds across the seam; without the halo you get visible ridges of malformed cells exactly on the tile boundaries. Trim the halo before stitching, keeping only cells whose seed falls inside the un-buffered tile.
Deterministic per-tile seeding. Derive each tile’s RNG seed from its integer (i, j) coordinates and the pipeline version hash, never from wall-clock time or process ID. A tile then produces identical geometry regardless of when or in what order it runs — which is what lets the async execution for large grids layer tessellate tiles concurrently without sacrificing reproducibility.
Streaming, not accumulation. Write each finished tile to a partitioned Parquet or FlatGeobuf sink and release it, rather than concatenating every cell in a single GeoDataFrame. Out-of-core spatial indexes and lazy evaluation keep peak heap flat as seed density climbs.
Bounded predicates over brute force. Replace O(n²) nearest-seed scans with an R-tree query, and prefer Shapely 2.x’s vectorized intersection/union_all over Python-level loops — the GEOS-backed batch ops are an order of magnitude faster and avoid the per-geometry object overhead that dominates large runs.
Why do thin sliver polygons keep appearing between adjacent cells?
Slivers come from floating-point imprecision: two cells that should share an edge instead have edges that differ in the last few decimal places, leaving a hair-thin gap polygon between them. Snap every geometry to a shared precision grid with Shapely's set_precision before merging, and dissolve any residual polygon below an area floor (calibrated to your spatial resolution) into the neighbour it shares the longest edge with. Snapping to a common grid is the fix; an area floor alone only hides the symptom.
My outermost zones have infinite or NaN area — what happened?
Raw Voronoi diagrams produce unbounded cells around the convex hull, so the outer zones have no finite extent. Never export the raw diagram. Clip every cell to the constraint envelope — Shapely's voronoi_polygons with extend_to set to the boundary returns finite cells directly, and a final intersection against the boundary guarantees no cell escapes the extent. After clipping, assert that every cell area is strictly positive and finite.
Two runs with the same seed produce different adjacency graphs. Why?
Non-reproducible topology almost always traces to an unstable input. Either the RNG drew from wall-clock time instead of an explicit integer seed, or the seed array reached the tessellator in a different row order between runs. Fix both: initialize numpy's default_rng with a version-controlled integer, sort or otherwise stabilize the seed order before tessellating, and add a CI test that runs the pipeline twice and asserts byte-identical well-known binary.
Cells cross a coastline or administrative boundary they should respect.
An unconstrained Voronoi diagram knows nothing about barriers, so cells span them freely. Either run a constrained Delaunay triangulation that forces the barrier edges into the mesh, or — simpler for clipping — intersect the finished tessellation with the multipart constraint geometry so each cell is split at the barrier. Use robust GEOS predicates for the overlay; naive coordinate comparison fails exactly at the boundary vertices where it matters most.
The partition passes geometry checks but downstream totals don't add up.
Geometric validity is necessary but not sufficient. If summed cell areas don't equal the covered extent, you have overlaps; if the covered extent is smaller than the input boundary, you have gaps. Add an area-conservation gate that compares the union of cells against the input extent within a relative tolerance, plus a no-overlap gate comparing the summed cell area against the union area. Both must pass before the partition is promoted, because any leak corrupts every event-to-zone aggregation.
