How do I handle outlets on the same stream segment?

Assign each outlet a unique reach ID and accumulation value. The upstream member will have a lower accumulation value and will be processed later in the sorted queue; its contributing area is then subtracted from the already-delineated downstream basin, producing the correct incremental drainage area between the two gauges.

Nested Catchment Delineation: Automated Python Workflows

Drainage networks are hierarchical by nature — every stream has tributaries, every gauging station sits inside a larger contributing basin, and flood forecasts depend on resolving those overlapping areas without double-counting. Nested catchment delineation is the computational workflow that extracts those hierarchical, non-overlapping drainage units programmatically from a flow-routed DEM. It is a core step within Watershed Delineation & Catchment Synchronization and produces the spatially consistent sub-basin geometries that distributed hydrologic models, water quality assessments, and infrastructure siting tools depend on.

Unlike flat basin partitioning (covered in basin partitioning strategies), nested delineation must maintain explicit parent-child topology: the full upstream area of a downstream outlet must contain all upstream sub-basins, and no cell may be claimed by more than one sub-basin at the same level. This topological constraint is the central engineering challenge. It demands strict downstream-to-upstream processing, outlet snapping to channel cells, and vector-level difference operations to subtract already-assigned areas. Getting any of these steps wrong produces overlapping polygons, missing slivers, or inflated contributing areas that propagate silently through every downstream model.

This guide provides a production-tested Python implementation, from hydro-conditioning through validation, together with an analysis of common failure modes and when to prefer alternative partitioning approaches.

Prerequisites & Environment Setup

Python stack

bash

pip install pysheds>=0.3.5 geopandas>=0.13.0 rasterio>=1.3.0 shapely>=2.0.0 numpy pandas

Input data specifications

Hydrologically conditioned DEM: 30 m or finer, sink-filled or breach-conditioned, FLOAT32, no-data = -9999. Unconditioned DEMs produce artificial sinks and fragmented flow paths. Apply DEM pit filling algorithms before this workflow.
Outlet points: CSV or GeoJSON with id, x, y columns. Coordinates must match the DEM’s CRS — ensure coordinate reference system alignment before snapping.
CRS: All layers in a metric equal-area or UTM projection (e.g., EPSG:32610). Geographic coordinates (WGS84 decimal degrees) introduce area distortion that breaks topological checks.
Memory: 16 GB RAM for domains up to ~1 000 km² at 30 m; 32 GB or out-of-core processing for larger extents.

Algorithm Mechanics

The non-overlap invariant

Nested catchment delineation must satisfy one invariant: the set of all sub-basin geometries at a given level is a partition of the parent basin — no gaps, no overlaps. This is enforced by the sequential difference operation. For each outlet processed after the first, the new upstream polygon is differenced against the union of all previously extracted polygons before it is recorded.

Formally, let A(i) be the full upstream contributing area of outlet i, and let P(k) be the incremental polygon assigned to outlet k. For outlets sorted from most-downstream (highest accumulation) to most-upstream:

text

P(k) = A(k)  –  UNION( P(1), P(2), ..., P(k-1) )

Each P(k) is the portion of the drainage basin that drains to outlet k but not to any already-processed downstream outlet.

D8 vs other routing models

D8 flow direction is strongly preferred for nested delineation because it assigns every cell to exactly one downslope neighbor, producing fully partitioned drainage basins with integer-valued boundaries. Multi-directional models such as D-Infinity routing distribute flow proportionally across two neighbors; this is physically realistic for hillslope processes but makes it impossible to assign a cell definitively to one sub-basin, creating fractional area leakage at boundaries.

Outlet snapping to stream cells

DEM resolution and minor GPS offsets mean that a nominal outlet coordinate rarely lands on the highest-accumulation cell in its neighbourhood. The snapping step searches a window of radius r cells centred on the nominal coordinate and relocates the outlet to the cell with the maximum accumulation value. Snapping radius selection is a balance: too small misses genuine channel cells; too large relocates the outlet to an adjacent tributary. For 30 m DEMs a radius of 5 cells (150 m) is a reasonable default; for 10 m LiDAR, 3 cells (30 m) is usually sufficient.

Parameter / encoding reference

Parameter	Typical value	Effect if wrong
`snap_radius_px`	3–7 cells	Too small: outlet stays on hillslope; too large: migrates to wrong tributary
`min_area_km2`	0.01	Below this, polygon fragments from raster edge effects are dropped
`area_tolerance_pct`	2 %	Maximum acceptable discrepancy between raster mask area and vector polygon area
D8 encoding (pysheds)	Powers of 2 (1=E … 128=NE)	Wrong encoding → entire flow direction grid misrouted
Flow fill strategy	`fill_depressions` then `resolve_flats`	Skipping resolve_flats leaves flat areas with undefined direction

Step-by-Step Workflow

Step 1 — Hydro-condition the DEM and compute flow routing

python

import logging
import numpy as np
from pysheds.grid import Grid

logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
logger = logging.getLogger(__name__)

logger.info("Loading hydrologically conditioned DEM")
grid = Grid.from_raster("dem_conditioned.tif")
dem = grid.read_raster("dem_conditioned.tif")

logger.info("Filling pits, depressions, and resolving flat areas")
pit_filled  = grid.fill_pits(dem)
flooded     = grid.fill_depressions(pit_filled)
inflated    = grid.resolve_flats(flooded)

logger.info("Computing D8 flow direction and accumulation")
fdir = grid.flowdir(inflated)       # pysheds power-of-2 encoding
acc  = grid.accumulation(fdir)

logger.info("Flow accumulation max cell = %.0f", float(acc.max()))

Step 2 — Load and snap outlet points

python

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

logger.info("Loading outlet coordinates")
outlets_df  = pd.read_csv("outlets.csv")
outlets_gdf = gpd.GeoDataFrame(
    outlets_df,
    geometry=gpd.points_from_xy(outlets_df["x"], outlets_df["y"]),
    crs="EPSG:32610",
)

acc_arr = np.array(acc)

def snap_to_stream(x: float, y: float, search_radius_px: int = 5) -> tuple[float, float]:
    """Snap a coordinate to the highest-accumulation cell within a pixel radius."""
    col_f, row_f = ~grid.affine * (x, y)
    row, col = int(row_f), int(col_f)
    r0 = max(0, row - search_radius_px)
    r1 = min(acc_arr.shape[0], row + search_radius_px + 1)
    c0 = max(0, col - search_radius_px)
    c1 = min(acc_arr.shape[1], col + search_radius_px + 1)
    window    = acc_arr[r0:r1, c0:c1]
    local_idx = np.unravel_index(np.argmax(window), window.shape)
    best_row, best_col = r0 + local_idx[0], c0 + local_idx[1]
    sx, sy = grid.affine * (best_col + 0.5, best_row + 0.5)
    return float(sx), float(sy)

logger.info("Snapping %d outlets to stream cells", len(outlets_gdf))
outlets_gdf[["snap_x", "snap_y"]] = outlets_gdf.apply(
    lambda r: pd.Series(snap_to_stream(r.geometry.x, r.geometry.y)), axis=1
)

# Retrieve accumulation value at snapped position for sort key
def acc_at(x: float, y: float) -> float:
    c, r = ~grid.affine * (x, y)
    return float(acc_arr[int(r), int(c)])

outlets_gdf["acc_value"] = outlets_gdf.apply(
    lambda r: acc_at(r["snap_x"], r["snap_y"]), axis=1
)

# Downstream-first order: highest accumulation value processed first
outlets_sorted = outlets_gdf.sort_values("acc_value", ascending=False).reset_index(drop=True)
logger.info("Outlets sorted; most-downstream outlet has acc = %.0f", outlets_sorted["acc_value"].iloc[0])

Step 3 — Iterative upstream extraction with topological subtraction

python

from rasterio.features import shapes as rio_shapes
from shapely.geometry import shape
from shapely.validation import make_valid

nested_catchments: list[dict] = []
processed_union = None

for idx, outlet in outlets_sorted.iterrows():
    oid   = outlet["id"]
    sx, sy = outlet["snap_x"], outlet["snap_y"]

    logger.info("Extracting upstream area for outlet %s (acc=%.0f)", oid, outlet["acc_value"])
    catch     = grid.catchment(x=sx, y=sy, fdir=fdir, xytype="coordinate")
    catch_arr = grid.view(catch, dtype=np.uint8)

    polys = [shape(s) for s, v in rio_shapes(catch_arr, transform=grid.affine) if v == 1]
    if not polys:
        logger.warning("Outlet %s produced no upstream cells; skipping", oid)
        continue

    raw_poly = polys[0] if len(polys) == 1 else gpd.GeoSeries(polys).unary_union
    raw_poly = make_valid(raw_poly)

    # Remove areas already assigned to downstream outlets
    if processed_union is not None:
        incremental = raw_poly.difference(processed_union)
        incremental = make_valid(incremental)
    else:
        incremental = raw_poly

    area_km2 = incremental.area / 1e6
    if incremental.is_empty or area_km2 < 0.01:
        logger.warning("Outlet %s incremental area too small (%.4f km²); skipping", oid, area_km2)
        continue

    nested_catchments.append({
        "id":       oid,
        "geometry": incremental,
        "area_km2": area_km2,
        "snap_x":   sx,
        "snap_y":   sy,
    })
    processed_union = (
        incremental if processed_union is None else make_valid(processed_union.union(incremental))
    )
    logger.info("Outlet %s → %.2f km² incremental area", oid, area_km2)

Step 4 — Vectorize and export

python

import geopandas as gpd

logger.info("Compiling %d nested catchments", len(nested_catchments))
final_gdf = gpd.GeoDataFrame(nested_catchments, crs="EPSG:32610")
final_gdf = final_gdf.sort_values("area_km2", ascending=False).reset_index(drop=True)
final_gdf["hierarchy_level"] = final_gdf.index + 1  # 1 = outermost/largest

output_path = "nested_catchments.gpkg"
final_gdf.to_file(output_path, driver="GPKG")
logger.info("Wrote %s — %d features, total area %.2f km²",
            output_path, len(final_gdf), final_gdf["area_km2"].sum())

Production-Ready Code

The function below integrates all steps with logging, error handling, and metadata export. Drop it into any pipeline that supplies a conditioned DEM and an outlet CSV.

python

import logging
import numpy as np
import pandas as pd
import geopandas as gpd
from pathlib import Path
from shapely.geometry import shape
from shapely.validation import make_valid
from rasterio.features import shapes as rio_shapes
from pysheds.grid import Grid

logger = logging.getLogger(__name__)


def delineate_nested_catchments(
    dem_path: str,
    outlets_csv: str,
    output_gpkg: str,
    epsg: int = 32610,
    snap_radius_px: int = 5,
    min_area_km2: float = 0.01,
) -> gpd.GeoDataFrame:
    """
    Delineate non-overlapping nested catchments from a conditioned DEM and outlet CSV.

    Parameters
    ----------
    dem_path : str
        Path to hydrologically conditioned DEM (FLOAT32, no-data = -9999).
    outlets_csv : str
        CSV with columns: id, x, y (in the same CRS as the DEM).
    output_gpkg : str
        Output GeoPackage path.
    epsg : int
        EPSG code for all spatial data (metric CRS required).
    snap_radius_px : int
        Search radius in pixels for outlet snapping to stream cells.
    min_area_km2 : float
        Minimum incremental area; polygons below this are discarded as artefacts.

    Returns
    -------
    GeoDataFrame with columns: id, geometry, area_km2, snap_x, snap_y, hierarchy_level.
    """
    # ── 1. Hydro-condition and route ──────────────────────────────────────────
    logger.info("Loading DEM: %s", dem_path)
    grid = Grid.from_raster(dem_path)
    dem  = grid.read_raster(dem_path)

    logger.info("Conditioning: fill pits → fill depressions → resolve flats")
    inflated = grid.resolve_flats(grid.fill_depressions(grid.fill_pits(dem)))

    logger.info("Computing D8 flow direction and accumulation")
    fdir = grid.flowdir(inflated)
    acc  = grid.accumulation(fdir)
    acc_arr = np.array(acc)
    logger.info("Accumulation range: %.0f – %.0f cells", acc_arr.min(), acc_arr.max())

    # ── 2. Snap outlets ───────────────────────────────────────────────────────
    df  = pd.read_csv(outlets_csv)
    gdf = gpd.GeoDataFrame(
        df,
        geometry=gpd.points_from_xy(df["x"], df["y"]),
        crs=f"EPSG:{epsg}",
    )

    def _snap(x: float, y: float) -> tuple[float, float, float]:
        c, r = ~grid.affine * (x, y)
        row, col = int(r), int(c)
        r0, r1 = max(0, row - snap_radius_px), min(acc_arr.shape[0], row + snap_radius_px + 1)
        c0, c1 = max(0, col - snap_radius_px), min(acc_arr.shape[1], col + snap_radius_px + 1)
        win = acc_arr[r0:r1, c0:c1]
        li  = np.unravel_index(np.argmax(win), win.shape)
        br, bc = r0 + li[0], c0 + li[1]
        sx, sy = grid.affine * (bc + 0.5, br + 0.5)
        return float(sx), float(sy), float(acc_arr[br, bc])

    logger.info("Snapping %d outlets (radius=%d px)", len(gdf), snap_radius_px)
    snapped = gdf.apply(lambda row: pd.Series(_snap(row.geometry.x, row.geometry.y),
                                               index=["snap_x", "snap_y", "acc_value"]), axis=1)
    gdf = pd.concat([gdf, snapped], axis=1)
    gdf = gdf.sort_values("acc_value", ascending=False).reset_index(drop=True)

    # ── 3. Iterative extraction with topological subtraction ──────────────────
    results: list[dict] = []
    cumulative_union = None

    for _, outlet in gdf.iterrows():
        oid = outlet["id"]
        sx, sy = outlet["snap_x"], outlet["snap_y"]
        logger.info("Processing outlet %s (acc=%.0f)", oid, outlet["acc_value"])

        catch_mask = grid.view(
            grid.catchment(x=sx, y=sy, fdir=fdir, xytype="coordinate"),
            dtype=np.uint8,
        )
        polys = [shape(s) for s, v in rio_shapes(catch_mask, transform=grid.affine) if v == 1]
        if not polys:
            logger.warning("Outlet %s: no upstream cells found; skipping", oid)
            continue

        raw = polys[0] if len(polys) == 1 else gpd.GeoSeries(polys).unary_union
        raw = make_valid(raw)

        incremental = raw.difference(cumulative_union) if cumulative_union else raw
        incremental = make_valid(incremental)

        area_km2 = incremental.area / 1e6
        if incremental.is_empty or area_km2 < min_area_km2:
            logger.warning("Outlet %s: incremental area %.5f km² below threshold; skipping", oid, area_km2)
            continue

        results.append({"id": oid, "geometry": incremental,
                         "area_km2": area_km2, "snap_x": sx, "snap_y": sy})
        cumulative_union = incremental if cumulative_union is None else make_valid(
            cumulative_union.union(incremental)
        )
        logger.info("  → %.3f km² incremental area recorded", area_km2)

    # ── 4. Compile and export ─────────────────────────────────────────────────
    if not results:
        raise RuntimeError("No valid nested catchments produced — check DEM conditioning and outlet coordinates")

    out_gdf = gpd.GeoDataFrame(results, crs=f"EPSG:{epsg}")
    out_gdf = out_gdf.sort_values("area_km2", ascending=False).reset_index(drop=True)
    out_gdf["hierarchy_level"] = out_gdf.index + 1

    out_gdf.to_file(output_gpkg, driver="GPKG")
    logger.info("Saved %d nested catchments to %s (total %.2f km²)",
                len(out_gdf), output_gpkg, out_gdf["area_km2"].sum())
    return out_gdf

Validation Protocol

A completed nested delineation run must pass three classes of checks before the output is used in any hydrologic model.

Area conservation check

The sum of all incremental polygon areas must equal the full upstream contributing area of the most-downstream outlet within an acceptable tolerance (typically ±2 % for 30 m DEMs; ±0.5 % for 1 m LiDAR).

python

full_upstream_km2 = float(acc_arr.max() * (grid.affine.a ** 2) / 1e6)  # approx from raster
sum_incremental   = out_gdf["area_km2"].sum()
pct_error = abs(sum_incremental - full_upstream_km2) / full_upstream_km2 * 100
logger.info("Area conservation error: %.2f%%", pct_error)
assert pct_error < 2.0, f"Area error {pct_error:.2f}% exceeds 2% tolerance"

Topological overlap check

No two incremental polygons may share interior area. Violations indicate either an incorrect processing order or a failing make_valid call.

python

from itertools import combinations

for i, j in combinations(range(len(out_gdf)), 2):
    overlap = out_gdf.geometry.iloc[i].intersection(out_gdf.geometry.iloc[j])
    if not overlap.is_empty and overlap.area > 1.0:  # 1 m² threshold for floating-point noise
        logger.error("Overlap between outlets %s and %s: %.2f m²",
                     out_gdf["id"].iloc[i], out_gdf["id"].iloc[j], overlap.area)

NHD/authoritative hydrography overlay

Dissolve all incremental polygons back to the full basin and compare with an authoritative boundary from USGS or NHD. For stream threshold tuning contexts, compare the extracted stream network against NHD flowlines to confirm the delineated divides align with mapped channel heads.

Geometry validity

python

invalid = out_gdf[~out_gdf.geometry.is_valid]
if not invalid.empty:
    logger.error("%d invalid geometries: %s", len(invalid), invalid["id"].tolist())

Common Failure Modes & Optimization

Double-counted area (overlapping polygons) Root cause: outlets processed in the wrong order, or a failing make_valid silently returning the original invalid geometry. Mitigation: assert incremental.is_valid after every difference operation; log cumulative union area after each step and verify monotonic growth.

Outlet snapped to wrong tributary Root cause: snap_radius_px too large relative to channel spacing. At 30 m resolution in a braided floodplain, a radius of 7 cells may jump across a channel divide. Mitigation: reduce radius and visually inspect snapped points against the flow accumulation raster in QGIS before running the full extraction. The outlet point mapping and validation workflow provides systematic QA for this step.

Memory exhaustion on large DEMs Root cause: grid.catchment() returns a full-grid mask for every outlet, allocated in RAM. For domains exceeding ~5 000 km² at 10 m resolution, a single mask can exceed 4 GB. Mitigation: use rasterio.windows to clip the DEM to a bounding box padded around each outlet’s known upstream extent before loading into pysheds, or switch to a tile-based whitebox workflow.

Slivers and zero-area fragments Root cause: raster-to-vector conversion of a binary catchment mask produces artefact polygons at raster edges and nodata boundaries. Mitigation: apply the min_area_km2 threshold (default 0.01 km²) and make_valid after every difference. Use shapely.ops.unary_union before recording the cumulative union to collapse adjacent fragments.

Flat-area stagnation producing undefined flow direction Root cause: depression-filling alone does not assign a unique flow direction to extended flat areas. Mitigation: apply grid.resolve_flats() after fill_depressions() — this assigns gradient-based pseudo-directions away from flat edges. See DEM pit filling algorithms for conditioning strategies.

Coordinate mismatch causing empty catchment masks Root cause: outlet coordinates in a different CRS than the DEM. The inverse affine transform places the point outside the raster extent, and grid.catchment() returns an empty mask. Mitigation: assert gdf.crs.to_epsg() == epsg before snapping; reproject outlets if needed.

When to Use This vs. Alternatives

Use nested catchment delineation when you need:

Incremental drainage areas between consecutive gauges on the same stream (inter-gauge sub-basins for routing models like HEC-HMS or SWAT)
Sub-basin geometries that are topologically children of a parent basin without overlap or gap
Attribution of water quality monitoring data to contributing areas at multiple scales simultaneously

Consider basin partitioning strategies instead when:

Sub-basins are independent (no parent-child relationship required), for instance when partitioning a region into computational units for a distributed rainfall-runoff model
Sub-basin boundaries should follow Strahler order or equal-area rules rather than monitoring point locations

Consider multi-directional routing (covered in the D-Infinity routing patterns guide) when hillslope dispersion matters more than clean sub-basin boundaries — for instance in wetland or fan-delta terrain where D8 over-concentrates flow into single-cell paths.

Watershed Delineation & Catchment Synchronization — parent section covering the full delineation and synchronization pipeline
Basin Partitioning Strategies — alternative partitioning approaches including Strahler-order and equal-area decomposition
Outlet Point Mapping & Validation — systematic QA for outlet coordinates before delineation
Boundary Topology Validation — geometry repair and topological integrity checks for delineated polygons
DEM Pit Filling Algorithms — conditioning strategies that determine flow routing quality upstream of this workflow
D8 Flow Direction Implementation — the flow routing algorithm underlying nested delineation

Frequently Asked Questions

Why must outlets be processed downstream-to-upstream?

Processing upstream outlets first causes their contributing areas to be subtracted from partially-built parent basins, producing geometries that are missing sub-basins or contain holes. Sorting by descending flow accumulation value guarantees each parent basin is established before any of its children are cut from it: the downstream outlet has the highest accumulation, so it is always processed first.

What snapping radius should I use?

A search radius of 3–7 DEM cells balances correction of minor GPS offsets against the risk of relocating an outlet to a different tributary. For 30 m DEMs this is 90–210 m; for 10 m LiDAR, 30–70 m. Always validate by overlaying snapped points on the flow accumulation raster before running the full extraction.

How do I handle two outlets on the same stream reach?

Assign each outlet a unique id and sort by accumulation value. The upstream member has a lower accumulation value and is processed later; its contributing area is then differenced against the already-recorded downstream basin, yielding the correct incremental drainage area between the two gauges.

Related topics