Detecting and Removing Sliver Polygons in GeoPandas #

Scenario / symptom: you overlay a solar parcel layer against a county boundary or union two administrative datasets, and the result carries hundreds of extra features — hair-thin, near-zero-area polygons hugging the shared edges. The feature count balloons, dissolve() produces ragged borders, and the total developable-area tally drifts by a fraction of a percent that a permit reviewer will later flag. These artefacts are sliver polygons, and they appear at the geometry-processing stage of the spatial data quality and validation workflow — the moment two layers whose shared borders were digitized independently get intersected or unioned. The parent workflow scores topological consistency as one of its four quality dimensions; slivers are the single most common way that dimension fails silently, because each artefact is a perfectly valid geometry that merely should not exist.

Slivers are not invalid in the Shapely sense — they pass is_valid, they have positive area, and they survive every predicate. That is exactly why a naive explode()-then-count pipeline never raises an exception: it inflates the record count, corrupts area-weighted metrics like MW density, and leaves the layer looking plausible. Removing them requires a geometric definition of “too thin to be real,” a threshold tuned to the source data, and a removal step that either drops the artefact or dissolves it back into the neighbour it was carved from — all executed in a metric CRS so that an area threshold means square metres, not square degrees.

Root-cause analysis #

Three compounding causes produce sliver artefacts, and each maps to a distinct stage of the fix below:

1. Vertex misregistration between layers. The shared boundary between a parcel layer and a county layer is digitized twice, from different sources, at different precisions. The two “identical” edges differ by a metre or two, so an overlay(how="intersection") or union slices the gap between them into long, thin polygons. This is the dominant source and it scales with the length of every shared border in the study area.
2. Coordinate precision and rounding. Reprojection, snapping, and float truncation nudge vertices by sub-metre amounts. A polygon that should share an edge exactly instead overlaps its neighbour by a razor-thin margin, and the overlay emits that margin as its own feature.
3. Absent geometric filtering after the overlay. GeoPandas returns every geometric intersection the algebra produces; it has no concept of “this piece is too small to matter.” Without an explicit thinness-and-area gate, the slivers flow straight into the feature count, the dissolve, and the downstream area sum.

The reliable test for a sliver is not area alone — a small-but-legitimate parcel would be dropped by a pure area filter. A sliver is defined by being both small and thin: it has a large perimeter relative to the area it encloses. The Polsby-Popper compactness score captures this precisely:

where $A$ is the polygon area and $P$ its perimeter. A perfect circle scores $\mathrm{PP}=1$; a long thin sliver scores near $0$. Combining a low compactness score with an absolute area ceiling isolates true artefacts while sparing small legitimate features, because a compact 300 m² parcel scores high on $\mathrm{PP}$ even though its area is tiny.

Pre-flight validation #

Run the detector before any removal so the decision is auditable and the thresholds can be tuned against real numbers. Work in an equal-area metric CRS — EPSG:5070 (NAD83 / Conus Albers) for contiguous-US energy work, or the local UTM zone such as EPSG:32610 for a single-region study — so that area_m2 is genuinely square metres and a metre-scale sliver is comparable across the layer. Computing area or perimeter in geographic EPSG:4326 returns degree-based numbers that make the thinness test meaningless.

The detector below annotates every feature with its area, its Polsby-Popper score, and a boolean is_sliver, without deleting anything. That separation is deliberate: it lets a reviewer inspect the flagged set before it is destroyed.

import numpy as np
import geopandas as gpd

TARGET_CRS = "EPSG:5070"  # NAD83 / Conus Albers Equal Area — area-true metric frame


def flag_slivers(
    parcels_gdf: gpd.GeoDataFrame,
    area_threshold_m2: float = 50.0,
    pp_threshold: float = 0.03,
) -> gpd.GeoDataFrame:
    """Annotate area_m2, Polsby-Popper compactness, and an is_sliver flag.

    A feature is a sliver only if it is BOTH tiny (area below the ceiling)
    AND thin (compactness below the limit) — so a small compact parcel survives.
    """
    if parcels_gdf.crs is None or parcels_gdf.crs.to_epsg() != 5070:
        parcels_gdf = parcels_gdf.to_crs(TARGET_CRS)

    geom = parcels_gdf.geometry
    area_m2 = geom.area
    perimeter_m = geom.length
    # Polsby-Popper: 4*pi*A / P^2 -> ~1 for a disc, ~0 for a hair-thin sliver
    pp = np.where(perimeter_m > 0, (4 * np.pi * area_m2) / np.square(perimeter_m), 0.0)

    out = parcels_gdf.copy()
    out["area_m2"] = area_m2
    out["pp_thinness"] = pp
    out["is_sliver"] = (area_m2 < area_threshold_m2) & (pp < pp_threshold)
    return out

Before committing to a threshold, inspect the distribution — flagged["area_m2"].describe() and a histogram of pp_thinness reveal whether the artefacts cluster far below your candidate cut-off or blur into legitimate small parcels. Genuine slivers from an overlay almost always sit in a tight spike near zero on both axes, well separated from real features, which is what makes the two-part gate robust.

Fix implementation #

The corrected removal function takes the flagged frame and either drops the slivers outright or dissolves each one back into the neighbouring polygon it shares the longest boundary with. Dropping is correct when the sliver is pure overlay noise that encloses no real land; dissolving is correct when the sliver was carved from a legitimate parcel and its area must be conserved. The dissolve_to_neighbour=True path preserves the total area exactly, which matters when the layer feeds a developable-area or capacity tally.

import geopandas as gpd
import pandas as pd


def remove_slivers(
    flagged_gdf: gpd.GeoDataFrame,
    dissolve_to_neighbour: bool = True,
) -> gpd.GeoDataFrame:
    """Remove flagged slivers by dropping them or merging into the best neighbour.

    Assumes flag_slivers() has populated is_sliver. Returns a frame with no
    features below the area/thinness gate and, when dissolving, conserved area.
    """
    slivers = flagged_gdf[flagged_gdf["is_sliver"]]
    keep = flagged_gdf[~flagged_gdf["is_sliver"]].copy()

    if slivers.empty:
        return keep.drop(columns=["is_sliver"], errors="ignore")

    if not dissolve_to_neighbour:
        # Pure noise: drop and log how much area was discarded
        return keep.drop(columns=["is_sliver"], errors="ignore")

    # Attach each sliver to the retained neighbour with the longest shared border
    retained = keep.reset_index(drop=True)
    sindex = retained.sindex
    for _, sliver in slivers.iterrows():
        candidates = retained.iloc[list(sindex.query(sliver.geometry, predicate="intersects"))]
        if candidates.empty:
            continue  # orphan sliver — routed to fallback below
        shared_len = candidates.geometry.intersection(sliver.geometry).length
        target = shared_len.idxmax()
        # Union the sliver geometry into its chosen neighbour, conserving area
        retained.at[target, "geometry"] = (
            retained.at[target, "geometry"].union(sliver.geometry)
        )
    return retained.drop(columns=["is_sliver"], errors="ignore")

The neighbour search uses the layer’s spatial index (sindex.query) to restrict candidates to bounding-box intersections before the exact border-length comparison runs, the same indexing discipline that keeps proximity and distance calculations tractable at scale. Choosing the neighbour by longest shared edge — rather than largest area or first match — is what makes the dissolve deterministic and defensible: the sliver rejoins the polygon it was most plausibly split from.

Fallback routing & performance tuning #

For national parcel layers or CI/CD runs, layer these strategies onto the core functions:

• Snap before you overlay. The cheapest sliver is the one never created. Apply shapely.set_precision(geom, grid_size=0.1) or gdf.geometry.snap() to align shared edges to a common grid before the overlay/union, so misregistered vertices collapse instead of slicing. This complements the geometry repair covered in cleaning messy shapefiles in geopandas.
• Tune the two thresholds independently. Loosen pp_threshold (toward 0.05) to catch fatter overlay chips; tighten area_threshold_m2 (toward 10) when legitimate micro-parcels exist. Never widen area alone — that is what deletes real features.
• Route orphan slivers, do not drop them silently. A sliver with no retained neighbour (an isolated overlay chip) should go to a quarantine layer with an error code, mirroring the remediation-queue pattern of the parent quality gate, not vanish from the record.
• Index once, reuse. Build retained.sindex a single time outside the loop (as above); rebuilding it per sliver turns an O(N) pass into O(N²).
• Chunk by geometry, not by row. For continental layers, partition by a coarse grid or state boundary so each remove_slivers call operates on a self-contained neighbourhood, keeping the union graph and spatial index in memory.

Downstream validation #

Before the cleaned layer feeds a siting model or a regulatory boundary overlay, gate it with a CI/CD assertion that proves the removal did its job without corrupting the total area. This is the check that catches a threshold regression or a dissolve that accidentally dropped real land.

import numpy as np
import geopandas as gpd


def assert_no_slivers(
    cleaned_gdf: gpd.GeoDataFrame,
    original_area_m2: float,
    area_threshold_m2: float = 50.0,
    pp_threshold: float = 0.03,
    area_tol: float = 1e-4,
) -> None:
    """CI/CD gate: no residual slivers, CRS intact, total area preserved."""
    assert cleaned_gdf.crs is not None and cleaned_gdf.crs.to_epsg() == 5070, \
        "cleaned layer must stay in the equal-area metric frame EPSG:5070"

    geom = cleaned_gdf.geometry
    area_m2 = geom.area
    pp = np.where(geom.length > 0, (4 * np.pi * area_m2) / np.square(geom.length), 0.0)
    residual = (area_m2 < area_threshold_m2) & (pp < pp_threshold)
    assert not residual.any(), f"{int(residual.sum())} sliver(s) survived removal"

    # Dissolve-to-neighbour must conserve total area within tolerance
    drift = abs(area_m2.sum() - original_area_m2) / original_area_m2
    assert drift <= area_tol, f"area drifted {drift:.4%} (> {area_tol:.4%}) during removal"

Passing original_area_m2 from the pre-removal frame turns the area check into a hard invariant: a dissolve conserves area exactly, so any drift beyond floating-point tolerance means a sliver was dropped when it should have been merged. Pair this with the geometry-validity and topology assertions from validating geometry topology with Shapely 2 predicates so a single CI step proves the layer is valid, sliver-free, and area-conserving before it ever reaches a feasibility study. Logging the flagged count and the discarded-versus-dissolved split keeps the operation auditable for the same permitting review the parent quality gate feeds.

Related #

• Spatial Data Quality & Validation — the parent workflow whose topological-consistency dimension slivers violate.
• Best practices for cleaning messy shapefiles in GeoPandas — geometry repair and snapping that prevents slivers upstream of the overlay.
• Validating geometry topology with Shapely 2 predicates — the validity and overlap checks that pair with sliver removal in one CI gate.
• Regulatory Boundary Mapping — the setback and jurisdictional overlays whose clean edges depend on sliver-free parcels.