Comparing D8 vs D-Infinity for Steep Terrain Hydrology

On steep terrain — slopes exceeding roughly 15° — D8’s constraint of routing all flow to a single downslope grid neighbor creates stair-step artifacts and artificial flow concentration that do not reflect physical overland-flow behavior. This page focuses on that specific failure mode and on the practical steps for switching to D-Infinity routing in a Python pipeline. It is part of D-Infinity Routing Patterns, which covers the full algorithm, facet geometry, and production implementation. Both pages sit within the broader Flow Routing Algorithms & Stream Network Extraction domain. If you are new to flow routing, work through D-Infinity Routing Patterns first before applying the terrain-specific decisions covered here.

Prerequisites

Beyond the standard prerequisites (richdem >= 0.3, a projected DEM in metres, Python 3.9+), this comparison requires:

• A DEM that has already been conditioned — sinks filled using DEM Pit Filling Algorithms — because un-filled depressions produce truncated accumulation in both algorithms.
• Enough RAM to hold two floating-point accumulation grids simultaneously. For a 10 000 × 10 000 cell DEM at float32, that is approximately 800 MB; plan for at least 2× that to account for intermediate arrays.
• numpy, logging, and optionally rasterio for export — all standard in a GIS Python environment.

Core Technique Explanation

How D8 routes flow on steep terrain

D8 maps each raster cell to exactly one of eight cardinal or diagonal neighbors — the steepest downslope direction. The encoded value (1, 2, 4, 8, 16, 32, 64, or 128) is deterministic and grid-aligned. On gentle slopes this is a reasonable approximation; on steep, convex hillslopes the algorithm forces flow along pixel edges and diagonals that rarely align with the true gradient direction. The result is a “stair-step” pattern in the flow direction raster and artificially high accumulation values along those stair-step paths — the watershed equivalent of aliasing.

How D-Infinity corrects for this

D-Infinity, formalized by Tarboton (1997), calculates the steepest descent direction as a continuous angle θ (in radians, 0 to 2π) across each of the eight triangular facets formed by a cell and its immediate neighbors. It then partitions the cell’s incoming flow proportionally between the two grid neighbors that bracket θ. If θ points exactly between two cells, flow splits 50/50; if it points nearly toward one cell, that cell receives almost all flow. This eliminates grid-alignment bias and allows flow to diverge naturally across micro-topographic features.

The diagram below contrasts D8 single-direction routing with D-Infinity proportional dispersion on a convex hillslope:

Annotated Code Example

The snippet below performs a side-by-side D8 and D-Infinity flow accumulation on the same conditioned DEM and computes a difference raster to identify where the two algorithms diverge most. Every non-obvious line is explained in an inline comment.

import logging
import numpy as np
import richdem as rd

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


def compare_d8_dinf(dem_path: str, output_diff_path: str | None = None) -> dict:
    """
    Run D8 and D-Infinity flow accumulation on a conditioned DEM and return
    comparison statistics. Optionally write the log-transformed difference
    raster to *output_diff_path* via richdem's SaveGDAL.

    Parameters
    ----------
    dem_path : str
        Path to a hydrologically conditioned GeoTIFF (sinks already filled,
        projected CRS in metres).
    output_diff_path : str or None
        If provided, the log10(D8) minus log10(Dinf) difference raster is
        written to this path.

    Returns
    -------
    dict with keys: max_d8, max_dinf, mean_abs_diff, n_cells_diverge_gt10pct
    """
    logger.info("Loading DEM from %s", dem_path)
    dem = rd.LoadGDAL(dem_path)

    # Confirm nodata is set — richdem uses it to mask border cells.
    # If no nodata is present, border cells produce spurious accumulation.
    if dem.no_data is None:
        logger.warning("DEM has no nodata value; border accumulation may be unreliable")

    # Fill any residual sinks. Using in_place=False preserves the original
    # array so we can inspect pre/post differences if needed.
    logger.info("Filling residual depressions (epsilon fill not applied here)")
    dem_filled = rd.FillDepressions(dem, in_place=False)

    # D8: every cell routes 100% of flow to its steepest single neighbor.
    # The method string must be exactly "D8" (case-sensitive in richdem).
    logger.info("Computing D8 flow accumulation")
    acc_d8 = rd.FlowAccumulation(dem_filled, method="D8")

    # D-Infinity: flow is partitioned between the two cells that bracket the
    # steepest facet angle. method="Dinf" (note: capital D, lowercase inf).
    logger.info("Computing D-Infinity flow accumulation")
    acc_dinf = rd.FlowAccumulation(dem_filled, method="Dinf")

    # Convert to numpy for arithmetic. richdem arrays are masked arrays;
    # .filled(0) replaces nodata cells with 0 before comparison.
    arr_d8 = np.array(acc_d8).astype(np.float32)
    arr_dinf = np.array(acc_dinf).astype(np.float32)

    # Log-transform both to compress the dynamic range for comparison.
    # log1p(x) = log(1+x) avoids log(0) on cells with zero accumulation.
    log_d8 = np.log1p(arr_d8)
    log_dinf = np.log1p(arr_dinf)

    diff = log_d8 - log_dinf  # positive = D8 concentrates more than D-Infinity

    # Count cells where the algorithms diverge by more than 10% of the D8 value.
    valid_mask = arr_d8 > 0
    pct_diff = np.abs(arr_d8 - arr_dinf) / (arr_d8 + 1e-6)
    n_diverge = int(np.sum((pct_diff > 0.10) & valid_mask))

    stats = {
        "max_d8": float(arr_d8.max()),
        "max_dinf": float(arr_dinf.max()),
        "mean_abs_diff": float(np.mean(np.abs(diff[valid_mask]))),
        "n_cells_diverge_gt10pct": n_diverge,
    }
    logger.info("Comparison stats: %s", stats)

    if output_diff_path:
        # Wrap the numpy array back into a richdem-compatible array by copying
        # the metadata from the original dem object before saving.
        diff_rd = rd.rdarray(diff, no_data=-9999.0)
        diff_rd.geotransform = dem.geotransform
        diff_rd.projection = dem.projection
        rd.SaveGDAL(output_diff_path, diff_rd)
        logger.info("Difference raster saved to %s", output_diff_path)

    return stats


if __name__ == "__main__":
    results = compare_d8_dinf(
        dem_path="steep_terrain_filled.tif",
        output_diff_path="d8_dinf_log_diff.tif"
    )
    print(results)

Parameter Reference Table

Worked Example: Output Interpretation

Running the code above on a 5 m LiDAR DEM of a steep volcanic catchment (mean slope ~22°, 4 km²) produces results like:

max_d8: 1 842 300    # cell count; multiply by cell_area (25 m²) = 46 km² — matches catchment area
max_dinf: 1 190 500  # lower peak because flow is spread across paths, not concentrated
mean_abs_diff: 0.41  # in log units; large divergence confirms steep-terrain D8 bias
n_cells_diverge_gt10pct: 182 440  # ~45% of valid cells differ by >10% between algorithms

What these numbers tell you:

• max_d8 >> max_dinf: D8 artificially concentrates accumulation along its stair-step paths. Applying a stream threshold calibrated on D8 data to a D-Infinity raster would extract far fewer channels.
• mean_abs_diff approaching 0.5 log units indicates that the two methods diverge significantly across the catchment — a clear signal that D8 is poorly suited to this terrain and D-Infinity should be used for stream threshold tuning.
• Spatial pattern of diff: Positive values (D8 > D-Infinity) concentrate along ridge lines and convex breaks-of-slope — exactly where D8 stair-stepping is worst. Negative values (D8 < D-Infinity) are rare and indicate areas where D8 happened to align with the true gradient.

When difference values are near zero everywhere, the terrain is gentle enough that algorithm choice has minimal impact and D8’s lower memory overhead is preferable.

The diagram below summarises the decision logic for choosing between the two methods based on terrain characteristics:

Gotchas & Edge Cases

• Stream threshold recalibration is mandatory. D-Infinity accumulation peaks are typically 20–50% lower than equivalent D8 peaks in steep terrain. If you copy a D8-calibrated threshold into a D-Infinity workflow, your extracted channel network will be sparse or absent. Re-derive thresholds using field-measured drainage density or high-resolution orthoimagery, as described in tuning flow accumulation thresholds for ephemeral streams.

• D-Infinity and flat areas interact poorly. After sink filling, large flat areas (lakes, agricultural fields) produce undefined descent angles for D-Infinity because no triangular facet has a positive slope. richdem handles this via flat-routing heuristics, but those heuristics can create unexpected accumulation patterns. Check the output for anomalous bull’s-eye patterns in flat regions and apply removing flat area artifacts from flow direction grids if needed.

• Memory doubles relative to D8. D-Infinity stores a floating-point angle per cell rather than an integer encoding. For a 50 000 × 50 000 DEM, D-Infinity direction storage alone requires ~10 GB. Process in tiles or use out-of-core methods for continental-scale datasets.

• CRS errors invalidate slope-based decisions. Computing a slope threshold to select between algorithms requires that your DEM uses a projected CRS in metres. Geographic CRS (degrees) produces incorrect slope magnitudes. Validate with rasterio before routing; see fixing CRS mismatches in watershed shapefiles for correction workflows.

• D-Infinity can oversmooth in deeply incised gorges. In narrow, high-relief canyons where the true flow path is strictly confined, D-Infinity’s dispersion spreads accumulation to adjacent walls, underestimating the main channel’s contributing area. In these zones, D8 or a multiple flow direction method that enforces channel confinement may better match physical reality.

Related Topics

• D-Infinity Routing Patterns — full D-Infinity algorithm, triangular facet geometry, and production implementation
• Stream Threshold Tuning — recalibrating accumulation thresholds after switching routing methods
• Tuning Flow Accumulation Thresholds for Ephemeral Streams — practical threshold derivation workflow
• Multiple Flow Direction Methods — Quinn MFD and other dispersive alternatives to D-Infinity
• D8 Flow Direction Implementation — D8 in depth, for cases where single-direction routing is appropriate