Coastal & Marine Spatial Analysis Pipelines

Repairing Invalid Coastline Polygons in Shapely

A national shoreline layer opens fine, plots fine, and then every overlay against it raises TopologyException: Input geom 1 is invalid: Self-intersection at .... Nothing about the layer looked wrong, and the coordinate in the message is one of several million. This page is the diagnosis-and-repair procedure for that situation, within the coastal vector processing workflow, and its central claim is that the reflexive fix — buffer(0) — is the wrong one for two of the three ways coastal geometry goes invalid, and that it fails silently in exactly the case that matters most.

Root cause: three different problems wearing one error message

GEOS reports invalidity through a single exception type, so three unrelated data problems arrive looking identical. Separating them is the whole of the diagnosis.

import geopandas as gpd
from shapely.validation import explain_validity

gdf = gpd.read_file("shoreline_mhw.gpkg").to_crs("EPSG:32618")
bad = gdf[~gdf.geometry.is_valid]
for idx, geom in bad.geometry.items():
    print(idx, explain_validity(geom))

# 1187 Self-intersection[454102.31 4293201.88]
# 1188 Self-intersection[454103.0000000001 4293202.0]
# 2044 Hole lies outside shell[451880.5 4290117.2]

Feature 1187 is a genuine self-intersection: a spit digitised so that the exterior ring crosses itself, describing two lobes of land joined at a point. Feature 1188 is a float-noise self-intersection: two vertices that should be identical differ in the last representable bit after a reprojection, and the ring grazes itself. Feature 2044 is a ring-nesting error: an interior ring — a lagoon — that is not contained by the exterior it was recorded against.

They need different treatment. The first is a data decision. The second should never have reached a human. The third is a producer error that a repair can only guess at.

Three Invalidity Classes, Three Different Repairs Panel one shows a bowtie polygon whose exterior ring crosses itself at a spit, labelled a genuine self-intersection requiring a split into two parts. Panel two shows two vertices separated by a fraction of a nanometre where a ring grazes itself, labelled float noise from reprojection and removed by a precision model rather than by a repair. Panel three shows an interior ring drawn partly outside the exterior shell it is recorded against, labelled a ring-nesting error that a geometric repair can only guess at and that should be referred back to the producer. 1 · genuine self-intersection two real lobes, joined at a point → make_valid: split, keep both 2 · float-noise near-touch Δ ≈ 1e-9 m two vertices that should be one → precision model, before repair 3 · hole outside shell a lagoon recorded against the wrong exterior → refer back; do not guess One exception type, three causes. A repair applied without the classification fixes the wrong thing correctly.

Step-by-step fix with production code

1. Eliminate the float-noise class first

A precision model snaps every coordinate to a grid. Choose the grid finer than any real measurement and coarser than float noise — a millimetre is comfortably both for coastal survey data. This removes the entire second class before any repair runs, which typically drops the invalid-feature count by an order of magnitude and leaves a set small enough to look at individually.

import logging

import shapely
from shapely import set_precision

logger = logging.getLogger(__name__)
GRID_M = 0.001  # 1 mm


def apply_precision(gdf: gpd.GeoDataFrame, grid: float = GRID_M) -> gpd.GeoDataFrame:
    """Snap to a fixed grid; returns a new layer, never mutates the input."""
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("apply a precision model in a projected CRS, not degrees")
    before = int((~gdf.geometry.is_valid).sum())
    out = gdf.copy()
    out.geometry = gpd.GeoSeries(
        shapely.set_precision(gdf.geometry.values, grid), crs=gdf.crs, index=gdf.index
    )
    after = int((~out.geometry.is_valid).sum())
    logger.info("precision grid %.4f m: invalid %d → %d", grid, before, after)
    return out

The CRS guard is not decorative. A grid size of 0.001 in a geographic CRS is a grid of one-thousandth of a degree — roughly 111 metres — which will destroy the layer while reporting success.

2. Repair what remains with make_valid

make_valid implements the OGC-defined repair: a self-intersecting ring becomes a multipolygon of the parts it actually bounds. buffer(0) implements something else entirely — it computes a zero-distance buffer, whose result for a bowtie is whichever lobe the algorithm considers dominant. Both return valid geometry; only one returns the same land.

def repair(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """make_valid across the layer, dropping anything that repairs to empty."""
    out = gdf.copy()
    out.geometry = gpd.GeoSeries(
        shapely.make_valid(gdf.geometry.values), crs=gdf.crs, index=gdf.index
    )
    empties = int(out.geometry.is_empty.sum())
    if empties:
        logger.warning("%d feature(s) repaired to empty geometry — dropping", empties)
        out = out[~out.geometry.is_empty].copy()
    return out

make_valid can also return a GeometryCollection when the input mixed dimensions — a polygon with a degenerate zero-area spur repairs to a polygon plus a linestring. That collection will fail a write to most formats. Filter to the dimension you asked for before handing the layer on:

from shapely.geometry import GeometryCollection, MultiPolygon, Polygon


def keep_polygons(geom):
    """Reduce a repair result to its polygonal parts only."""
    if isinstance(geom, (Polygon, MultiPolygon)):
        return geom
    if isinstance(geom, GeometryCollection):
        parts = [g for g in geom.geoms if isinstance(g, (Polygon, MultiPolygon))]
        if not parts:
            raise ValueError("repair produced no polygonal part")
        return shapely.union_all(parts)
    raise TypeError(f"unexpected geometry type after repair: {type(geom).__name__}")

3. Audit the area, and raise

This is the step that separates a repair from a deletion. Area before and after must agree to within a tolerance; anything larger means the repair made a decision that a human should have made.

AREA_TOLERANCE = 1e-4  # 0.01%


def audit_area(before: gpd.GeoDataFrame, after: gpd.GeoDataFrame) -> float:
    """Relative area change across a repair. Raises if it exceeds tolerance."""
    a0 = float(sum(shapely.make_valid(g).area for g in before.geometry))
    a1 = float(after.geometry.area.sum())
    delta = (a0 - a1) / a0 if a0 else 0.0
    logger.info("area %.1f m² → %.1f m² (%.6f%%)", a0, a1, delta * 100.0)
    if abs(delta) > AREA_TOLERANCE:
        raise ValueError(
            f"repair changed layer area by {delta:.4%}; inspect the features "
            "that changed rather than accepting the result"
        )
    return delta

Note that a0 is computed on repaired copies. The area of an invalid polygon is undefined in GEOS — comparing it against anything is meaningless — so both sides of the comparison have to be measured the same way for the number to mean anything.

The Precision Model Does Most of the Work A funnel of invalid feature counts through four stages of the repair sequence for a national shoreline layer of 1.4 million polygons. On load, 41,800 features are invalid. After the millimetre precision model, 3,100 remain. After make_valid, 12 remain. Those 12 are reviewed individually and turn out to be ring-nesting errors that are referred back to the data producer. A caption notes that the first stage removes 93 per cent of the problem without any geometric decision being taken, which is what makes the residue small enough to look at. 1.4 M polygon national shoreline — invalid features per stage on load 41,800 precision 1 mm 3,100 — float noise gone, no decision taken make_valid 12 — small enough to open individually reviewed 12 ring-nesting errors → referred back to the producer, not repaired 93% of the problem was never a geometry problem, and nothing about the coastline was decided to remove it.

Verification and acceptance test

The acceptance test asserts three things at once: the layer is valid, the area survived, and — the check people forget — the feature count did not change in a way nobody asked for.

def test_shoreline_repair_preserves_land():
    raw = gpd.read_file("shoreline_mhw.gpkg").to_crs("EPSG:32618")
    snapped = apply_precision(raw)
    fixed = repair(snapped)
    fixed.geometry = fixed.geometry.map(keep_polygons)

    assert bool(fixed.geometry.is_valid.all()), "layer still invalid after repair"
    assert abs(audit_area(raw, fixed)) <= AREA_TOLERANCE
    # A bowtie legitimately becomes one MultiPolygon, not two rows, so the row
    # count must be unchanged — a drop means features repaired to empty.
    assert len(fixed) == len(raw), f"{len(raw) - len(fixed)} feature(s) vanished"

Run it against the real layer in CI rather than against a synthetic fixture. Synthetic bowties are easy and the interesting failures in coastal data are not the ones anyone thinks to construct.

Why buffer(0) became the folk remedy, and what it costs

buffer(0) predates make_valid by roughly a decade. Before GEOS 3.8 it was the only repair available from Shapely, it is one call, and it works — in the narrow sense that its output is always valid. That combination made it the answer on every mailing list, and the advice outlived the reason for it.

What it actually does is compute the zero-distance offset curve of the input. For a polygon with a hole nested correctly, that is the identity and the call is harmless. For a bowtie, the offset construction produces the region the algorithm considers to be on the interior side of the ring, and for a self-crossing ring that region is one lobe. Which lobe depends on the winding, not on which one is the real headland.

The failure has three properties that make it unusually dangerous in coastal work. It is silent: no exception, no warning, and the output passes every validity assertion. It is plausible: a coastline missing one headland still looks like a coastline. And it is not idempotent-detectable: running the repair again changes nothing, so a re-run offers no signal that something was lost the first time.

The concrete test is the area audit above, and it takes one line. On a national shoreline the difference is typically small in aggregate — tenths of a per cent — which is exactly why it survives review. The features it removes are small, and small coastal features are disproportionately the ones that matter: spits, breakwaters, dredged channel margins, the sand bar that closed last winter.

There is one legitimate remaining use. buffer(0) is still a fast way to dissolve a set of overlapping polygons into their union footprint when you genuinely want the footprint and do not care about the parts. Written as union_all it says so; written as buffer(0) it reads as a repair and the next person to touch the code will treat it as one.

Which Repair, and When A four-row decision table. Row one: symptom is invalid features numbering in the thousands with coordinates differing in the last bits, cause is reprojection float noise, action is a precision model, and the risk if the wrong tool is used is that make_valid churns through the whole layer for nothing. Row two: symptom is a self-intersection at a spit or jetty, cause is a genuine digitising crossing, action is make_valid, and the risk with buffer zero is that a lobe of land is deleted silently. Row three: symptom is hole lies outside shell, cause is a producer error in ring nesting, action is to refer it back rather than repair, and the risk of repairing is that the guess is recorded as fact. Row four: symptom is many overlapping polygons where only the footprint is wanted, action is union_all, and the note is that buffer zero does the same thing while reading as a repair. symptomcausedo thisrisk if you don’t thousands invalid, last-bit deltasreprojection noiseset_precisionrepair churns for nothing self-intersection at a spitreal digitising crossingmake_validbuffer(0) deletes a lobe hole lies outside shellproducer ring errorrefer backa guess becomes a fact overlaps, footprint wantednot an error at allunion_allreads as a repair later Only the second row is a repair in the ordinary sense. The first is hygiene, the third is a data-quality referral, and the fourth is an intentional operation that has been misfiled as a repair for a decade. Classifying first is what makes each row cheap; skipping the classification is what makes all four expensive.

Edge cases and gotchas

  • make_valid has two methods in GEOS 3.10+. The default "linework" preserves every input vertex and can produce collections; "structure" produces a cleaner polygonal result but discards degenerate parts. Shapely 2.1 exposes the choice as make_valid(geom, method="structure"). For shoreline data the default is usually right, because a degenerate spur is often a real jetty digitised badly and worth seeing rather than dropping.
  • Repairing before reprojecting wastes the work. A reprojection re-introduces float noise at the new coordinates, so a layer repaired in EPSG:4326 and then projected to UTM can arrive invalid again. Reproject first, then apply the precision model in the working CRS, then repair.
  • A valid layer can still be a broken partition. is_valid is a per-feature property. Two features that overlap are each perfectly valid and together describe two square kilometres of land in one square kilometre of space. That needs the self-overlap gate from the parent topic, not a validity check.

Up: Vector Shoreline & Coastal Geometry Processing