When a single scipy.spatial.Voronoi call over a metropolitan-scale seed set exhausts heap memory and emits unclosed, self-intersecting zoning cells, the fix is to replace monolithic tessellation with bounded incremental construction over a half-edge mesh.
Part of Polygon Tessellation Algorithms: that page covers the end-to-end seed-to-partition workflow; this page drills into the specific failure where tessellation scales correctly in a notebook but collapses — on memory, on precision, or on topology — the moment the seed count and extent reach production size.
Three independent mechanisms converge into one failure chain, and all three are deterministic rather than random.
First, unbounded cells. A raw Voronoi diagram produces infinite rays around every seed on the convex hull. Clipping those rays to an irregular zoning envelope is where the geometry library is asked to intersect a ray with a polygon edge — and when seed coordinates come straight off a uniform grid or an unjittered spatial Poisson point process, three seeds frequently end up collinear. Collinear seeds yield degenerate circumcircles with effectively infinite radius, so the Delaunay dual that underpins the diagram becomes ill-conditioned and the clipped ring fails to close.
Second, memory fragmentation. The naive pattern builds the full triangulation, then the full diagram, then clips every cell, holding all three representations in memory at once. Above roughly 10⁵ seeds the intermediate edge lists, vertex arrays, and per-cell coordinate buffers fragment the heap; peak resident memory runs several multiples of the final output size, and the process is killed long before it writes a partition.
Third, silent topological corruption. The unclosed rings and duplicate vertices produced by the first two mechanisms do not always raise. They serialize to GeoJSON as plausible-looking polygons and propagate into downstream Spatial Distribution & Pattern Generation consumers, where they corrupt density estimates and invalidate any feature engineered from zone areas. Because synthetic zoning frequently feeds privacy-sensitive aggregation, a leaky partition also breaks the area-conservation invariant that the differential privacy mechanisms downstream depend on for their re-identification bounds.
The optimization, therefore, is not a faster diagram — it is a construction strategy that never materializes the whole diagram at once and never lets a degenerate circumcircle form in the first place.
Reproduce both failure modes deterministically. Seeds drawn from an unjittered grid guarantee collinearity; a large count drives the memory profile.
python
import numpy as np
import tracemalloc
from scipy.spatial import Voronoi
# Unjittered grid -> guaranteed collinear seeds -> degenerate circumcircles.
side =360# 360 x 360 = 129,600 seeds
xs, ys = np.meshgrid(np.linspace(0,50_000, side),
np.linspace(0,50_000, side))
seeds = np.column_stack([xs.ravel(), ys.ravel()])
tracemalloc.start()
vor = Voronoi(seeds)# builds the FULL diagram in RAM
_, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()# Count the unbounded cells you would have to clip by hand.
unbounded =sum(1for region in vor.regions if-1in region ornot region)print(f"peak RAM: {peak /1e6:.0f} MB, unbounded regions: {unbounded}")# peak RAM: ~900 MB, unbounded regions: 1436 -> scales super-linearly with `side`
The unbounded-region count is non-zero by construction, and peak memory grows faster than the seed count because the intermediate edge structures dominate. Both numbers get worse, not better, as the extent approaches a real metropolitan footprint.
The production fix has two parts: stabilize the seeds so no degenerate circumcircle can form, then build the partition tile-by-tile against a fixed memory ceiling using a half-edge (doubly-connected edge list) layout that discards triangulation scratch as soon as each cell’s boundary is extracted.
Apply a hash-derived jitter to break collinearity without introducing statistical bias, and enforce a minimum inter-seed distance tied to the target cell radius. This mirrors the seeding contract used by the point process simulation models so the same seeds reproduce across machines.
python
import numpy as np
defstabilize_seeds(seeds: np.ndarray,*, seed:int, jitter:float)-> np.ndarray:"""Break collinearity deterministically; jitter is small vs. cell radius."""
rng = np.random.default_rng(seed)# version-controlled integer seed# Hash-based offset: reproducible per-point, not wall-clock dependent.
offset = rng.uniform(-jitter, jitter, size=seeds.shape)
stabilized = seeds + offset
# Sort to fix row order -> identical input to the tessellator every run.
order = np.lexsort((stabilized[:,1], stabilized[:,0]))return stabilized[order]
Partition the envelope with an R-tree, process one tile at a time inside a 1.5x buffer halo that captures cross-tile edges, clip every cell to the constraint envelope so no infinite ray survives, and free the per-tile triangulation immediately. Reproject to a metric equal-area CRS first — EPSG:6933 keeps areas true, which the area-conservation gate later depends on; the default EPSG:4326 would make every area comparison meaningless on a sphere.
python
import geopandas as gpd
from shapely import voronoi_polygons, set_precision, MultiPoint
from shapely.geometry import box
EQUAL_AREA ="EPSG:6933"# metric, area-true; motivated -> not the 4326 default
GRID =1e-6# shared precision grid -> snaps coincident edgesdeftessellate_tile(tile_seeds, envelope,*, buffer_frac=1.5):"""Build clipped, snapped cells for one tile against a fixed memory budget."""
halo = box(*tile_seeds.total_bounds).buffer(
buffer_frac * tile_seeds.geometry.distance(tile_seeds.geometry.centroid).max())
local = tile_seeds.clip(halo)# only seeds touching the tile
mp = MultiPoint(local.geometry.tolist())# extend_to forces finite cells -> no unbounded rays to clip by hand.
cells = voronoi_polygons(mp, extend_to=envelope.geometry.iloc[0])
cells = gpd.GeoSeries(list(cells.geoms), crs=EQUAL_AREA)
cells = cells.intersection(envelope.geometry.iloc[0])# clip to zoning extent
cells = set_precision(cells, GRID)# snap to shared griddel mp, local # release tile scratchreturn cells[cells.area >0]# drop emptiesdeftessellate(seeds_gdf, envelope, tiles):
parts =[tessellate_tile(seeds_gdf.clip(t), envelope)for t in tiles]return gpd.GeoSeries(gpd.pd.concat(parts), crs=EQUAL_AREA).reset_index(drop=True)
Because each tile holds only its own seeds plus the halo, peak memory is bounded by the largest tile rather than the whole extent — the per-worker ceiling becomes a tuning knob instead of a cliff. set_precision on a shared grid is what makes the cross-tile seams merge cleanly: two cells that should share an edge snap to identical coordinates instead of leaving a sliver. For external robustness guarantees on the underlying predicates, the clip and overlay steps rely on the exact-arithmetic routines documented in the GEOS C++ Library, and the circumcircle handling aligns with the CGAL 2D Triangulation documentation.
A partition is only correct if it conserves area against the input extent and reproduces byte-for-byte across runs. Wire both as CI assertions so a regression cannot be promoted into the simulation registry — the same gate discipline the scoping rules and data contracts require for every generated artifact.
python
import numpy as np
from shapely import to_wkb
from shapely.ops import unary_union
deftest_partition_is_complete_and_reproducible(seeds_gdf, envelope, tiles):
cells = tessellate(seeds_gdf, envelope, tiles)# 1. Area conservation: union of cells == input extent, within tolerance.
covered = unary_union(cells.values).area
extent = envelope.geometry.iloc[0].area
assertabs(covered - extent)/ extent <1e-4,"gap or overlap in partition"# 2. No overlap: summed cell area must not exceed the union area.assert cells.area.sum()<= covered *(1+1e-9),"cells overlap"# 3. Validity: every geometry closed and non-self-intersecting.assert cells.is_valid.all(),"invalid ring produced"# 4. Reproducibility: a second run is byte-identical.
again = tessellate(seeds_gdf, envelope, tiles)assert[to_wkb(g)for g in cells]==[to_wkb(g)for g in again],"non-deterministic"
The 1e-4 relative tolerance on area conservation is the single most important gate: it catches both gaps (covered extent too small) and overlaps (summed area too large) that geometric is_valid checks alone will pass straight through.
Antimeridian-spanning extents. A zoning envelope that crosses ±180° longitude produces cells whose coordinates wrap, and the equal-area reprojection to EPSG:6933 will fold them onto themselves, collapsing area to near zero and tripping the conservation gate falsely. Split the envelope at the antimeridian before tessellating and reassemble after export, or shift the whole extent into a local projected CRS whose central meridian sits inside the region.
Floating-point precision at the precision-grid boundary. Setting GRID too coarse relative to the smallest legitimate cell will snap distinct vertices together and create the slivers you were trying to dissolve; setting it too fine leaves the original cross-tile gaps unsnapped. Calibrate GRID to roughly one order of magnitude below your minimum zoning resolution, and assert that no cell’s area collapses below the floor after snapping.
Duplicate or near-coincident seeds (null-island artifacts). Seeds that collapse to (0, 0) from an upstream join failure, or two seeds within the precision grid of each other, produce zero-area Voronoi cells that survive as empty geometries. Deduplicate within the grid tolerance before tessellating and assert the post-clip cell count equals the unique seed count, so a silent upstream NaN-to-zero coercion cannot quietly drop a zone.
Why not just raise the memory limit and keep the single Voronoi call?
Peak memory for monolithic construction grows super-linearly with seed count because the triangulation, the diagram, and the clipped cells all coexist in RAM. Raising the ceiling postpones the failure by one extent size; tiled construction removes it by bounding peak memory to the largest single tile, which you control directly with the tile grid.
How big should each tile and its buffer halo be?
Size tiles so the seeds in one tile plus its 1.5x halo fit comfortably under the per-worker memory budget — typically tens of thousands of seeds per tile. The 1.5x buffer must exceed the maximum expected Voronoi radius so every cross-tile edge is captured; too small a halo leaves gaps at tile seams that the area-conservation gate will catch but not locate.
Does the hash-based jitter bias the synthetic zoning statistics?
No, provided the jitter magnitude is small relative to the inter-seed distance — a few orders below the cell radius. It exists only to break exact collinearity so no three seeds share a circumcircle of infinite radius. Because it is drawn from a seeded default_rng, it is fully reproducible and contributes no run-to-run variance.
Why reproject to EPSG:6933 instead of staying in EPSG:4326?
Area conservation is the central correctness gate, and area is undefined for degrees on a sphere. EPSG:6933 is a metric equal-area projection, so cell areas in square metres are comparable and the union-versus-extent tolerance is meaningful. Reproject in, tessellate and gate in the equal-area CRS, then reproject back only for export if a geographic CRS is required downstream.
