Coastal & Marine Spatial Analysis Pipelines

Simplifying Shorelines Without Breaking Topology

A shoreline digitised at 1:10,000 carries far more vertices than a 50-metre analysis grid can use. Dropping them is obviously worthwhile — every spatial predicate, every render, every serialization pays for vertices — and the obvious call, gdf.geometry.simplify(25), produces a layer that is smaller, looks identical at review zoom, and is quietly no longer a partition of the coast. This page explains why, within the coastal vector processing workflow, and gives the shared-edge procedure that gets the same reduction without the damage.

Root cause: two features, two independent copies of one boundary

In a simple-features layer, adjacency is not represented. Two polygons that share a boundary each store their own complete copy of it, vertex for vertex. They agree exactly because the digitiser snapped them, not because the format knows they are related.

Douglas-Peucker is deterministic but not local: which vertices it keeps on a run of coordinates depends on the endpoints of the run and on the vertices around it. Polygon A’s copy of the shared boundary is embedded in A’s ring, with A’s other vertices before and after it; polygon B’s copy sits in a different ring with different neighbours. The two copies are therefore simplified against different context and come out different.

The result is not a small error. Wherever the two simplified lines cross, you get alternating overlaps and gaps — slivers — a few metres wide and running the length of every shared boundary in the layer.

import geopandas as gpd

zones = gpd.read_file("management_zones.gpkg").to_crs("EPSG:32618")
print(zones.geometry.area.sum())          # 4_182_119_884.2  (the true total)

naive = zones.copy()
naive.geometry = zones.geometry.simplify(25.0)
print(naive.geometry.area.sum())          # 4_182_004_711.9  (looks fine)

# The damage is not in the total; it is in the pairwise relationships.
pairs = naive.sindex.query(naive.geometry, predicate="overlaps")
print(len(pairs[0]))                      # 3_884 overlapping pairs, from 0

The area total barely moves, because the slivers gain and lose roughly equally. Every check that looks at the layer as a whole passes. The check that fails is the one nobody runs.

One Boundary, Simplified Twice, Becomes Two Boundaries Left panel: a dense zigzag boundary between two zones is drawn, with zone A's simplified copy as a solid line and zone B's simplified copy as a dashed line. The two diverge, producing four shaded sliver regions where they overlap and three unshaded gaps where neither covers the space. Right panel: the same boundary extracted once as an arc, simplified once and used by both zones, so a single line separates them with no slivers or gaps. Below, counts are given: 3,884 overlapping pairs and 41,200 square metres of gap on the left, zero and zero on the right, for the same 62 per cent vertex reduction. simplified per feature zone A copy — solid · zone B copy — dashed shaded: both zones claim it unshaded between them: neither does 3,884 overlapping pairs · 41,200 m² of gap simplified as a shared arc one arc, referenced by both zones the partition is preserved by construction, not by luck 0 overlapping pairs · 0 m² of gap Both panels achieved the same 62% vertex reduction. Only one of them still describes a partition of the coast.

Step-by-step fix with production code

1. Choose the tolerance from the analysis, not the file size

The tolerance is a statement about what spatial detail the downstream analysis can act on. A habitat model on a 100 m grid cannot use detail below roughly half a cell; a dredging boundary that determines a pay volume cannot lose a metre. Deriving it:

def tolerance_for_grid(cell_size_m: float, safety: float = 0.4) -> float:
    """Simplification tolerance for a layer feeding a raster analysis.

    Keeps displacement well inside a cell so a boundary never crosses into a
    neighbouring cell as a result of simplification alone.
    """
    if cell_size_m <= 0:
        raise ValueError("cell_size_m must be positive")
    return cell_size_m * safety


tolerance_for_grid(100.0)   # 40.0 m  — habitat model
tolerance_for_grid(2.0)     # 0.8 m   — dredge design surface

Where one layer feeds analyses at both ends of that range, publish it twice with the tolerance in the filename. A compromise tolerance is wrong for both consumers and, unlike two files, cannot be corrected downstream.

2. Build a topology, simplify the arcs, rebuild

The Python route is the topojson package, which constructs the shared-arc representation, applies the tolerance per arc, and returns a GeoDataFrame. The arcs are the unit of simplification, so a boundary shared by two features is simplified exactly once.

import logging

import geopandas as gpd
import topojson as tp

logger = logging.getLogger(__name__)


def simplify_preserving_topology(
    gdf: gpd.GeoDataFrame, tolerance_m: float
) -> gpd.GeoDataFrame:
    """Shared-edge simplification. Input must be in a projected CRS."""
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("tolerance is in CRS units; use a projected CRS")
    if not bool(gdf.geometry.is_valid.all()):
        raise ValueError("simplify a valid layer only; repair it first")

    before = int(gdf.geometry.count_coordinates().sum())
    topo = tp.Topology(gdf, prequantize=False, shared_coords=True)
    out = topo.toposimplify(tolerance_m).to_gdf()
    out.crs = gdf.crs
    after = int(out.geometry.count_coordinates().sum())

    logger.info(
        "toposimplify %.2f m: %d → %d vertices (%.1f%% reduction)",
        tolerance_m, before, after, 100.0 * (1.0 - after / before),
    )
    return out

prequantize=False matters. Quantization snaps coordinates to an integer lattice before simplification and is a second, hidden simplification with its own tolerance; leaving it on means the effective displacement is larger than the tolerance you asked for and you cannot say by how much.

Outside Python, mapshaper -simplify and PostGIS ST_SimplifyPreserveTopology do the same job; PostGIS’s version preserves topology within a single geometry but not across features, so for a multi-feature partition it is ST_CoverageSimplify (PostGIS 3.4+) that is wanted, not the older function whose name suggests otherwise.

3. Verify the partition survived 

def assert_still_a_partition(gdf: gpd.GeoDataFrame, max_sliver_m2: float = 1.0) -> None:
    """No feature pair may overlap, beyond a numerical-noise tolerance."""
    left, right = gdf.sindex.query(gdf.geometry, predicate="overlaps")
    total = 0.0
    for i, j in zip(left, right):
        if i >= j:
            continue  # each pair once
        total += gdf.geometry.iloc[i].intersection(gdf.geometry.iloc[j]).area
    if total > max_sliver_m2:
        raise ValueError(f"simplification opened {total:,.1f} m² of overlap slivers")
    logger.info("partition intact: %.3f m² total overlap", total)
Most of the Reduction Is Bought by a Very Small Tolerance A curve of remaining vertex count against simplification tolerance in metres for a national shoreline. At zero tolerance the layer holds 41 million vertices. By 1 metre it holds 14 million, by 5 metres 6 million, by 25 metres 3.9 million, and by 100 metres 3.1 million, so the curve is steep at the start and nearly flat past about 25 metres. Two vertical markers show the tolerance limits set by two downstream analyses: 0.8 metres for a dredge design surface and 40 metres for a habitat model. A caption observes that the dredge limit still buys a two-thirds reduction, so the cautious tolerance is not an expensive one. simplification tolerance (m) vertices retained 0.8 m — dredge surface still a 66% reduction 40 m — habitat model 90% reduction, no further gain past here 040100250 41 M vertices at zero tolerance; 14 M at one metre. Choosing the tolerance the strictest consumer allows costs surprisingly little, because the curve does most of its work early.

Verification and acceptance test 

def test_simplification_is_lossless_for_the_partition():
    zones = gpd.read_file("management_zones.gpkg").to_crs("EPSG:32618")
    simple = simplify_preserving_topology(zones, tolerance_for_grid(100.0))

    assert_still_a_partition(simple)
    assert len(simple) == len(zones)
    rel = abs(simple.geometry.area.sum() - zones.geometry.area.sum()) / zones.geometry.area.sum()
    assert rel < 1e-3, f"area changed by {rel:.4%}"
    # Displacement bound: no simplified boundary may move further than the tolerance.
    assert simple.geometry.hausdorff_distance(zones.geometry).max() <= 40.0

The Hausdorff assertion is the one worth keeping. Area and count can both survive a simplification that moved a boundary a kilometre in one place and compensated elsewhere; the Hausdorff distance is the direct statement of “nothing moved further than I authorised”.

What Simplification Costs Downstream, and How to Say So

A simplified shoreline is a different shoreline, and the difference propagates into every product built on it. The three effects worth quantifying before publishing are length, area and the position of the boundary itself.

Length falls, always, and by a lot. Coastlines are close to fractal at survey scales, so their measured length is a function of the ruler. Simplifying at 40 m typically removes 15–30% of the measured length of a rocky coast and much less of a straight sandy one. Any product that reports “kilometres of shoreline” — habitat extent, erosion exposure, permit frontage — is reporting a number that depends on a parameter buried three stages upstream. Publish the tolerance next to the figure or the figure is not reproducible.

Area is nearly conserved and locally is not. Douglas-Peucker cuts corners in both directions, so a long boundary loses and gains area in roughly equal measure and the total barely moves. Locally the picture is different: a headland can lose a hectare while the bay beside it gains one. Zonal statistics over small polygons — a single berth, a monitoring quadrat — feel this directly even when the layer-wide area check passes.

Boundary position moves by up to the tolerance, by construction. That is what the tolerance means, and it is why the Hausdorff assertion above is the honest test. A downstream consumer computing distance-to-shore is inheriting that displacement as a systematic error, so it belongs in their uncertainty budget, which means it has to reach them.

The mechanism for all three is the same and it is cheap: write simplification_tolerance_m, vertices_before, vertices_after and max_displacement_m into the layer metadata at the moment the simplification runs. Four numbers, generated rather than typed, and every question above becomes answerable without re-running anything.

The Same Tolerance Removes Very Different Amounts of Coast Two curves of measured shoreline length against simplification tolerance, normalised to the unsimplified length. The rocky-coast curve falls steeply, losing about 28 per cent of its measured length by a 40 metre tolerance. The sandy-coast curve falls gently, losing about 6 per cent at the same tolerance. A vertical marker at 40 metres shows the two curves widely separated there. A caption notes that both segments are in the same layer and were simplified with the same parameter, so any per-segment length statistic is now incomparable between them unless the tolerance is published alongside. simplification tolerance (m) length retained rocky coast — 28% gone at 40 m sandy coast — 6% gone at 40 m 40 m 040100200 One parameter, two very different edits to the same layer. Report frontage from a simplified layer without stating the tolerance and the number is not reproducible, even by you.

Edge cases and gotchas

  • The outer coastline is not shared with anything, and still must not move. Where a layer’s outer boundary is also the boundary of a published national shoreline, simplifying it independently makes the two disagree. Either clip to the authoritative line after simplification, or exclude the outer ring from simplification and accept the vertex cost.
  • Islands smaller than the tolerance disappear. Topology-preserving simplification will not create invalid geometry, but a 30 m islet under a 40 m tolerance collapses and is dropped. Check the feature count, and if small islands matter, set a minimum-area floor and exclude those features from the simplification pass entirely.
  • prequantize defaults to True in topojson. Left on, it snaps to a lattice sized from the data extent before your tolerance is applied, so the effective displacement is larger than requested and varies with the layer’s bounding box. Set it off, or measure the Hausdorff distance and stop trusting the tolerance parameter.

Up: Vector Shoreline & Coastal Geometry Processing