Implementing D8 Flow Routing with RichDEM Python Bindings

rd.FlowDirection() and rd.FlowAccumulation() with method='D8' are the two RichDEM calls that take you from a conditioned DEM to a routable drainage network without leaving Python. This page details exactly what those calls accept, what they return, and where they fail — as a focused implementation reference under the D8 Flow Direction Implementation workflow, which sits within the broader Flow Routing Algorithms & Stream Network Extraction coverage on this site.

The richdem package wraps highly optimized C++ terrain-analysis routines behind a thin Python API. The result is a library that executes watershed-scale D8 routing in seconds without requiring external GIS applications or manual GDAL driver configuration.

Prerequisites

This page assumes you have already addressed upstream data-quality concerns:

• Your DEM is free of projection inconsistencies — see coordinate reference system alignment if your source rasters arrive in mixed CRS.
• You understand why DEM pit filling algorithms must precede any flow-routing step. RichDEM applies its own FillDepressions() (described below), but understanding what it does prevents surprises when sink counts are unexpectedly high.
• Your DEM is a single-band float32 or float64 GeoTIFF in a projected coordinate system with linear units in metres.

Environment

RichDEM requires compiled C++ extensions and GDAL I/O, which introduces strict version dependencies:

• Python: 3.8–3.11. Python 3.12 may require a conda-forge rebuild due to NumPy ABI changes.
• GDAL: 3.4–3.7. Mismatched GDAL/PROJ versions produce ImportError: libgdal.so or, worse, silent coordinate corruption.
• OS: Linux and macOS have native wheel support. On Windows, install via conda-forge to avoid MSVC runtime DLL failures.
• RAM: A 10,000 × 10,000 float32 DEM needs roughly 400 MB to load; flow accumulation nearly doubles the footprint. Size your environment accordingly before running regional analyses.

conda create -n hydro-d8 python=3.10 richdem rasterio numpy -c conda-forge
conda activate hydro-d8

Core Technique: the RichDEM D8 Binding Calls

The three function calls that do the actual work are rd.LoadGDAL(), rd.FlowDirection(), and rd.FlowAccumulation(). Understanding exactly what each call does — and what it does not do — prevents the most common integration mistakes.

Pipeline: DEM to routable network

D8 direction encoding: how values 1–8 map to neighbours

Each output cell in the direction grid receives an integer that identifies which of the eight surrounding neighbours is the steepest descent target. RichDEM uses clockwise encoding starting from east.

Cells that drain off the DEM edge, or that sit on a flat area RichDEM cannot resolve, receive value 0. An output dominated by a single numeric value — for example 90 % of cells encoded as 1 (east) — is a strong signal that the DEM’s projection uses geographic degrees rather than projected metres, causing near-zero east–west gradients relative to diagonal slopes.

rd.LoadGDAL() — what it returns

rd.LoadGDAL(path) returns an rdarray, which is a NumPy ndarray subclass that carries extra attributes: geotransform (a six-element GDAL affine tuple), projection (a WKT string), and no_data. These attributes travel with the array through the fill and direction calculations so that rd.SaveGDAL() can reconstruct a fully georeferenced output without any manual GDAL driver configuration.

rd.FlowDirection() — method='D8'

flow_dir = rd.FlowDirection(dem_filled, method='D8')

This assigns each cell to exactly one of its eight immediate neighbours based on the steepest downhill gradient. Flat areas are handled by a secondary slope-routing step that RichDEM applies automatically after the primary gradient pass.

rd.FlowAccumulation() — passing a conditioned DEM vs. a direction grid

There are two valid call signatures:

# Option A: let RichDEM derive direction internally (convenient for one-shot scripts)
flow_acc = rd.FlowAccumulation(dem_filled, method='D8')

# Option B: supply a precomputed direction grid (preferred when you need both outputs)
flow_acc = rd.FlowAccumulation(flow_dir)   # no method arg needed

Option B is preferable when you need to inspect or export flow_dir as a separate product, because it avoids computing direction twice. If you supply a direction grid and a method string simultaneously, the method string is silently ignored — worth noting in log output.

Annotated Code Example

The script below is a complete, production-ready D8 routing function. Every non-obvious line carries an inline comment.

import logging
import richdem as rd
import numpy as np
from pathlib import Path

logger = logging.getLogger(__name__)


def run_d8_routing(
    dem_path: str,
    output_dir: str,
    fill_in_place: bool = False,
) -> tuple[rd.rdarray, rd.rdarray]:
    """
    Compute D8 flow direction and accumulation using RichDEM.

    Parameters
    ----------
    dem_path     : path to a single-band GeoTIFF DEM (float32 or float64)
    output_dir   : directory for output rasters (created if absent)
    fill_in_place: if True, FillDepressions modifies the loaded array;
                   set False to keep the raw DEM unmodified for diffing

    Returns
    -------
    (flow_dir, flow_acc) as rdarray objects with geotransform/projection intact
    """
    out = Path(output_dir)
    out.mkdir(parents=True, exist_ok=True)
    logger.info("Starting D8 routing pipeline for %s", dem_path)

    # --- 1. Load ---
    dem = rd.LoadGDAL(dem_path)
    # rdarray stores .geotransform, .projection, .no_data automatically;
    # log them now so downstream CRS issues are immediately traceable
    logger.info(
        "Loaded DEM: shape=%s  dtype=%s  nodata=%s",
        dem.shape, dem.dtype, dem.no_data
    )

    # --- 2. Validate ---
    if dem.ndim != 2:
        raise ValueError(
            f"RichDEM expects a 2-D array; got shape {dem.shape}. "
            "Split multi-band DEMs before loading."
        )
    nan_count = int(np.sum(np.isnan(dem)))
    if nan_count > 0:
        raise ValueError(
            f"DEM contains {nan_count} NaN cells. "
            "Fill or interpolate them before routing — the priority-flood "
            "algorithm stalls on NaN boundaries and produces zero-accumulation rings."
        )
    # Verify the no_data sentinel is set; a missing sentinel routes edge pixels as terrain
    if dem.no_data is None:
        logger.warning(
            "no_data sentinel is None. If edge pixels use a fill value (e.g. -9999), "
            "set dem.no_data = <value> before FillDepressions."
        )
    logger.info("Validation passed (no NaNs, 2-D shape confirmed, nodata=%s).", dem.no_data)

    # --- 3. Fill depressions ---
    # priority-flood raises each sink cell to the lowest surrounding spill-point elevation,
    # ensuring every interior pixel has a valid downstream path to a DEM boundary
    logger.info("Filling depressions (in_place=%s)...", fill_in_place)
    dem_filled = rd.FillDepressions(dem, in_place=fill_in_place)
    # With in_place=False RichDEM allocates a full copy — ~400 MB extra for a 10k×10k float32

    # --- 4. Compute D8 flow direction ---
    # Output: integers 1-8 (clockwise from east) or 0 for outlets / unresolvable flats
    logger.info("Computing D8 flow direction...")
    flow_dir = rd.FlowDirection(dem_filled, method='D8')
    unique_vals = sorted(np.unique(np.asarray(flow_dir)).tolist())
    logger.info("Direction grid unique values: %s", unique_vals)
    zero_fraction = float(np.mean(np.asarray(flow_dir) == 0))
    if zero_fraction > 0.05:
        logger.warning(
            "%.1f%% of cells have direction=0 (outlet/flat). "
            "If this fraction is unexpectedly high, check depression filling and CRS.",
            zero_fraction * 100
        )

    # --- 5. Compute flow accumulation from precomputed direction grid ---
    # Passing flow_dir (not dem_filled + method) avoids recomputing direction internally
    logger.info("Computing flow accumulation from direction grid...")
    flow_acc = rd.FlowAccumulation(flow_dir)
    max_acc = int(np.max(np.asarray(flow_acc)))
    logger.info("Accumulation max = %d cells (primary outlet or main channel)", max_acc)

    # --- 6. Export --- rd.SaveGDAL copies geotransform + projection from the rdarray
    dir_out = out / "d8_flow_direction.tif"
    acc_out = out / "d8_flow_accumulation.tif"
    logger.info("Saving flow direction  → %s", dir_out)
    rd.SaveGDAL(str(dir_out), flow_dir)
    logger.info("Saving flow accumulation → %s", acc_out)
    rd.SaveGDAL(str(acc_out), flow_acc)

    logger.info("D8 routing complete.")
    return flow_dir, flow_acc


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    run_d8_routing("input_dem.tif", "./hydro_outputs")

Parameter Reference Table

Worked Example: Reading and Interpreting Outputs

Direction grid

After calling rd.FlowDirection(), log the unique values to check routing completeness:

import numpy as np

dir_array = np.asarray(flow_dir)
unique_vals = np.unique(dir_array)
logger.info("Direction encoding values present: %s", unique_vals)
# Expected: a subset of {0, 1, 2, 3, 4, 5, 6, 7, 8}
# 0 = outlet or unresolvable flat
# If value 0 covers more than ~5% of cells, depressions were not fully resolved

The encoding uses clockwise integers starting from east (see diagram above). An output dominated by a single direction value is a strong signal of CRS problems — typically a DEM left in geographic degrees so that diagonal gradients appear artificially steeper than cardinal ones.

Accumulation raster

The flow accumulation raster stores cell counts. To derive drainage area in square kilometres:

# Index 1 of geotransform is x-pixel size in projection units (metres for a projected CRS)
cell_size_m = abs(dem.geotransform[1])
area_km2 = np.asarray(flow_acc) * (cell_size_m ** 2) / 1e6
logger.info("Max drainage area at outlet: %.2f km²", float(np.max(area_km2)))

Apply a stream threshold tuning cutoff to isolate channel cells from hillslopes:

# 1 km² contributing area — adjust based on terrain type and stream permanence
threshold_cells = int(1e6 / cell_size_m ** 2)
stream_mask = np.asarray(flow_acc) >= threshold_cells
logger.info(
    "Stream network: %d cells above %.0f-cell threshold",
    int(stream_mask.sum()), float(threshold_cells)
)

For ephemeral streams and semi-arid basins the appropriate threshold is often much lower — see tuning flow accumulation thresholds for ephemeral streams for threshold-selection strategies that account for drainage density and land cover.

Validation checks

Before using the outputs in downstream delineation, run these three sanity checks:

1. Direction coverage: unique_vals should contain at least five distinct integers from 1–8. A result containing only values 1 and 5 (east and west) suggests the DEM’s row/column orientation is swapped.
2. Zero-fraction threshold: If more than 5 % of direction cells are 0, the depression-fill did not reach all sinks. Re-examine the no_data sentinel and rerun FillDepressions().
3. Accumulation maximum vs. basin area: max_acc * cell_area_km2 should approximate the known basin area in km². A value ten times larger than expected indicates the DEM edges are not behaving as outlets — set a correct no_data value or clip the DEM to the true watershed boundary before routing.

Gotchas & Edge Cases

• NaN cells stall the priority-flood. rd.FillDepressions() uses a heap-based flood-fill that stops at NaN boundaries. A single NaN island inside a basin splits the routing graph and produces a ring of zero-accumulation cells around it. Interpolate or mask NoData regions before calling FillDepressions().

• Geographic CRS produces nonsense gradients. If the DEM’s projection is in degrees (EPSG:4326 or similar), RichDEM computes rise/run using degree units, making diagonal gradients artificially steep relative to cardinal ones. Always reproject to a projected CRS before routing. If your source data is in geographic degrees, see coordinate reference system alignment for a rasterio reprojection workflow.

• Passing dem_filled to FlowAccumulation() with method='D8' is not the same as passing flow_dir. When you supply a DEM and a method string, RichDEM recomputes flow direction internally. When you supply a direction grid with no method argument, it uses the precomputed grid. Supplying both a direction grid and a method string causes the method string to be silently ignored.

• no_data sentinel mismatch causes edge routing. If your GeoTIFF uses -9999 as NoData but rdarray.no_data reads as None (possible when the GDAL metadata band does not set NODATA), edge pixels will be routed as terrain rather than masked. Always check dem.no_data immediately after LoadGDAL() and set it explicitly — dem.no_data = -9999 — before filling.

• Memory doubles during FillDepressions(in_place=False). With in_place=False, RichDEM allocates a full copy of the DEM array before modifying it. For a 20,000 × 20,000 float32 DEM (1.6 GB) that copy costs another 1.6 GB. If RAM is constrained, use in_place=True and save the raw DEM to disk before filling so you can still compute a difference raster later.

Frequently Asked Questions

Does rd.FlowAccumulation require a precomputed flow direction grid?

No. When you call rd.FlowAccumulation(dem_filled, method='D8'), RichDEM derives flow direction internally from the conditioned DEM. If you need to inspect the direction grid separately, call rd.FlowDirection(dem_filled, method='D8') first and pass the result to rd.FlowAccumulation(flow_dir) without a method argument. The two-call pattern is recommended for any workflow that exports both outputs.

What do the direction encoding values 1 through 8 represent?

RichDEM encodes the eight cardinal and diagonal neighbours as integers 1 through 8 in clockwise order starting from east: 1=E, 2=NE, 3=N, 4=NW, 5=W, 6=SW, 7=S, 8=SE. Cells with no valid downstream neighbour — outlets and unresolvable flats — are encoded as 0. The diagram in the “Core Technique” section above maps these integers visually to the 3×3 neighbourhood.

How do I convert cell-count accumulation to drainage area in square kilometres?

Multiply the flow accumulation raster by the cell area: area_km2 = flow_acc * (cell_size_m ** 2) / 1e6. Retrieve cell_size_m from the geotransform stored in dem.geotransform[1] (the x-pixel size in the DEM’s projection units). For a 10-metre DEM this gives area_km2 = flow_acc * 100 / 1e6, so a threshold of 10,000 accumulation cells corresponds to 1 km² of contributing area.

Related Topics

• D8 Flow Direction Implementation — parent page covering the algorithm mechanics, direction encoding conventions, and validation protocol for D8 in detail
• Flow Routing Algorithms & Stream Network Extraction — overview comparing D8, D-Infinity, and multiple-flow-direction methods with trade-off guidance
• D-Infinity Routing Patterns — when to move beyond single-direction routing on divergent hillslopes and how RichDEM’s Tarboton method differs
• Stream Threshold Tuning — how to select the contributing-area cutoff that separates channels from hillslopes in your accumulation raster
• DEM Pit Filling Algorithms — the broader theory behind depression removal, including comparisons of priority-flood vs. iterative filling approaches