⚙️ How-To Hydrology & watersheds Terrain & elevation Batch processing Hydrology & Watersheds · step 10 of 20 Updated

How to Split a Catchment into Sub-basins at Many Outlets

Problem statement

A single watershed answers one question: what drains to this point? Monitoring networks, water-quality programmes and hydrological models ask the same question at dozens of points at once, and they usually want the answer split two ways. Nested catchments give everything upstream of each outlet, so they overlap. Incremental sub-basins give only the land between an outlet and the next outlets upstream, so they tile the catchment without overlap and add up to the whole.

Done for the 87 USGS monitoring locations with a published drainage area of at least 1 mi² in the Esopus Creek basin, New York, on a 10 m breached 3DEP DEM:

  • Snapping each gauge to the nearest channel cell with an accumulated area within a factor of two of its published area placed 86 of 87 gauges; the median area error was 1.1%, 80 were within 5%, and one gauge could not be matched at all.
  • Tracing 87 nested catchments one by one took 5.3 s; labelling every cell with its first downstream outlet in one pure-Python pass over 8.3 million cells took 8.1 s, and WhiteboxTools' watershed did the same in 0.30 s, agreeing on 100% of cells.
  • The 87 incremental sub-basins summed to 492.02 km², exactly the nested catchment at the basin outlet, with a median of 2.47 km² and one of 0.00 km² where two gauges sat almost on top of each other.
  • The outlets formed a tree 20 levels deep: the Coldbrook gauge at the bottom, 20 headwater gauges with nothing upstream.

Quick answer 

import whitebox

wbt = whitebox.WhiteboxTools()
wbt.d8_pointer("/data/esopus_breached.tif", "/data/esopus_d8.tif")
wbt.watershed("/data/esopus_d8.tif", "/data/outlets_snapped.shp", "/data/subbasins.tif")

Give watershed all the snapped outlets at once and it labels each cell with the first outlet downstream of it: incremental sub-basins in one call. Snap the outlets first, and check the result adds up.

Two panels contrasting nested catchments, which overlap and each include everything upstream, with incremental sub-basins, which tile the basin and sum to the outlet area.
Nested areas answer "what drains here"; incremental areas answer "what drains here and nowhere earlier".

Step-by-step solution

1. Condition the DEM once

Every outlet uses the same flow directions, so condition and route the whole DEM once. This guide starts from the breached 10 m DEM made in conditioning a DEM.

2. Collect outlets with their published areas

Gauge coordinates are rarely on the modelled channel. The published drainage area is the best evidence of where a gauge belongs: of the 500 USGS monitoring locations in the bounding box, 182 had a drainage area, and 87 with at least 1 mi² lay within 200 m of the basin boundary.

3. Snap each outlet by area, not by maximum accumulation

Snapping to the highest accumulation nearby works for one gauge on a main river. With many gauges it moves tributary gauges onto the main stem beside them. Instead, look within 150 m for cells whose accumulated area is between half and twice the published area, and take the nearest. That matched 86 of 87 gauges; the remaining one fell back to the maximum and was 100% out (Example 1).

4. Drop outlets that snap to the same cell

Two gauges on the same cell would create an empty sub-basin and an ambiguous label. None of the 87 collided here, but in dense networks it is common; keep one per cell.

5. Label every cell with its first downstream outlet

Walking downstream from each cell until an outlet is reached — and remembering the answer for every cell on the path — labels the whole grid in one pass (Example 2). The alternative, one catchment trace per outlet, gives nested areas: 5.3 s for 87 outlets, or 0.06 s each, and the traces overlap.

6. Build the outlet tree

The cell just downstream of each outlet carries the label of the next outlet down. That gives a drains_to table: a tree with Coldbrook at the root, 20 headwater outlets with nothing upstream, and two outlets receiving three or more direct upstream outlets. The longest chain was 20 outlets deep, down the Esopus main stem.

7. Check the areas add up

The incremental areas must sum to the nested area at the lowest outlet: 492.02 km² for both. A nested area can also be rebuilt by summing an outlet's incremental area and those of everything upstream of it in the tree. If the sums differ, some cells drained out of the grid or to a pit before reaching an outlet.

8. Polygonise and store both views

Polygonise the label raster once, dissolve by label, and store the incremental polygons with nested area, incremental area and the downstream outlet (Example 3). A quarter of the sub-basins came out multipart because diagonal D8 connections join cells that share only a corner.

Bar chart of the time to delineate 87 outlets: a per-outlet catchment loop, a one-pass Python labelling, WhiteboxTools watershed, WhiteboxTools unnest_basins and polygonising the labels.
One labelling pass beats 87 traces, and compiled code beats both.

Code examples

Example 1 — snap many gauges by their published area 

import numpy as np

if not hasattr(np, "in1d"):                                     # pysheds 0.5 on NumPy 2.4+
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
import geopandas as gpd
import pandas as pd
from pysheds.grid import Grid

grid = Grid.from_raster("esopus_breached.tif")
fdir = grid.flowdir(grid.read_raster("esopus_breached.tif"))
F = np.asarray(fdir)
acc = np.asarray(grid.accumulation(fdir))
cell_m = abs(grid.affine.a)
cell_km2 = cell_m * abs(grid.affine.e) / 1e6

basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(grid.crs.srs)
gauges = gpd.read_file("usgs_gauges_esopus.geojson").to_crs(grid.crs.srs)
gauges = gauges[gauges.drainage_area.ge(1.0) & gauges.within(basin.geometry.iloc[0].buffer(200))]


def snap_by_area(point, area_km2, radius_m=150):
    col, row = ~grid.affine * (point.x, point.y)
    row, col = int(row), int(col)
    k = int(radius_m / cell_m)
    window = acc[row - k:row + k + 1, col - k:col + k + 1] * cell_km2
    dy, dx = np.mgrid[-k:k + 1, -k:k + 1]
    distance = np.hypot(dy, dx) * cell_m
    plausible = (window > area_km2 / 2) & (window < area_km2 * 2) & (distance <= radius_m)
    if plausible.any():
        j = np.unravel_index(np.where(plausible, distance, np.inf).argmin(), window.shape)
    else:
        j = np.unravel_index(np.where(distance <= radius_m, window, -1).argmax(), window.shape)
    return row - k + j[0], col - k + j[1], bool(plausible.any())


rows = []
for site, name, point, mi2 in zip(gauges.monitoring_location_number, gauges.monitoring_location_name,
                                   gauges.geometry, gauges.drainage_area):
    row, col, plausible = snap_by_area(point, mi2 * 2.589988)
    rows.append({"site": site, "name": name, "published_km2": mi2 * 2.589988, "row": row, "col": col,
                 "nested_km2": acc[row, col] * cell_km2, "matched": plausible})
outlets = pd.DataFrame(rows).drop_duplicates(["row", "col"]).reset_index(drop=True)
error = (outlets.nested_km2 / outlets.published_km2 - 1).abs()
print(f"{len(rows)} gauges, {len(outlets)} distinct outlets, {outlets.matched.sum()} matched by area; "
      f"median error {error.median():.1%}, within 5% {(error <= 0.05).sum()}, worst {error.max():.0%}")
87 gauges, 87 distinct outlets, 86 matched by area; median error 1.1%, within 5% 80, worst 100%

Example 2 — label incremental sub-basins in one pass 

import time

# pysheds' default D8 codes: N 64, NE 128, E 1, SE 2, S 4, SW 8, W 16, NW 32
STEP = {64: (-1, 0), 128: (-1, 1), 1: (0, 1), 2: (1, 1), 4: (1, 0), 8: (1, -1), 16: (0, -1), 32: (-1, -1)}


def label_subbasins(F, outlets):
    """Label each cell with the first outlet downstream of it (0 if it reaches none)."""
    nrow, ncol = F.shape
    labels = np.full(F.size, -1, dtype=np.int32)
    for k, (r, c) in enumerate(outlets, start=1):
        labels[r * ncol + c] = k
    codes = F.ravel()
    for start in range(F.size):
        path, cur = [], start
        while labels[cur] == -1:
            path.append(cur)
            dr, dc = STEP.get(int(codes[cur]), (0, 0))
            r, c = divmod(cur, ncol)
            if (dr, dc) == (0, 0) or not (0 <= r + dr < nrow and 0 <= c + dc < ncol):
                labels[path] = 0
                break
            cur = (r + dr) * ncol + c + dc
        else:
            labels[path] = labels[cur]
    return labels.reshape(F.shape)


start = time.perf_counter()
labels = label_subbasins(F, list(zip(outlets.row, outlets.col)))
outlets["incremental_km2"] = np.bincount(labels.ravel(), minlength=len(outlets) + 1)[1:] * cell_km2
print(f"labelled {F.size:,} cells in {time.perf_counter() - start:.1f} s")


def next_outlet(r, c):
    dr, dc = STEP.get(int(F[r, c]), (0, 0))
    inside = (dr, dc) != (0, 0) and 0 <= r + dr < F.shape[0] and 0 <= c + dc < F.shape[1]
    return int(labels[r + dr, c + dc]) if inside else 0


outlets["drains_to"] = [next_outlet(r, c) for r, c in zip(outlets.row, outlets.col)]
depth = []
for k in range(1, len(outlets) + 1):
    d = 0
    while outlets.drains_to.iloc[k - 1] > 0:
        k, d = outlets.drains_to.iloc[k - 1], d + 1
    depth.append(d)
outlets["depth"] = depth
root = outlets.nested_km2.idxmax()
print(f"incremental total {outlets.incremental_km2.sum():.2f} km2 vs nested at {outlets.site[root]} "
      f"{outlets.nested_km2[root]:.2f} km2; median {outlets.incremental_km2.median():.2f} km2; "
      f"headwaters {int((~outlets.index.isin(outlets.drains_to - 1)).sum())}; deepest {max(depth)}")
labelled 8,348,820 cells in 8.1 s
incremental total 492.02 km2 vs nested at 01362500 492.02 km2; median 2.47 km2; headwaters 20; deepest 20

Example 3 — cross-check with WhiteboxTools and save polygons 

import os

import rasterio
import whitebox
from rasterio import features
from shapely.geometry import Point, shape

here = os.path.abspath(".")
wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(here)
centres = [Point(*(grid.affine * (c + 0.5, r + 0.5))) for r, c in zip(outlets.row, outlets.col)]
gpd.GeoDataFrame({"fid": range(1, len(outlets) + 1)}, geometry=centres, crs=grid.crs.srs).to_file("outlets_snapped.shp")
wbt.d8_pointer("esopus_breached.tif", "esopus_d8.tif")
start = time.perf_counter()
wbt.watershed("esopus_d8.tif", "outlets_snapped.shp", "subbasins.tif")
seconds = time.perf_counter() - start
with rasterio.open("subbasins.tif") as src:
    wb = src.read(1)
either = (labels > 0) | (wb > 0)
print(f"whitebox watershed {seconds:.2f} s; cells agreeing with the Python labels {np.mean(labels[either] == wb[either]):.2%}")

pieces = [{"fid": int(v), "geometry": shape(g)}
          for g, v in features.shapes(labels, mask=labels > 0, transform=grid.affine)]
subbasins = gpd.GeoDataFrame(pieces, crs=grid.crs.srs).dissolve("fid").reset_index()
subbasins = subbasins.merge(outlets.assign(fid=outlets.index + 1)[["fid", "site", "name", "published_km2", "nested_km2",
                                                                    "incremental_km2", "drains_to", "depth"]], on="fid")
subbasins.to_file("esopus_subbasins.gpkg")
print(f"{len(subbasins)} sub-basins from {len(pieces)} pieces, {subbasins.area.sum() / 1e6:.2f} km2, "
      f"{(subbasins.geom_type == 'MultiPolygon').sum()} multipart")
print(subbasins.nlargest(5, "nested_km2")[["site", "published_km2", "nested_km2", "incremental_km2", "depth"]].round(2).to_string(index=False))
whitebox watershed 0.30 s; cells agreeing with the Python labels 100.00%
87 sub-basins from 149 pieces, 492.02 km2, 25 multipart
    site  published_km2  nested_km2  incremental_km2  depth
01362500         497.28      492.02            11.87      0
01362430         372.96      372.70             8.19      1
01362420         365.19      364.51             7.30      2
01362405         357.42      357.20             3.89      3
01362250         215.75      216.10            11.09      4

Explanation

Nested and incremental are two views of one tree

Each outlet's nested catchment is its own incremental sub-basin plus the incremental sub-basins of every outlet upstream. Store the incremental polygons and the drains_to tree, and any nested catchment is a dissolve away; store only nested polygons and you have 87 overlapping shapes that are hard to combine.

Why one pass is faster than one trace per outlet

A catchment trace starts at an outlet and searches upstream, visiting every cell in the nested catchment. Downstream outlets repeat the work of upstream ones: the 87 traces visited the Coldbrook catchment's 4.9 million cells once for Coldbrook and again, in part, for every gauge above it. A downstream walk with memoisation visits each cell once, whatever the number of outlets. Compiled, that is the 0.30 s of watershed; in pure Python it is 8.1 s, independent of the number of outlets.

Why area-based snapping matters more with many outlets

With one outlet on a large river, the highest accumulation nearby is almost always correct. Many gauges sit on small tributaries a few tens of metres from a larger channel, and there the maximum is the wrong river. The published drainage area picks the channel of the right size; see pour points explained.

What unnest_basins adds

WhiteboxTools' unnest_basins writes nested catchments without overlap problems by putting each nesting level in its own raster: 21 rasters for a 20-level tree, in 4.28 s. It is useful when a model needs nested catchment rasters; for polygons, the incremental labels and the tree are more compact.

Table of the five largest Esopus Creek gauges with published, nested and incremental areas and depth in the outlet tree.
The main stem gauges have large nested areas and small incremental ones: most of their water arrives through the gauge above.

Edge cases or notes

  • Gauges very close together give incremental sub-basins of a few cells; Stony Clove Creek at Phoenicia had 0.00 km² because a gauge just upstream took its whole catchment.
  • Outlets off the DEM or in NoData get no label; check that each outlet cell has a label equal to its own index.
  • Cells that drain to the grid edge or a pit get 0; a large 0 area inside the basin means poor conditioning.
  • Published areas that include non-contributing land will not match any cell; widen the plausibility factor for those gauges.
  • Diagonal D8 links make multipart polygons; buffer by a tiny distance and dissolve if a model needs single parts.
  • Lakes and reservoirs between gauges distort incremental areas when the DEM routes flow across the water surface incorrectly.
  • Large regions make the pure-Python pass slow; use the WhiteboxTools tool or a compiled loop.

FAQ

What is the difference between nested and incremental sub-basins?

A nested catchment is everything upstream of an outlet, so catchments overlap. An incremental sub-basin is only the area between an outlet and the next outlets upstream; the 87 Esopus sub-basins summed to exactly the 492.02 km² catchment at the lowest gauge.

How do I delineate watersheds for many points at once?

Pass all snapped outlets to one labelling step, such as WhiteboxTools' watershed tool with a multi-point shapefile. It labelled 87 sub-basins in 0.30 s.

Why do my sub-basins not add up to the whole catchment?

Some cells drain to the grid edge, to a pit or to an outlet that is not in your set. Check the area labelled 0 inside the basin and the conditioning of the DEM.

How should I snap many gauges to a DEM stream network?

Use the published drainage area: choose the nearest cell whose accumulated area is within a factor of two of it. That matched 86 of 87 gauges with a median area error of 1.1%.

How do I get nested catchments from incremental sub-basins?

Follow the drains_to tree upstream from an outlet and dissolve its incremental sub-basin with all those above it.

Why are some sub-basin polygons multipart?

D8 routes flow diagonally between cells that touch only at a corner, and polygonising treats those as separate parts. 25 of 87 sub-basins were multipart in the measured run.