Validating Catchment Boundary Topology with GeoPandas and Shapely

A set of delineated catchment polygons is only trustworthy once you can prove three things about it: every polygon is individually valid, no two polygons overlap, and together they leave no gaps inside the modelled area. This guide builds those three proofs in Python, serving as a focused reference under the Boundary Topology Validation topic, which belongs to the wider Watershed Delineation & Catchment Synchronization coverage on this site.

Vector catchments produced by tracing a flow-accumulation raster inherit the staircase geometry of the grid, and every stair step is an opportunity for a self-touching ring, a duplicated vertex, or a hairline gap against a neighbour. Left unchecked these defects propagate: an overlap double-counts rainfall, a gap loses runoff entirely, and an invalid ring crashes the first spatial join that touches it. Catching them early is far cheaper than debugging a mass-balance discrepancy later.

Two distinct classes of defect

It helps to separate the problem into geometry validity, which is a property of a single feature, and coverage topology, which is a relationship between features.

The distinction matters because the tools differ. Validity is tested per feature with shapely.is_valid and repaired with shapely.make_valid. Coverage is tested across features with a geopandas overlay for overlaps and a unary_union for gaps. A dataset can pass every validity check and still be riddled with overlaps, so both stages are mandatory. Slivers — the narrow overlap and gap polygons that dominate raster-traced coverages — get their own dedicated treatment in detecting and removing sliver polygons in catchment boundaries; this page focuses on detecting the full spectrum of defects first.

Prerequisites

conda create -n topo python=3.11 geopandas shapely pyproj numpy -c conda-forge
conda activate topo

Work in a projected CRS with metre units so that area and buffer tolerances are physically meaningful. If your catchment shapefiles arrive in inconsistent projections, reconcile them first following fixing CRS mismatches in watershed shapefiles. Shapely 2.0 or newer is required for the vectorised make_valid.

Stage 1: per-feature validity

The first pass flags and explains every invalid geometry. explain_validity returns a human-readable reason and the offending coordinate, which is invaluable when triaging hundreds of catchments.

import logging
import geopandas as gpd
from shapely.validation import explain_validity, make_valid

logger = logging.getLogger(__name__)


def validate_and_repair(catchments: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Flag invalid catchment geometries, repair them, and log area drift."""
    invalid_mask = ~catchments.geometry.is_valid
    n_invalid = int(invalid_mask.sum())
    logger.info("Validity scan: %d of %d catchments are invalid", n_invalid, len(catchments))

    repaired = catchments.copy()
    for idx in catchments.index[invalid_mask]:
        geom = catchments.geometry.loc[idx]
        reason = explain_validity(geom)
        area_before = geom.area
        fixed = make_valid(geom)
        area_after = fixed.area
        drift = abs(area_after - area_before)
        logger.warning(
            "Catchment %s invalid (%s); area drift after make_valid = %.1f m^2",
            idx, reason, drift
        )
        repaired.loc[idx, "geometry"] = fixed

    still_bad = int((~repaired.geometry.is_valid).sum())
    logger.info("After repair: %d catchments remain invalid", still_bad)
    return repaired

A repair that changes area by more than a couple of cells is a red flag: make_valid closed a self-intersection by discarding a lobe, which usually means the delineation itself produced a bowtie where the boundary crossed a flat divide. Investigate those cases rather than silently accepting the repair.

Stage 2: detecting overlaps

Two valid catchments must not share interior area. A geopandas self-overlay with how="intersection" returns exactly the shared regions; anything with positive area is an overlap to resolve.

import logging
import geopandas as gpd

logger = logging.getLogger(__name__)


def find_overlaps(catchments: gpd.GeoDataFrame, area_tol: float = 1.0) -> gpd.GeoDataFrame:
    """Return overlap polygons between distinct catchments above area_tol (m^2)."""
    left = catchments.reset_index(names="cid_left")
    right = catchments.reset_index(names="cid_right")
    inter = gpd.overlay(left, right, how="intersection", keep_geom_type=True)

    # Drop self-pairs and mirror duplicates
    inter = inter[inter.cid_left < inter.cid_right].copy()
    inter["overlap_area"] = inter.geometry.area
    overlaps = inter[inter.overlap_area > area_tol]
    logger.info(
        "Overlap scan: %d catchment pairs share more than %.1f m^2",
        len(overlaps), area_tol
    )
    if not overlaps.empty:
        worst = overlaps.overlap_area.max()
        logger.warning("Largest overlap = %.1f m^2", worst)
    return overlaps

keep_geom_type=True is important: without it, catchments that merely touch along an edge produce zero-area LineString results that pollute the overlap set. Filtering by an area tolerance derived from the DEM cell size then separates real overlaps from touching-edge noise.

Stage 3: detecting gaps and asserting coverage

Gaps are the inverse problem. Dissolve the whole coverage into one geometry with unary_union, then compare it against the region the coverage is supposed to fill. Any difference is a hole.

import logging
import geopandas as gpd
from shapely.ops import unary_union

logger = logging.getLogger(__name__)


def find_gaps(catchments: gpd.GeoDataFrame, study_area, area_tol: float = 1.0) -> gpd.GeoDataFrame:
    """Return gap polygons inside study_area not covered by any catchment."""
    dissolved = unary_union(catchments.geometry.values)
    holes = study_area.difference(dissolved)

    gap_gdf = gpd.GeoDataFrame(geometry=[holes], crs=catchments.crs).explode(index_parts=False)
    gap_gdf["gap_area"] = gap_gdf.geometry.area
    gap_gdf = gap_gdf[gap_gdf.gap_area > area_tol]

    summed = float(catchments.geometry.area.sum())
    covered = float(dissolved.area)
    logger.info("Coverage: summed=%.1f m^2, dissolved=%.1f m^2", summed, covered)
    if summed - covered > area_tol:
        logger.warning(
            "Summed area exceeds dissolved area by %.1f m^2 -> overlaps present",
            summed - covered
        )
    logger.info("Gap scan: %d holes above %.1f m^2", len(gap_gdf), area_tol)
    return gap_gdf


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    gdf = gpd.read_file("catchments.gpkg").to_crs(5070)
    logger.info("Loaded %d catchments", len(gdf))
    repaired = validate_and_repair(gdf)
    overlaps = find_overlaps(repaired)
    study = unary_union(repaired.geometry.values).convex_hull
    gaps = find_gaps(repaired, study)

The mass-balance check inside find_gaps is the cheapest topology test you can run: if the summed area of the individual catchments exceeds the area of their dissolved union, some polygons must be overlapping, because overlapping regions are counted twice in the sum but once in the union. This single comparison catches overlap problems even before you run the full overlay.

Worked example: reading a mass-balance discrepancy

Suppose a coverage of 412 catchments reports a summed area of 5,004.2 km^2 while the dissolved union measures 4,998.7 km^2. The 5.5 km^2 excess is the fingerprint of overlaps, because overlapping regions are counted in the per-feature sum but only once in the union. Running find_overlaps on the same coverage returns eleven overlap polygons, of which nine are hairline slivers under 50 m^2 and two are substantial: a 3.1 km^2 and a 2.3 km^2 shared region between adjacent catchments whose divides were traced from two different flow-direction tiles that met along a mosaic seam.

The nine slivers are rasterisation noise, resolved by snapping coordinates to the DEM grid with shapely.set_precision before re-dissolving. The two large overlaps are real errors: they trace back to a tiling boundary where the flow-direction grids disagreed by a few cells. The fix is upstream — re-route the seam tile with a shared buffer — but the topology check is what surfaced a defect that no per-feature validity test would ever have caught. That is the practical value of the mass-balance comparison: a single subtraction pointed straight at a two-kilometre error hiding inside a valid dataset.

Assertion reference

Gotchas & edge cases

• unary_union is not free on large coverages. Dissolving tens of thousands of high-vertex catchments is memory-heavy. Tile the study area and dissolve per tile, then stitch, rather than dissolving a national coverage in one call.
• Floating-point coordinates create phantom slivers. Two catchments traced from the same raster edge can differ by nanometres and produce a hairline overlap. Snap coordinates to a grid with shapely.set_precision before the overlay so identical edges compare equal.
• A valid MultiPolygon can hide a topology error. make_valid may return a MultiPolygon whose parts are disjoint. If a catchment should be a single connected area, assert geom_type == "Polygon" and treat any MultiPolygon as suspect.
• The study-area boundary defines what counts as a gap. Using the convex hull as the reference invents gaps in concave coastlines. Prefer the true modelled boundary, or the dissolved union buffered inward, when concavity matters.

Related Topics

• Boundary Topology Validation — parent topic covering the full validation workflow for delineated catchment coverages
• Detecting and Removing Sliver Polygons in Catchment Boundaries — the follow-on that measures and eliminates the narrow slivers this page detects
• Watershed Delineation & Catchment Synchronization — the overview connecting delineation, validation, and reconciliation