Delineating Nested Catchments from Multiple Pour Points

When several gauges sit along the same river, the catchment above the downstream gauge fully contains the catchment above each upstream gauge, and that nesting is exactly the structure a routing model needs. This guide delineates every catchment from a shared flow-direction grid and then rebuilds the containment tree in Python, as a focused reference within the Nested Catchment Delineation topic, which sits inside the wider Watershed Delineation & Catchment Synchronization coverage on this site.

The efficient pattern is to condition the DEM and compute flow direction once, then trace a catchment from each pour point against that single grid. Delineating each outlet against its own freshly computed direction grid is both slower and dangerous: tiny numerical differences produce divides that do not nest cleanly, so the downstream catchment fails to contain its upstream siblings. One grid, many traversals, then a spatial join to recover the hierarchy is the reliable route.

The two-phase workflow

Delineation and hierarchy reconstruction are separate concerns. The first phase is raster: mask cells that drain to each outlet. The second phase is vector: discover which mask contains which.

The concentric result on the right is the whole objective: a strict nesting where P3’s catchment lies inside P2’s, which lies inside P1’s. That structure lets a model compute incremental drainage — the area between two nested outlets — simply by subtracting the inner geometry from the outer. Reconstructing it correctly depends on getting both the delineation and the containment test right, and the same synchronisation concerns arise when the nested divides come from DEMs of different resolution, covered in synchronizing nested catchment boundaries across DEM resolutions.

Prerequisites

conda create -n nested python=3.11 pysheds geopandas rasterio shapely numpy -c conda-forge
conda activate nested

The DEM must already be hydrologically conditioned. If sinks remain, catchments leak across false divides. Fill depressions following best practices for filling sinks in high-resolution LiDAR data before you begin, and confirm your pour points are close to the channel so snapping succeeds.

Phase 1: delineate every catchment from one grid

pysheds computes flow direction once and then traces a catchment per pour point. Snapping each point to the nearest high-accumulation cell before tracing is what keeps a gauge from delineating a hillslope sliver instead of the true basin.

import logging
import numpy as np
import geopandas as gpd
from pysheds.grid import Grid
from shapely.geometry import shape

logger = logging.getLogger(__name__)


def delineate_pour_points(dem_path: str, pour_points: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Delineate one catchment per pour point against a single flow-direction grid."""
    grid = Grid.from_raster(dem_path)
    dem = grid.read_raster(dem_path)

    # Condition once, reuse for all points
    filled = grid.fill_depressions(dem)
    flooded = grid.resolve_flats(filled)
    fdir = grid.flowdir(flooded)
    acc = grid.accumulation(fdir)
    logger.info("Flow direction and accumulation computed; %d pour points to process", len(pour_points))

    records = []
    for pid, point in pour_points.iterrows():
        # Snap to the nearest large-accumulation cell so we trace the channel, not a hillslope
        x_snap, y_snap = grid.snap_to_mask(acc > 1000, (point.geometry.x, point.geometry.y))
        catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, xytype="coordinate")
        grid.clip_to(catch)
        geoms = [shape(geom) for geom, val in grid.polygonize() if val == 1]
        if not geoms:
            logger.warning("Pour point %s produced no catchment cells; check snapping", pid)
            continue
        merged = max(geoms, key=lambda g: g.area)  # dominant polygon
        records.append({"point_id": pid, "geometry": merged, "area_m2": merged.area})
        logger.info("Pour point %s -> catchment area %.1f m^2", pid, merged.area)
        grid.clip_to(dem)  # reset clip for the next point

    return gpd.GeoDataFrame(records, crs=pour_points.crs)

Because flow direction is computed a single time, every catchment is traced against identical divides. That consistency is precisely what makes the downstream containment test reliable in phase 2. If you prefer the RichDEM toolchain, the equivalent is to compute one D8 direction grid and label upstream cells per outlet; the nesting logic that follows is identical.

Phase 2: reconstruct the nesting hierarchy

With every catchment vectorised, a spatial join using the contains predicate reveals which catchment encloses which. The direct parent of a catchment is the smallest catchment that contains it.

import logging
import geopandas as gpd

logger = logging.getLogger(__name__)


def build_nesting(catchments: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Assign each catchment its direct parent: the smallest catchment containing it."""
    # Test containment on representative points to avoid shared-boundary false negatives
    pts = catchments.copy()
    pts["geometry"] = catchments.geometry.representative_point()

    joined = gpd.sjoin(
        pts, catchments, how="left", predicate="within", lsuffix="child", rsuffix="parent"
    )
    joined = joined[joined.index != joined.index_parent]  # drop self-containment
    logger.info("Containment join produced %d enclosing relationships", len(joined))

    parent_of = {}
    for child_id, group in joined.groupby(level=0):
        # Smallest enclosing catchment is the direct parent
        enclosing = catchments.loc[group.index_parent.values]
        if enclosing.empty:
            parent_of[child_id] = -1  # a root catchment
            continue
        direct = enclosing.area.idxmin()
        parent_of[child_id] = direct
        logger.debug("Catchment %s -> parent %s", child_id, direct)

    result = catchments.copy()
    result["parent_id"] = result.index.map(lambda i: parent_of.get(i, -1))
    n_roots = int((result.parent_id == -1).sum())
    logger.info("Hierarchy built: %d catchments, %d roots", len(result), n_roots)
    return result


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    points = gpd.read_file("gauges.gpkg").to_crs(5070)
    logger.info("Loaded %d pour points", len(points))
    catchments = delineate_pour_points("conditioned_dem.tif", points)
    nested = build_nesting(catchments)
    nested.to_file("nested_catchments.gpkg", driver="GPKG")

Testing containment on each catchment’s representative_point instead of its full polygon sidesteps the shared-outlet problem: nested catchments share their downstream boundary, so a strict polygon-in-polygon contains can return False along that coincident edge. A representative point is guaranteed interior, so it lands unambiguously inside every true ancestor.

Worked example: incremental drainage areas

Three gauges sit along one river: an upstream headwater gauge P3, a midstream gauge P2, and a downstream gauge P1. Delineation returns 42.1 km^2 for P3, 118.6 km^2 for P2, and 305.4 km^2 for P1. The containment join assigns P3 a parent of P2 and P2 a parent of P1, giving the nesting P1 contains P2 contains P3.

The incremental areas fall straight out of the hierarchy by subtraction. The local drainage that enters between P3 and P2 is 118.6 - 42.1 = 76.5 km^2, and the local drainage between P2 and P1 is 305.4 - 118.6 = 186.8 km^2. These increments are what a routing model actually needs: they represent the ungauged area contributing to each reach, and they must sum with the headwater area back to the total, 42.1 + 76.5 + 186.8 = 305.4 km^2, a closure check worth asserting in code. If the increments do not close, either a catchment failed to nest or a pour point snapped to the wrong channel. The nesting tree therefore does double duty: it orders the catchments for routing and it provides a built-in consistency check on the delineation, both from a single spatial join rather than repeated graph walks over the flow-direction grid.

Hierarchy field reference

Gotchas & edge cases

• Unsnapped pour points delineate the wrong basin. A gauge coordinate a few cells off the channel traces the tiny hillslope catchment it happens to land on. Always snap to a high-accumulation cell and log the snap distance so large moves are visible.
• Coincident outlets break strict contains. Nested catchments share a downstream edge, so polygon containment can fail there. Use representative points or centroids for the join, never raw geometries.
• Non-nested gauges on separate tributaries have no parent. Two gauges on different tributaries of the same river produce sibling catchments that do not contain each other; both are roots until a downstream gauge encloses them. A -1 parent is legitimate, not an error.
• Polygonising leaves stray one-cell fragments. Tracing can emit isolated diagonal cells as separate polygons. Keep the largest-area polygon per pour point, or dissolve fragments that touch the main mask, before building the hierarchy.

Related Topics

• Nested Catchment Delineation — parent topic covering how nested drainage areas are defined and used in modelling
• Synchronizing Nested Catchment Boundaries Across DEM Resolutions — how to reconcile the same nested divides delineated from DEMs of differing resolution
• Watershed Delineation & Catchment Synchronization — the overview tying delineation, nesting, and reconciliation together