Fixing Sliver Polygons in Tessellation Output

A synthetic tessellation that should partition an area cleanly instead carries hundreds of hair-thin sliver polygons along shared boundaries, each a few square centimetres of area wedged between cells that were supposed to be coincident.

Part of Polygon Tessellation Algorithms for Synthetic Spatial Data: that workflow covers generating tessellated zoning and coverage layers; this page isolates the sliver defect — the near-zero-area artifacts that appear when two boundaries that should be identical differ by a floating-point epsilon — and the two-stage fix that removes them while conserving total area.

Root Cause: Coordinates That Should Be Equal Differ by an Epsilon

A tessellation is supposed to be edge-matched: adjacent cells share an identical boundary, vertex for vertex. Slivers appear when two boundaries that represent the same line carry coordinates that differ in their least significant bits. When you overlay, union, or difference those cells, the geometry engine treats the two almost-coincident lines as distinct, and the microscopic gap or overlap between them becomes a valid polygon with a valid — but tiny — area.

The epsilon has three common origins. Reprojection through pyproj transforms coordinates with transcendental functions whose results are rounded to float64; the same input point routed through two transform paths lands at two slightly different outputs. Independent construction — computing a shared edge once from cell A’s parameters and once from cell B’s — accumulates different rounding. Repeated boolean operations in GEOS introduce new intersection vertices computed in floating point, so each union nudges shared boundaries apart.

Cell-area distribution showing slivers as a separate population, with the threshold in the gap Cell area runs along a logarithmic horizontal axis spanning about seven orders of magnitude. Two populations are plotted. The larger, in the primary colour, is the real cells: a single log-normal hump centred around a few thousand square metres. The smaller, in the accent colour, sits far to the left, centred around a few thousandths of a square metre — the slivers produced where two boundaries that should have been identical differ by a floating-point epsilon. Between the two populations is a wide empty region containing no cells at all, spanning several orders of magnitude. A vertical line marks the geometric midpoint of that gap, which is where an automatic sliver threshold belongs, and the value is printed. The point the figure makes is that this is why the threshold can be automated safely: slivers are not merely small cells on a continuum, they are a distinct population separated by a gap so wide that any threshold inside it gives the same answer. A note adds the corollary — if a tessellation's area histogram has no gap, the small cells are real and dissolving them is data loss rather than repair. Slivers are a separate population, not the small end of one 10⁻⁴ 10⁻³ 10⁻² 10⁻¹ 1 10¹ 10² 10³ 10⁴ 0 50 100 150 200 250 300 cell area (m², log scale) cells threshold anywhere in the gap → 11.6 m² real cells (900) slivers (140) Simulated from a tessellation with a known sliver defect, fixed seed. The gap spans 3.9 orders of magnitude, so every threshold inside it gives the same answer — which is what makes the rule safe to automate. If your own area histogram has no gap, the small cells are real and dissolving them is data loss.
The gap between the two populations spans nearly four orders of magnitude, which is what makes an automatic area threshold safe.

Floating-point coordinates carry a resolution that varies with magnitude. A float64 has 52 bits of mantissa, so near a coordinate of magnitude xx the spacing between representable values is on the order of 2−52 ∣x∣2^{-52}\,|x|. In a projected CRS with coordinates in the millions of metres, that is roughly

Δ≈2−52×106≈2×10−10 m\Delta \approx 2^{-52} \times 10^{6} \approx 2 \times 10^{-10}\ \text{m}

per representable step near the origin, but the accumulated error from a chain of transforms and boolean operations is far larger — routinely nanometres to micrometres — and that is more than enough to open a sliver. Two boundaries agreeing to 10 significant figures still differ, and the engine faithfully reports the difference as area. Enforcing a single explicit precision instead of trusting float64’s full resolution is what collapses the two lines into one, and it composes with the CRS contracts a generation pipeline enforces so precision is pinned alongside the projection.

Two-stage sliver removal: precision-grid snapping then dissolve-by-longest-shared-edge Three panels connected by arrows. Panel one, labelled epsilon-doubled boundary: two large adjacent quadrilateral cells A and B whose shared edge is drawn as two almost-parallel lines a tiny distance apart, with a thin shaded sliver polygon between them annotated area approximately zero. Panel two, labelled snap to precision grid: a faint dotted grid overlays the cells and the two boundary lines have collapsed onto one line of shared grid vertices, closing the sliver. Panel three, labelled dissolve sub-threshold cell: a small residual cell is merged into the neighbour with which it shares its longest edge, shown by an arrow into cell B, with a note that total area is conserved. Epsilon-doubled boundary A B sliver area ≈ 0 snap Snap to precision grid one shared edge, sliver closed dissolve Dissolve sub-threshold cell A B merged into longest-edge neighbour total area conserved
Snapping every vertex to a shared precision grid collapses the doubled boundary into one line; a dissolve step then absorbs any residual sub-threshold cell into its longest-shared-edge neighbour without losing area.

Minimal Reproducer: Manufacture a Sliver

Build two rectangles that share an edge, then perturb one shared vertex by a sub-micrometre epsilon and overlay. Pin the toolchain so the GEOS version and precision behaviour match across machines.

geopandas==0.14.*
shapely==2.*         # GEOS 3.11+; set_precision, overlay
numpy==1.26.*

python
import geopandas as gpd
from shapely.geometry import Polygon

eps = 1e-7  # ~100 nm in a metric CRS — well below any real feature, above float noise

a = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
# B shares the x=10 edge, but its two shared vertices are nudged by eps.
b = Polygon([(10 + eps, 0), (20, 0), (20, 10), (10 - eps, 10)])

gdf = gpd.GeoDataFrame(geometry=[a, b], crs="EPSG:6933")   # equal-area metric CRS
overlay = gpd.overlay(gdf, gdf.iloc[[0]].rename_axis("i"), how="union", keep_geom_type=True)

areas = overlay.geometry.area.sort_values()
print("smallest pieces (m^2):", areas.head(3).tolist())
# a hair-thin sliver of area ~1e-6 m^2 appears between A and B

Any piece whose area is orders of magnitude below the smallest legitimate cell — here anything near 1e-6 m^2 against 100 m² cells — is a sliver, not a feature.

Fix: Snap to a Shared Precision Grid, Then Dissolve Sub-Threshold Cells

The fix is two stages. First, snap every geometry to one explicit precision grid with shapely.set_precision, forcing coordinates that agree within the grid size onto identical values so doubled boundaries collapse. Second, dissolve any cell whose area still falls below a threshold into the neighbour with which it shares the longest edge, so no area is discarded.

Surviving slivers and destroyed detail against snapping grid step, with the safe window The horizontal axis is the precision-grid step in degrees, on a logarithmic scale spanning from a billionth of a degree to about a thousandth. Two curves are plotted against a shared percentage axis, and the safe window between them is just over an order of magnitude wide. The first, falling, is the share of sliver-causing coordinate disagreements that survive snapping: at a very fine step almost all of them survive, and it reaches zero once the step exceeds the largest disagreement. The second, rising, is the share of genuine small features destroyed by snapping: it is zero while the step is smaller than the finest real detail and climbs steeply once it is not. Between the point where the first curve reaches zero and the point where the second leaves zero is a shaded safe window, and its width in orders of magnitude is printed. The note underneath gives the rule: set the step from the geometry rather than by feel — an order of magnitude above the coordinate disagreement you actually measure, and an order of magnitude below the smallest feature the contract says must survive — and if those two constraints do not leave a window, the fix is upstream in how the boundaries were constructed, not in the snapping. Set the step from the geometry, not by feel 1e-9 1e-8 1e-7 1e-6 1e-5 1e-4 0% 25% 50% 75% 100% snapping grid step (degrees, log scale) share of cells (%) safe window — 0.6 orders of magnitude slivers surviving the snap genuine features destroyed by the snap 700 cells, coordinate disagreement ~N(0, 1e-7°), finest genuine feature ~3e-5°, fixed seed. The rule: an order of magnitude above the disagreement you measure, an order below the smallest feature the contract protects. If those two leave no window, the bug is in how the boundaries were built.
Set the grid step an order of magnitude above the coordinate disagreement you measure and an order below the smallest feature the contract protects.
python
import numpy as np
import geopandas as gpd
from shapely import set_precision

def snap_to_grid(gdf: gpd.GeoDataFrame, grid_size: float) -> gpd.GeoDataFrame:
    """Stage 1: snap every vertex to a shared precision grid.
    grid_size is the smallest coordinate distinction kept, in CRS units (metres).
    Pick it well below the smallest real feature but above accumulated float error."""
    snapped = gdf.copy()
    # set_precision rounds coordinates to the grid and re-noding removes the doubled edge.
    snapped["geometry"] = set_precision(gdf.geometry.values, grid_size=grid_size)
    snapped = snapped[~snapped.geometry.is_empty & snapped.geometry.notna()]
    return snapped.set_geometry("geometry")

python
import geopandas as gpd
from shapely.ops import unary_union

def dissolve_slivers(gdf: gpd.GeoDataFrame, min_area: float) -> gpd.GeoDataFrame:
    """Stage 2: merge each sub-threshold cell into its longest-shared-edge neighbour.
    Longest shared edge — not nearest centroid — keeps the merge topologically sensible."""
    gdf = gdf.reset_index(drop=True)
    keep = gdf[gdf.geometry.area >= min_area].copy()
    slivers = gdf[gdf.geometry.area < min_area]
    sindex = keep.sindex
    for _, sliver in slivers.iterrows():
        candidates = keep.iloc[list(sindex.query(sliver.geometry, predicate="touches"))]
        if candidates.empty:
            continue
        # shared-edge length = length of the boundary intersection with each candidate
        shared = candidates.geometry.intersection(sliver.geometry.boundary).length
        target = shared.idxmax()                                  # longest shared edge wins
        keep.at[target, "geometry"] = unary_union(
            [keep.at[target, "geometry"], sliver.geometry]        # absorb the sliver's area
        )
    return keep

python
def clean_tessellation(gdf, grid_size=1e-3, min_area=1e-2):
    """Full pass: snap to a 1 mm grid, then dissolve cells under 0.01 m^2."""
    snapped = snap_to_grid(gdf, grid_size=grid_size)     # closes epsilon-doubled edges
    return dissolve_slivers(snapped, min_area=min_area)  # absorbs any residual splinters

Grid size and area threshold are the two dials. Set grid_size an order of magnitude below your smallest meaningful coordinate distinction — 1 mm in a metric CRS is generous for zoning-scale work — and set min_area to the smallest area any legitimate cell could have. Snapping alone closes the epsilon-doubled boundaries; the dissolve stage only catches genuine splinters left by irregular input, which is exactly the kind of artifact that also appears in Voronoi zoning maps near collinear seed points.

Verification Step: Assert No Slivers and Conserved Area

Gate the cleaned output on two invariants: every cell exceeds the minimum area, and total area is preserved to within the snapping tolerance. The area-conservation check is what proves the dissolve absorbed slivers rather than dropping them.

python
import numpy as np

def test_tessellation_has_no_slivers():
    grid_size, min_area = 1e-3, 1e-2
    raw = build_tessellation()                       # the pipeline under test
    cleaned = clean_tessellation(raw, grid_size, min_area)

    # 1. No cell below the sliver threshold survives.
    assert cleaned.geometry.area.min() >= min_area, "sliver survived the clean pass"

    # 2. Total area conserved: dissolve moved area, never discarded it.
    #    Tolerance scales with grid snapping (grid_size * total boundary length).
    tol = grid_size * raw.geometry.length.sum()
    assert abs(cleaned.geometry.area.sum() - raw.geometry.area.sum()) < tol, "area not conserved"

    # 3. No self-intersections or invalid geometries introduced by snapping.
    assert cleaned.geometry.is_valid.all(), "snapping produced invalid geometry"

    # 4. Cell count dropped by exactly the sliver count (nothing spurious created).
    assert len(cleaned) <= len(raw), "clean pass created new cells"

The tolerance on the area check is not a magic number: snapping can legitimately move each boundary vertex by up to half a grid cell, so total area may shift by about grid_size times the total boundary length. Deriving tol from the geometry rather than hard-coding it keeps the gate correct when the grid size or the layout changes.

Edge Cases & Gotchas

Snapping that collapses a thin-but-real cell. A legitimately narrow feature — a road median, a drainage easement — can be thinner than the grid size and will be snapped out of existence or merged away with the slivers. Set grid_size below the width of the narrowest real feature, and log every dissolved cell’s pre-merge area so a genuine thin feature disappearing is caught in review rather than silently absorbed.

Antimeridian-spanning tessellations. Snapping in EPSG:4326 degrees near ±180° longitude applies a uniform grid to a coordinate that wraps, so a cell straddling the seam has vertices at +179.999° and −179.999° that snapping treats as far apart, tearing the cell. Reproject into a projected CRS whose central meridian sits inside the region before snapping, or split the layer at the antimeridian, clean each side, and re-merge.

Precision mismatch across a distributed pipeline. When tiles are cleaned independently by parallel workers, a grid_size that differs even slightly between workers re-opens slivers at tile seams — the same class of reproducibility hazard that asynchronous large-grid execution must guard against. Pin grid_size as a single pipeline-level constant, derived from the CRS contract, so every worker snaps to the identical grid and seam boundaries stay coincident.

Frequently Asked Questions

Why do slivers appear even when my cells are supposed to share exact boundaries?
Because a shared boundary computed twice — once from each cell, or once through each of two reprojection paths — accumulates different floating-point rounding, so the two lines differ in their least significant bits. The geometry engine treats the near-coincident lines as distinct and reports the microscopic gap between them as a valid polygon. Snapping both to one precision grid forces the coordinates to be identical and closes the gap.
How do I choose the precision grid size?
Set the grid size an order of magnitude below the smallest coordinate distinction that is meaningful for your data, but well above the accumulated floating-point error — in a metric CRS, 1 mm is generous for zoning-scale tessellations. Too coarse and you collapse legitimately thin features; too fine and it fails to close the epsilon-doubled boundaries you are trying to merge.
Why dissolve into the longest-shared-edge neighbour rather than the nearest one?
A sliver is a boundary artifact, so the neighbour it shares the most boundary with is almost always the cell it was split from. Merging by longest shared edge keeps the result topologically sensible and preserves the intended adjacency, whereas nearest-centroid can attach the sliver to a cell it barely touches, distorting the partition.
How do I know the fix conserved area instead of deleting slivers?
Assert that the total area after cleaning equals the total before, within a tolerance of the grid size times the total boundary length. Because the dissolve merges each sliver into a neighbour rather than dropping it, the sum of areas is preserved up to the small shift snapping introduces. If total area falls by roughly the summed sliver area, cells are being deleted, not dissolved.