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.

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

$$\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.

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.

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

pythonCopy
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.

pythonCopy
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")

pythonCopy
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

pythonCopy
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.

pythonCopy
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

Related

• Polygon Tessellation Algorithms for Synthetic Spatial Data — the parent workflow for generating tessellated layers this clean pass protects.
• Optimizing Voronoi Tessellation for Synthetic Zoning Maps — a generator whose collinear-seed artifacts produce the same slivers this page removes.
• Enforcing CRS Contracts in a Generation Pipeline — where the precision grid should be pinned alongside the projection.
• Async Execution for Large Grids — the distributed setting where a mismatched grid size re-opens slivers at tile seams.