Partitioning large watersheds into sub-basins with RichDEM

RichDEM’s Python bindings provide a C+±backed, OpenMP-parallelized route from raw elevation data to labelled sub-basin rasters without touching a proprietary GIS application. The technique on this page — conditioning the DEM, computing flow direction, thresholding the accumulation surface, and calling rd.Watersheds — sits squarely inside the Basin Partitioning Strategies workflow, which is itself a core stage of the Watershed Delineation & Catchment Synchronization pipeline. Understanding the internal mechanics of each RichDEM call helps you tune thresholds, anticipate memory limits, and validate that the resulting partitions are hydrologically sound.

Prerequisites

This page assumes you already have a hydrologically plausible DEM. If your elevation model still contains sensors noise, vegetation returns, or unprocessed sinks, start with DEM Pit Filling Algorithms before continuing here. The specific requirements for this technique are:

OpenMP thread control: set OMP_NUM_THREADS to the number of physical cores (not logical threads) before importing richdem. Hyperthreading rarely helps raster routing and can cause cache thrashing.

How rd.Watersheds partitions a basin

The diagram below shows how RichDEM’s four main calls chain into a complete sub-basin extraction. The flow direction grid is the shared dependency: it drives both the accumulation surface (used for stream extraction) and the watershed transform itself.

rd.Watersheds performs a reverse-flow traversal: starting from every cell that the stream mask marks as 1, it walks upward through the flow direction grid and stamps every upstream cell with the label of the stream cell it drains to. Cells that drain to a tie point (equally distant from two stream cells in the direction graph) receive the label of whichever stream cell is reached first in the traversal order. This is deterministic for D8 and approximately so for D∞.

Annotated code example

The function below implements the complete pipeline with structured logging. Every non-obvious call is explained inline.

import os
import logging
import numpy as np
import rasterio
import richdem as rd

# Set thread count BEFORE importing richdem submodules that initialise OpenMP
os.environ["OMP_NUM_THREADS"] = "8"

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s [%(levelname)s] %(message)s"
)
log = logging.getLogger(__name__)


def partition_subbasins(
    dem_path: str,
    out_dir: str,
    accumulation_threshold: int = 1000,
    method: str = "D8",
) -> dict:
    """
    Partition a large watershed into sub-basins using RichDEM.

    Parameters
    ----------
    dem_path : str
        Path to the input DEM GeoTIFF (projected CRS, nodata=-9999).
    out_dir : str
        Directory for output rasters (created if absent).
    accumulation_threshold : int
        Minimum upstream cells to define a stream initiation point.
        Lower values = finer sub-basins; higher values = fewer, larger ones.
    method : str
        Routing algorithm: 'D8' (single-direction) or 'Dinf' (divergent).

    Returns
    -------
    dict with paths to 'subbasins', 'flow_accum', and 'stream_network' outputs.
    """
    os.makedirs(out_dir, exist_ok=True)
    routing = "Dinf" if method.upper() in ("DINF", "D-INF", "DINFINITY") else "D8"
    log.info("Starting sub-basin partition: method=%s, threshold=%d", routing, accumulation_threshold)

    # Step 1: capture spatial reference BEFORE handing the file to richdem.
    # RichDEM's rdArray does not expose the full rasterio CRS/transform objects,
    # so we read them independently and use them when writing outputs.
    with rasterio.open(dem_path) as src:
        ref_transform = src.transform
        ref_crs = src.crs
        log.info("DEM extent: %s  CRS: %s", src.bounds, ref_crs)

    dem = rd.LoadGDAL(dem_path, no_data=-9999)
    if dem is None:
        raise RuntimeError(f"rd.LoadGDAL returned None for {dem_path!r}. "
                           "Verify GDAL driver support and nodata value.")

    # Step 2: fill depressions using the Priority-Flood algorithm (Barnes 2014).
    # in_place=False keeps the original array intact for difference-raster validation.
    log.info("Filling depressions (Priority-Flood)...")
    dem_filled = rd.FillDepressions(dem, in_place=False)

    # Step 3: compute flow accumulation.
    # rd.FlowAccumulation internally calls FlowDirection then accumulates.
    # Passing the *filled* DEM here is critical; passing the raw DEM
    # produces accumulation errors wherever sinks remain.
    log.info("Computing %s flow accumulation...", routing)
    acc = rd.FlowAccumulation(dem_filled, method=routing)
    acc_np = np.array(acc)  # copy to numpy; avoids modifying the richdem buffer

    # Step 4: threshold the accumulation surface to produce a binary stream mask.
    # Every cell with acc >= threshold is treated as a stream initiation point.
    # rd.Watersheds uses the *positions* of 1-valued cells as outlet seeds.
    log.info("Extracting streams with threshold >= %d cells...", accumulation_threshold)
    streams = (acc_np >= accumulation_threshold).astype(np.int32)
    stream_count = int(streams.sum())
    log.info("Stream cells identified: %d", stream_count)
    if stream_count == 0:
        raise ValueError("Threshold too high — no stream cells found. "
                         "Lower accumulation_threshold or check DEM extent.")

    # Step 5: compute flow direction separately for rd.Watersheds.
    # Note: FlowAccumulation above derived flow direction internally but does
    # not expose it. We must call FlowDirection explicitly to obtain the grid
    # that Watersheds requires as its first argument.
    log.info("Computing %s flow direction for watershed labelling...", routing)
    flow_dir = rd.FlowDirection(dem_filled, method=routing)

    # Step 6: label sub-basins. rd.Watersheds traverses flow_dir in reverse,
    # starting from every non-zero cell in streams, and stamps each upstream
    # cell with the label of its stream-seed cell.
    log.info("Labelling sub-basins via rd.Watersheds...")
    subbasins = rd.Watersheds(flow_dir, streams)
    subbasins_np = np.array(subbasins)
    unique_basins = len(np.unique(subbasins_np[subbasins_np > 0]))
    log.info("Sub-basins identified: %d", unique_basins)

    # Step 7: export all three output rasters with the CRS and transform
    # captured in step 1. Tiled LZW output keeps files seekable for downstream
    # rasterio.windows reads without loading the entire array again.
    outputs = {
        "subbasins": os.path.join(out_dir, "subbasins.tif"),
        "flow_accum": os.path.join(out_dir, "flow_accum.tif"),
        "stream_network": os.path.join(out_dir, "stream_network.tif"),
    }
    for key, arr in [
        ("subbasins", subbasins_np),
        ("flow_accum", acc_np),
        ("stream_network", streams),
    ]:
        _write_raster(arr, outputs[key], ref_transform, ref_crs)
        log.info("Wrote %s → %s", key, outputs[key])

    log.info("Pipeline complete.")
    return outputs


def _write_raster(
    data: np.ndarray,
    out_path: str,
    transform,
    crs,
    nodata: float = -9999.0,
) -> None:
    """Write a numpy array to a tiled, LZW-compressed GeoTIFF."""
    with rasterio.open(
        out_path, "w",
        driver="GTiff",
        height=data.shape[0],
        width=data.shape[1],
        count=1,
        dtype=data.dtype,
        crs=crs if crs is not None else "EPSG:4326",
        transform=transform,
        compress="LZW",
        tiled=True,
        blockxsize=512,
        blockysize=512,
        nodata=nodata,
    ) as dst:
        dst.write(data, 1)

Parameter reference

For a detailed discussion of threshold selection methodology, see stream threshold tuning and the worked comparison of D8 versus D∞ accuracy in comparing D8 vs D∞ for steep terrain hydrology.

Worked example and output interpretation

Running partition_subbasins("basin_10m.tif", "output/", accumulation_threshold=800, method="D8") on a 1 200 × 900 cell DEM (10 m resolution, ~108 km²) typically produces:

flow_accum.tif — values range from 1 (headwater cells) to the basin area in cells (here ≈1 080 000 for the outlet cell). A correct accumulation raster shows a clear branching pattern when rendered; flat regions or concentric rings indicate unfilled sinks.

stream_network.tif — binary mask (0/1). At threshold 800 the network should show hierarchically connected channels. Spot-check that channel heads appear on hillslopes rather than in flat valley bottoms (a sign of threshold too low) and that main channels are continuous (broken channels indicate accumulation artifacts or a threshold that is too high).

subbasins.tif — integer labels 1…N, where N equals the number of stream cells (seeds). Each label covers all cells that ultimately drain to that seed. Correct output shows compact, contiguous regions whose boundaries trace topographic divides. Fragmented patches (small isolated islands of a single label) indicate that the DEM still contains unfilled pits that interrupt the reverse-flow traversal.

When integrating sub-basin outputs into distributed models (SWAT, HEC-HMS), verify that the BASIN_ID attribute joins cleanly on the integer label, that -9999 nodata is excluded from raster statistics, and that the CRS matches the model’s project CRS — consult coordinate reference system alignment if reprojection is needed.

Gotchas and edge cases

• rd.FlowAccumulation does not return the flow direction grid. The internal direction computation is not exposed by the Python binding. You must call rd.FlowDirection separately to obtain the grid required by rd.Watersheds. Omitting this step and passing the accumulation array as the first argument to rd.Watersheds will raise a silent type mismatch or produce nonsensical labels.

• Mismatched nodata between richdem and rasterio. If the source GeoTIFF declares nodata as None or a value other than -9999, rd.LoadGDAL may treat border pixels as valid terrain. Always set nodata explicitly in both tools, or reproject/rewrite the DEM with gdal_translate -a_nodata -9999 before loading.

• Memory footprint is 2–3× the raw DEM size. RichDEM holds the input rdArray, the filled copy, the flow direction grid, and the accumulation surface in RAM simultaneously. For DEMs exceeding 8–12 GB uncompressed, pre-tile with gdal_translate -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 and process each tile independently with a boundary buffer of ≥500 m to capture cross-tile upstream contributions. Mosaic sub-basin rasters afterward, adding a per-tile integer offset to label IDs so that labels remain globally unique.

• D∞ fractional routing and rd.Watersheds compatibility. rd.Watersheds was designed for single-direction flow grids. When method="Dinf", each cell’s flow is split between two neighbors, so the reverse traversal may assign a cell to whichever neighbor is processed first rather than the dominant-flow neighbor. For strict partition accuracy on D∞ outputs, consider using the flow direction only to determine the primary (larger-proportion) neighbor, effectively falling back to D8 semantics at the watershed step.

• Zero-area sub-basins at DEM edges. Cells on the raster boundary that receive no upstream flow get a label from the nearest stream seed. If the DEM was clipped tightly to a watershed boundary, edge cells may receive incorrect labels from outside-basin stream seeds. Apply a one-cell inward buffer mask on the sub-basin output and set edge labels to nodata.

• in_place=True prevents pre/post comparison. If you call rd.FillDepressions(dem, in_place=True), the original elevation values are overwritten and you cannot compute a difference raster to validate how many cells were modified. Keep in_place=False during development and switch to True only when RAM is genuinely constrained and the fill step has been independently validated.

Related topics

• Basin Partitioning Strategies — parent page covering the full partitioning workflow, topology validation, and production automation patterns
• Watershed Delineation & Catchment Synchronization — grandparent covering the end-to-end pipeline from DEM conditioning to catchment synchronization
• DEM Pit Filling Algorithms — prerequisite conditioning step; use before any RichDEM routing call
• Stream Threshold Tuning — detailed methodology for selecting accumulation_threshold values
• D-Infinity Routing Patterns — when to switch from D8 to D∞ and how divergent routing affects sub-basin geometry