D-Infinity Routing Patterns for Automated Watershed Modeling in Python

D-Infinity routing solves a fundamental limitation of deterministic single-direction algorithms: the tendency to concentrate all upslope area into one downstream cell and artificially channelize flow across hillslopes. By evaluating the steepest descent across eight triangular facets formed by each grid cell and its neighbors, the algorithm selects a continuous flow angle and distributes flow proportionally to the two adjacent downslope cells whose combined angle brackets it. This makes D-Infinity the method of choice for divergent hillslope hydrology, ridge dispersion, and high-resolution LiDAR-derived terrain analysis. This page is part of the Flow Routing Algorithms & Stream Network Extraction domain. For a direct runtime comparison of the two most common single-path approaches, see D8 Flow Direction Implementation; for methods that split flow across more than two receivers, see Multiple Flow Direction Methods.

Prerequisites & Environment Setup

Python stack

Install via conda-forge to guarantee binary compatibility with GDAL:

conda create -n hydro-dinf python=3.11 richdem rasterio numpy geopandas shapely scikit-image -c conda-forge
conda activate hydro-dinf

Input data specifications

• CRS: A metric projected coordinate system (UTM, State Plane, or national grid). Geographic coordinates (degrees) distort slope calculations and produce invalid flow angles.
• Vertical units: Must match horizontal units (both metres or both feet). Mismatched datums cause systematic slope errors.
• Data type: float32 or float64 for elevation values; nodata must be explicitly defined in the file metadata.
• Preprocessing: Hydrological conditioning (depression filling, flat-area resolution) must happen before D-Infinity direction computation. See DEM Pit Filling Algorithms for detail on sink-filling strategies.
• Memory budget: D-Infinity accumulation stores float64 arrays. Plan for 3–4× the DEM file size in RAM during peak computation. Datasets exceeding available RAM require tiled processing (see Production-Ready Code below).

Algorithm Mechanics

Triangular facet evaluation

Around each grid cell, D-Infinity constructs eight triangular facets. Each facet is defined by the center cell and two adjacent neighbors: one cardinal (N, E, S, W) and one diagonal (NE, SE, SW, NW). For each facet the algorithm computes the slope of the plane through the three elevation values. The plane’s steepest descent direction is expressed as a continuous angle α (radians, 0–2π measured from east, counter-clockwise).

The algorithm selects the facet with the maximum downward slope. If the steepest descent falls within the facet’s angular span, α is used directly. If it falls outside (common in convex terrain), α is clamped to the facet’s nearest bounding edge, which results in flow directed entirely to one neighbor — functionally equivalent to D8 for that cell.

Proportional partitioning

Once α is determined, flow is split between the two receivers that define the selected facet:

p  = 1 − (α − edge_angle) / facet_width    → fraction to receiver R1
1−p =    (α − edge_angle) / facet_width    → fraction to receiver R2

Where edge_angle is the angle of the facet’s first bounding edge and facet_width is the angular width of the facet (π/4 for cardinal-diagonal pairs). This proportional assignment propagates upslope contributing area as floating-point values — every cell accumulates a weighted sum of all upslope contributions, not a simple count.

Direction encoding

richdem stores D-Infinity direction as a float32 raster with values in [0, 2π]. A value of -1 (or the nodata value) marks cells that cannot drain: flat areas without resolved gradients, and true pits that remain after conditioning. The no-data flag for direction rasters must be set before passing to accumulation.

Parameter and encoding table

Step-by-Step Workflow

Step 1 — Load and validate the DEM

import logging
import rasterio
import richdem as rd

logger = logging.getLogger(__name__)

def load_and_validate_dem(path: str) -> rd.rdarray:
    """Load a DEM from *path*, assert projected CRS and defined nodata."""
    logger.info("Loading DEM from %s", path)
    with rasterio.open(path) as src:
        if not src.crs.is_projected:
            raise ValueError(
                f"DEM at {path} uses a geographic CRS ({src.crs}). "
                "Re-project to a metric CRS before routing."
            )
        if src.nodata is None:
            raise ValueError(
                f"DEM at {path} has no nodata value defined. "
                "Set nodata in the file metadata before loading."
            )
        nodata = src.nodata
        logger.info("CRS: %s | nodata: %s | shape: %s", src.crs, nodata, src.shape)

    dem = rd.LoadGDAL(path)
    dem.no_data = nodata
    return dem

Step 2 — Hydrological conditioning

D-Infinity cannot produce valid directions across unresolved depressions. Two conditioning steps are required.

def condition_dem(dem: rd.rdarray) -> rd.rdarray:
    """Fill sinks then enforce gradient across flat areas."""
    logger.info("Filling depressions with priority-flood algorithm")
    filled = rd.FillDepressions(dem, in_place=False)

    logger.info("Resolving flat areas to enforce drainage gradient")
    rd.ResolveFlats(filled, in_place=True)
    return filled

FillDepressions uses a priority-flood approach that minimally raises depression cells to the pour-point elevation, preserving surrounding terrain. ResolveFlats then imposes a tiny manufactured gradient across zero-gradient zones so the D-Infinity facet evaluator always finds a downslope direction. If you are working with coastal DEMs or tidal flats, mask marine areas before this step to prevent the ocean from acting as a giant sink that absorbs all coastal drainage.

Step 3 — Compute D-Infinity flow direction

def compute_dinfinity_direction(filled_dem: rd.rdarray) -> rd.rdarray:
    """Return per-cell flow angle raster (float32, 0–2π radians)."""
    logger.info("Computing D-Infinity flow direction")
    dinf_dir = rd.FlowDirection(filled_dem, method='Dinf')
    logger.info(
        "Direction range: %.4f to %.4f rad (nodata = %s)",
        float(dinf_dir[dinf_dir != dinf_dir.no_data].min()),
        float(dinf_dir[dinf_dir != dinf_dir.no_data].max()),
        dinf_dir.no_data,
    )
    return dinf_dir

Step 4 — Compute flow accumulation

import numpy as np

def compute_accumulation(filled_dem: rd.rdarray) -> rd.rdarray:
    """Return upslope contributing area per cell as float64."""
    logger.info("Computing D-Infinity flow accumulation")
    acc = rd.FlowAccumulation(filled_dem, method='Dinf')
    arr = np.array(acc)
    logger.info(
        "Accumulation stats — min: %.1f, max: %.1f, mean: %.1f cells",
        arr[arr >= 0].min(), arr.max(), arr[arr >= 0].mean(),
    )
    return acc

For guidance on translating the accumulation raster into a channel network, see stream threshold tuning for systematic calibration of flow accumulation thresholds.

Step 5 — Threshold and vectorize streams

import geopandas as gpd
import rasterio.features as feats
from shapely.geometry import shape
from skimage.morphology import skeletonize

def extract_stream_network(
    acc: rd.rdarray,
    dem_path: str,
    threshold_m2: float = 1_000_000,
) -> gpd.GeoDataFrame:
    """Threshold accumulation, skeletonize, and vectorize as a GeoDataFrame."""
    with rasterio.open(dem_path) as src:
        cell_area_m2 = abs(src.transform.a * src.transform.e)
        transform = src.transform
        crs = src.crs

    threshold_cells = threshold_m2 / cell_area_m2
    logger.info(
        "Threshold: %.0f m² → %.1f cells (cell area = %.2f m²)",
        threshold_m2, threshold_cells, cell_area_m2,
    )

    acc_arr = np.array(acc)
    stream_mask = (acc_arr >= threshold_cells).astype(np.uint8)

    # Skeletonize to single-pixel-wide centre lines before vectorizing
    skeleton = skeletonize(stream_mask).astype(np.uint8)
    logger.info("Stream mask: %d cells → %d skeleton cells", stream_mask.sum(), skeleton.sum())

    shapes = list(feats.shapes(skeleton, mask=skeleton, transform=transform))
    geometries = [shape(s) for s, v in shapes if v == 1]
    gdf = gpd.GeoDataFrame(geometry=geometries, crs=crs)
    logger.info("Vectorized %d stream segments", len(gdf))
    return gdf

Production-Ready Code

The function below integrates all steps, includes logging throughout, preserves raster metadata, and writes a vectorized stream network to GeoPackage. Pass tile_window to process large DEMs in tiles.

import logging
import numpy as np
import rasterio
import rasterio.features as feats
import richdem as rd
import geopandas as gpd
from pathlib import Path
from shapely.geometry import shape
from skimage.morphology import skeletonize

logger = logging.getLogger(__name__)


def run_dinfinity_pipeline(
    dem_path: str,
    output_gpkg: str,
    threshold_m2: float = 1_000_000,
) -> gpd.GeoDataFrame:
    """
    End-to-end D-Infinity flow routing pipeline.

    Parameters
    ----------
    dem_path : str
        Path to a hydrologically conditioned, projected DEM (GeoTIFF).
    output_gpkg : str
        Output path for the vectorized stream network (GeoPackage).
    threshold_m2 : float
        Upslope area threshold in square metres for stream initiation.

    Returns
    -------
    gpd.GeoDataFrame
        Stream network with CRS matching the input DEM.
    """
    logger.info("=== D-Infinity pipeline start: %s ===", dem_path)

    # --- 1. Validate input ---
    with rasterio.open(dem_path) as src:
        if not src.crs.is_projected:
            raise ValueError(
                f"DEM must use a projected CRS. Found: {src.crs}"
            )
        if src.nodata is None:
            raise ValueError("DEM must define a nodata value.")
        nodata = src.nodata
        cell_area_m2 = abs(src.transform.a * src.transform.e)
        transform = src.transform
        crs = src.crs
        logger.info("Input CRS: %s | cell area: %.2f m²", crs, cell_area_m2)

    # --- 2. Load into richdem ---
    dem = rd.LoadGDAL(dem_path)
    dem.no_data = nodata

    # --- 3. Condition ---
    logger.info("Filling depressions")
    filled = rd.FillDepressions(dem, in_place=False)
    logger.info("Resolving flats")
    rd.ResolveFlats(filled, in_place=True)

    # --- 4. D-Infinity direction ---
    logger.info("Computing D-Infinity flow direction")
    _dinf_dir = rd.FlowDirection(filled, method='Dinf')

    # --- 5. Flow accumulation ---
    logger.info("Computing D-Infinity flow accumulation")
    acc = rd.FlowAccumulation(filled, method='Dinf')
    acc_arr = np.array(acc, dtype=np.float64)
    logger.info(
        "Accumulation max: %.0f cells (%.2f km²)",
        acc_arr.max(),
        acc_arr.max() * cell_area_m2 / 1e6,
    )

    # --- 6. Stream mask and skeletonization ---
    threshold_cells = threshold_m2 / cell_area_m2
    stream_mask = (acc_arr >= threshold_cells).astype(np.uint8)
    skeleton = skeletonize(stream_mask).astype(np.uint8)
    logger.info(
        "Stream cells before/after skeletonize: %d / %d",
        stream_mask.sum(), skeleton.sum(),
    )

    # --- 7. Vectorize ---
    raw_shapes = list(feats.shapes(skeleton, mask=skeleton, transform=transform))
    geometries = [shape(s) for s, v in raw_shapes if v == 1]
    if not geometries:
        logger.warning("No stream segments extracted. Consider lowering threshold_m2.")
    gdf = gpd.GeoDataFrame(geometry=geometries, crs=crs)
    gdf["source_dem"] = str(Path(dem_path).name)
    gdf["threshold_m2"] = threshold_m2

    # --- 8. Write output ---
    out_path = Path(output_gpkg)
    out_path.parent.mkdir(parents=True, exist_ok=True)
    gdf.to_file(str(out_path), layer="streams_dinf", driver="GPKG")
    logger.info("Stream network written to %s (%d segments)", out_path, len(gdf))

    logger.info("=== D-Infinity pipeline complete ===")
    return gdf


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    run_dinfinity_pipeline(
        dem_path="data/dem_filled.tif",
        output_gpkg="outputs/streams_dinf.gpkg",
        threshold_m2=500_000,
    )

Validation Protocol

After running the pipeline, verify correctness at three levels before using the output in modeling.

Direction raster checks

import numpy as np
import richdem as rd

def validate_direction_raster(dinf_dir: rd.rdarray) -> dict:
    """Assert D-Infinity direction values fall in [0, 2π] and report nodata fraction."""
    arr = np.array(dinf_dir, dtype=np.float64)
    nd = dinf_dir.no_data
    valid = arr[arr != nd]
    assert valid.min() >= 0.0, "Negative direction values found — possible unresolved flats."
    assert valid.max() <= 2 * np.pi + 1e-6, "Direction values exceed 2π — check nodata masking."
    nodata_frac = (arr == nd).mean()
    logger.info("Valid direction cells: %.1f%% | nodata: %.1f%%",
                (1 - nodata_frac) * 100, nodata_frac * 100)
    return {"valid_pct": (1 - nodata_frac) * 100, "nodata_pct": nodata_frac * 100}

A nodata fraction above 5% typically signals unresolved flat areas; re-run rd.ResolveFlats or inspect the DEM for embedded void regions.

Accumulation sanity checks

• The maximum accumulation cell should correspond to the basin outlet or a major trunk channel. Verify its location by overlaying on a hillshade or reference imagery.
• Sum the total accumulation across all cells and compare to the expected contributing area (DEM cell count × cell area). D-Infinity totals will be slightly lower than D8 totals because fractional routing can leave residual area on ridge cells.
• Compare extracted stream length against National Hydrography Dataset (NHD) polylines. A discrepancy greater than 15% suggests the threshold needs adjustment. See stream threshold tuning for systematic calibration.

Topological validation

def check_stream_topology(gdf: gpd.GeoDataFrame) -> None:
    """Warn about disconnected segments that indicate routing errors."""
    from shapely.ops import unary_union
    merged = unary_union(gdf.geometry)
    if hasattr(merged, 'geoms'):
        n_parts = len(list(merged.geoms))
        logger.warning(
            "Stream network has %d disconnected components. "
            "Check for unresolved flats or excessive nodata near divides.",
            n_parts,
        )
    else:
        logger.info("Stream network is fully connected.")

Common Failure Modes & Optimization

Flat-area direction undefined. Cells in zero-gradient zones receive nodata direction values. The fix is rd.ResolveFlats applied after rd.FillDepressions, never before. Running them in the wrong order means filled depressions re-introduce flats that go unresolved.

Artificial interior depressions from road embankments or bridges. These appear in LiDAR-derived DEMs as enclosed pits that fill but do not drain naturally. Apply a breaching algorithm (rd.BreachDepressions) before filling to cut a channel through the obstruction and preserve the natural flow path. See best practices for filling sinks in high-resolution LiDAR data for a combined breaching-and-filling strategy.

Memory exhaustion on large DEMs. A 1-metre DEM covering a large watershed can exceed 50 GB uncompressed. Use tiled processing: split the DEM into overlapping tiles (≥500 m buffer overlap to prevent edge artifacts), run the full pipeline on each tile independently, and mosaic accumulation results by summing shared boundary cells. Store each tile’s direction raster as float32 and accumulation as float64.

Coastal and tidal boundaries acting as infinite sinks. Mask marine areas using a shoreline polygon before calling FillDepressions. Without masking, the entire coastal fringe drains into the ocean cell and accumulation values along the coast become meaningless.

CRS mismatch between DEM and output vector. Always re-read the CRS from the original raster file when creating the output GeoDataFrame rather than hardcoding an EPSG code. See coordinate reference system alignment for a systematic CRS validation workflow.

Threshold too high hides headwater channels. If the extracted network shows only trunk channels, lower threshold_m2. For steep terrain and ephemeral stream detection, thresholds in the range 10 000–100 000 m² are often appropriate; for low-relief agricultural catchments 500 000–2 000 000 m² is typical.

Accumulation stored as int32 truncates large watersheds. The rd.FlowAccumulation return type from richdem is float32 by default when using the D-Infinity method. Cast to float64 before applying thresholds to avoid precision loss on basins exceeding ~16 million cells.

When to Use D-Infinity vs. Alternatives

Use D-Infinity when:

• Your analysis area includes divergent hillslope flow (ridges, alluvial fans, talus slopes) where routing all flow to a single neighbor would be physically incorrect.
• You are working with high-resolution LiDAR (≤2 m) where subtle slope gradients create real bifurcated flow paths.
• You need continuous upslope area estimates for terrain wetness index, RUSLE slope-length factor, or similar indices that require fractional area distribution.

Consider D8 Flow Direction Implementation when:

• Your target output is a discrete drainage network topology (stream nodes and links) that requires integer-count accumulation for threshold logic.
• Computational speed is critical and the terrain is predominantly convergent (valley-dominated) — D8 runs significantly faster and produces equivalent stream networks in most valley settings.
• Downstream modeling tools (HEC-RAS geometry, SWMM connectivity) expect single-receiver topology.

Consider Multiple Flow Direction Methods when:

• You need to distribute flow across all downslope neighbors, not just two, for applications such as dispersion modeling or soil moisture redistribution.
• You are removing flat-area artifacts from an existing flow direction grid without recomputing the full routing from scratch.

Related Topics

• Comparing D8 vs D-Infinity for Steep Terrain Hydrology — side-by-side benchmarks and terrain-type decision criteria
• D8 Flow Direction Implementation — single-direction deterministic routing in Python
• Multiple Flow Direction Methods — splitting flow across all downslope neighbors
• Stream Threshold Tuning — calibrating flow accumulation thresholds for ephemeral and perennial channels
• DEM Pit Filling Algorithms — sink resolution strategies before routing

Frequently Asked Questions

When should I choose D-Infinity over D8?

Choose D-Infinity when your DEM captures divergent hillslope flow, on ridges, or in complex terrain where routing all upslope area to a single downstream cell creates artificial channelization. D8 remains appropriate for clearly convergent drainage networks and when computational speed is the primary constraint.

What data type should the D-Infinity accumulation array use?

Always cast to float64 after accumulation. Fractional flow partitioning means cell values are non-integer, and large watersheds can produce contributing areas whose total cell count exceeds int32 limits (~2.1 billion).

Why does D-Infinity produce undefined directions in flat areas?

The algorithm selects the facet with the maximum downward slope; zero-gradient zones have no steepest direction and the result is undefined. Apply rd.ResolveFlats after depression filling to impose a minimal drainage gradient across flat regions before computing direction.