Vector Shoreline & Coastal Geometry Processing #

Operational intent: this topic covers the vector half of a coastal pipeline — the shoreline polygons, channel centrelines, habitat boundaries and survey footprints that every raster product is eventually clipped, masked or attributed against. It sits inside Marine Spatial Data Fundamentals & Architecture, downstream of CRS alignment and upstream of every masking and zonal-statistics stage in the fleet. The scale is awkward rather than large: a national shoreline is a few hundred megabytes of geometry, small enough to hold in memory and complex enough that a single invalid ring will abort a spatial join four stages later, after twenty minutes of raster work has already been done.

Coastal vector data has a property that inland vector data does not: its boundary is a modelled surface, not an observed line. A shoreline is the intersection of the land with a tidal datum, so “the coastline” is only meaningful once you say which datum — mean high water, mean lower low water, the highest astronomical tide — and the same physical beach yields lines tens or hundreds of metres apart under different choices. Everything below assumes that choice has been made explicitly and recorded; if it has not, no amount of geometric rigour downstream will make the result mean anything.

Reference Configuration #

The Three Ways Coastal Geometry Is Invalid #

GEOS defines validity precisely, and coastal data violates it in three characteristic ways that have different causes and different correct repairs. Conflating them is why buffer(0) has such a reputation: it “fixes” all three by producing output that is valid, and for two of them the output is wrong.

Self-intersecting rings. The exterior ring crosses itself, typically at a narrow spit, a jetty digitised as a zero-width sliver, or a river mouth where two survey lines were joined. GEOS reports Self-intersection with a coordinate. The correct repair splits the ring into the polygons it actually describes, which make_valid does and buffer(0) does not — buffer(0) discards whichever lobe it considers the smaller, silently deleting a real headland.

Ring order and nesting errors. An interior ring — a lagoon, a harbour basin — that is not contained by its exterior, or that is wound in the same direction as the exterior. Shapefiles encode nesting by winding order, and a producer that got the winding wrong emits a file that opens without complaint and describes land where there is water. The repair is to recompute containment from geometry rather than to trust the recorded order.

Duplicate and near-duplicate vertices. Two vertices at the same location, or a few nanometres apart after a reprojection. These are legal in isolation but make almost every predicate fragile, and they are the direct cause of the “worked yesterday, fails today” class of GEOS errors, because the failure depends on floating-point rounding in the specific operation. The repair is a precision model applied once, at ingestion, before any spatial work.

The distinction between the first and third matters operationally. A self-intersection is a data error that a human should look at. A near-duplicate vertex is float noise that should never have reached a human. Applying a precision grid at ingestion removes the entire third category, which in practice is most of the volume, and leaves a much smaller set of genuine geometry problems that are worth reviewing individually.

Production Repair: Precision, Then Validity, Then Audit #

The order matters. Snapping to a precision grid first removes the noise-driven invalidities so that the validity pass deals only with real geometry problems; auditing area afterwards catches the case where a repair succeeded and lost land. Each step logs rather than prints, and the function raises rather than returning a silently degraded layer.

"""Coastal vector repair: precision model, validity, and an area audit."""
from __future__ import annotations

import logging
from dataclasses import dataclass

import geopandas as gpd
import shapely
from shapely import set_precision
from shapely.geometry.base import BaseGeometry

logger = logging.getLogger(__name__)

PRECISION_GRID_M = 0.001      # 1 mm; below any real survey resolution
AREA_LOSS_TOLERANCE = 1e-4    # 0.01% of layer area


@dataclass(frozen=True)
class RepairReport:
    """What the repair actually did, so the caller can gate on it."""
    n_features: int
    n_invalid_in: int
    n_repaired: int
    area_before_m2: float
    area_after_m2: float

    @property
    def relative_area_loss(self) -> float:
        if self.area_before_m2 <= 0.0:
            return 0.0
        return (self.area_before_m2 - self.area_after_m2) / self.area_before_m2


def _valid_area(geom: BaseGeometry | None) -> float:
    """Area of a geometry that may be invalid — measured on a repaired copy."""
    if geom is None or geom.is_empty:
        return 0.0
    return shapely.make_valid(geom).area


def repair_coastal_layer(
    gdf: gpd.GeoDataFrame,
    *,
    grid_size: float = PRECISION_GRID_M,
    max_area_loss: float = AREA_LOSS_TOLERANCE,
) -> tuple[gpd.GeoDataFrame, RepairReport]:
    """Snap to a precision grid, repair validity, and refuse to lose land.

    Raises
    ------
    ValueError
        If the layer has no projected CRS (areas would be meaningless), or if
        the repair removed more area than `max_area_loss`.
    """
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError(
            "repair_coastal_layer requires a projected CRS; got "
            f"{gdf.crs}. Reproject before repairing so areas are in metres."
        )

    geoms = gdf.geometry
    n_invalid_in = int((~geoms.is_valid).sum())
    area_before = float(sum(_valid_area(g) for g in geoms))

    # 1. Precision model first: removes float-noise self-intersections so the
    #    validity pass below only sees genuine geometry problems.
    snapped = shapely.set_precision(geoms.values, grid_size)

    # 2. make_valid, never buffer(0): it splits a bowtie into its parts instead
    #    of discarding the smaller one.
    repaired = shapely.make_valid(snapped)

    out = gdf.copy()
    out.geometry = gpd.GeoSeries(repaired, crs=gdf.crs, index=gdf.index)
    out = out[~out.geometry.is_empty].copy()

    area_after = float(out.geometry.area.sum())
    report = RepairReport(
        n_features=len(gdf),
        n_invalid_in=n_invalid_in,
        n_repaired=n_invalid_in - int((~out.geometry.is_valid).sum()),
        area_before_m2=area_before,
        area_after_m2=area_after,
    )

    if not bool(out.geometry.is_valid.all()):
        raise ValueError("repair left invalid geometries; inspect them individually")

    if report.relative_area_loss > max_area_loss:
        raise ValueError(
            f"repair removed {report.relative_area_loss:.4%} of layer area "
            f"({area_before - area_after:,.0f} m²), above the "
            f"{max_area_loss:.4%} tolerance. This is a data problem, not a "
            "geometry problem — review the failing features before proceeding."
        )

    logger.info(
        "repaired %d/%d invalid feature(s); area delta %.6f%%",
        report.n_repaired, report.n_invalid_in, report.relative_area_loss * 100.0,
    )
    return out, report

Two design choices in that function are deliberate and worth stating. It raises on excessive area loss rather than warning, because a repair that deletes a square kilometre of coast is not a repair, and a warning in a log is not a control. And it measures area before repair using a repaired copy, because the area of an invalid polygon is undefined in GEOS and comparing it to anything is meaningless — the copy is thrown away immediately and exists only so the two numbers are commensurable.

Simplification Without Breaking the Topology #

Coastal polygons are dense. A national shoreline digitised at 1:10,000 carries far more vertices than a 50-metre analysis grid can use, and carrying them costs time in every predicate and every render. Simplification is the obvious answer and the obvious implementation is wrong.

Per-feature Douglas-Peucker — gdf.geometry.simplify(tolerance) — treats each polygon independently. Where two polygons share an edge, as adjacent management zones or a shoreline and the estuary polygon inside it always do, the two copies of that edge are simplified with different vertex neighbourhoods and no longer coincide. The result is a layer full of slivers and gaps a few metres wide, invisible at the zoom level anyone reviews it at, and fatal to every subsequent overlay: a point-in-polygon test in a sliver returns no match, and a zonal statistic over the pair double-counts the overlap.

Topology-preserving simplification builds a shared-edge graph first, simplifies each edge once, and reassembles the polygons from the simplified edges. mapshaper, PostGIS’s ST_SimplifyPreserveTopology and GRASS v.generalize all implement it; Shapely does not, which is why the correct answer to “how do I simplify a coastline in Shapely” is usually “not in Shapely”.

The other constraint is the tolerance itself. It has to be chosen against the analysis, not against file size: a tolerance of 25 metres is negligible against a 100-metre habitat grid and catastrophic against a 5-metre dredging boundary, and the same layer often feeds both. Where that happens, publish the layer at two tolerances with the tolerance in the name, rather than picking a compromise that is wrong for both. The shoreline simplification guide works through both halves with runnable code.

Validation Gates #

Three checks, all cheap, all run before a vector layer is allowed downstream.

Gate 1 — validity and emptiness. Every geometry is_valid, no geometry empty, no null geometries. This is a one-liner and it catches the majority of problems that would otherwise surface as an exception in an unrelated stage an hour later.

def assert_layer_clean(gdf: gpd.GeoDataFrame) -> None:
    if gdf.geometry.isna().any():
        raise ValueError(f"{int(gdf.geometry.isna().sum())} null geometries")
    if gdf.geometry.is_empty.any():
        raise ValueError(f"{int(gdf.geometry.is_empty.sum())} empty geometries")
    bad = gdf.index[~gdf.geometry.is_valid].tolist()
    if bad:
        raise ValueError(f"{len(bad)} invalid geometries at index {bad[:10]}")

Gate 2 — extent sanity. The layer’s bounds fall inside the survey’s expected extent, in the working CRS. This is the vector equivalent of the bounds assertion that catches axis-order inversion, and it catches exactly the same class of error: a layer that reprojected “successfully” into the wrong hemisphere.

Gate 3 — self-overlap. Within a layer that is supposed to partition space — management zones, survey blocks, habitat classes — no two features may overlap by more than a tolerance. gdf.sindex.query(gdf.geometry, predicate="overlaps") gives the candidate pairs in one call, and the total overlapping area is the number to threshold on. A layer that fails this gate will produce zonal statistics that sum to more than the whole, which is the kind of error that gets noticed in a briefing rather than in a log.

Spatial Indexing Is the Difference Between Seconds and Hours #

Every coastal workflow eventually joins one layer against another: soundings against survey blocks, AIS positions against port polygons, habitat cells against management zones. Done naively that is a nested loop, and at coastal data volumes the nested loop does not finish.

Shapely 2.0 and GeoPandas expose an R-tree through gdf.sindex, and the important part of the API is query with a predicate argument. Without a predicate it returns bounding-box candidates, which still need refining; with one — "intersects", "within", "contains" — it returns the exact pairs, performing the bounding-box filter internally and the exact predicate only on survivors. The cost difference against a manual loop is three to four orders of magnitude on a million-point join, and the code is shorter:

import numpy as np


def points_in_zones(points: gpd.GeoDataFrame, zones: gpd.GeoDataFrame) -> np.ndarray:
    """Zone index for every point, -1 where a point falls outside all zones."""
    if points.crs != zones.crs:
        raise ValueError(f"CRS mismatch: {points.crs} vs {zones.crs}")

    tree_idx, point_idx = zones.sindex.query(points.geometry, predicate="within")
    assignment = np.full(len(points), -1, dtype="int64")
    assignment[tree_idx] = point_idx
    logger.info("assigned %d/%d points to a zone", int((assignment >= 0).sum()), len(points))
    return assignment

Two details bite people. The index is built lazily on first access and cached on the GeoDataFrame, so building it inside a loop over chunks rebuilds it every iteration — build it once, outside. And the returned arrays are index positions, not DataFrame index labels; on a layer with a non-default index, using them as labels produces a silent misalignment rather than an error.

For repeated joins against a layer that does not change — a national shoreline, a fixed set of management zones — the index can be persisted by keeping the GeoDataFrame alive across the run, which is the single largest speedup available in most coastal pipelines and requires no new dependency.

Common Failure Modes #

A GeoDataFrame with a valid CRS and geometries in a different one. Assigning gdf.crs labels; it does not transform. A layer relabelled rather than reprojected passes every validity check and sits thousands of kilometres from where it belongs. Guard with Gate 2.

Mixed geometry types in one column. A shapefile of “coastline” containing both Polygon and MultiPolygon, or worse, a stray LineString from a digitising error. Most operations handle the first pair transparently and then one does not, usually a format writer. Normalise to multi-part on ingestion.

unary_union on a large layer, in one call. Memory grows super-linearly with the number of parts and the operation is not interruptible. Union in batches — union each thousand features, then union the results — which is both faster and bounded in memory.

Silent CRS loss through a pandas operation. pd.concat on GeoDataFrames returns a DataFrame with an object column, not a GeoDataFrame, and the CRS is gone. The result still plots, because matplotlib does not care, and fails at the next spatial join. Re-wrap explicitly after any pandas-level operation.

Pipeline Integration and Downstream Handoff #

The output of this stage is a GeoPackage per layer, in the working projected CRS, with the shoreline datum and epoch recorded in the layer metadata and in the filename. Alongside it goes the RepairReport as JSON — feature counts, how many were repaired, the area delta — which becomes part of the run manifest described in the parent section. Downstream, the masking stage in bathymetric processing consumes the shoreline directly, and the AIS density products consume the same layer to clip traffic rasters to water.

Publish the vertex count and the simplification tolerance in that metadata too. They are the two numbers a downstream consumer needs to decide whether this layer is precise enough for their analysis, and they are the two numbers that are impossible to recover from the geometry alone once several layers have been merged.

Related #

• Repairing Invalid Coastline Polygons in Shapely — the three invalidity classes, diagnosed and repaired individually
• Simplifying Shorelines Without Breaking Topology — shared-edge simplification and choosing a tolerance
• Clipping Marine Rasters to a Coastline Mask — the raster side of the same boundary
• CRS Alignment for Coastal GIS Projects — the reprojection stage this topic assumes has already run
• Tidal Datum Transformations in Python — where the datum a shoreline is defined against comes from

Up: Marine Spatial Data Fundamentals & Architecture