Removing flat area artifacts from flow direction grids

Zero-gradient cells are the single most common preprocessing failure in Python hydrology workflows. When a flow direction algorithm encounters a cell with no downslope neighbor it assigns a null direction code, breaking accumulation matrices and corrupting watershed boundaries. This page covers the specific richdem calls — rd.FillDepressions followed by rd.ResolveFlats — that eliminate those artifacts before any routing is computed. It sits within the Multiple Flow Direction Methods workflow, which is part of the broader Flow Routing Algorithms & Stream Network Extraction topic area.

Prerequisites

This page assumes familiarity with the Multiple Flow Direction Methods stack: Python 3.9+, richdem>=0.3.4, and rasterio>=1.3.0. Beyond those base requirements, you need:

• A single-band float32 or float64 GeoTIFF projected to a metric CRS (UTM or an equal-area projection). Angular projections produce distorted slope values that undermine flat-detection heuristics.
• A consistent nodata value set in the raster metadata. richdem reads this via dem.no_data; mismatched or absent nodata causes edge cells to be treated as valid terrain.
• Optional: a permanent water mask as a boolean raster aligned to the DEM, to exclude lakes and reservoirs before resolution.

pip install richdem>=0.3.4 rasterio>=1.3.0 numpy>=1.24.0

Why flat areas break routing

Flat areas enter a DEM through three routes: LiDAR pulse interpolation across smooth surfaces, DEM resampling that clips decimal precision, and genuine near-level terrain such as floodplains, playas, and constructed detention ponds. In each case the result is a contiguous region where every cell shares the same elevation value, leaving the routing algorithm with no slope signal to follow.

For D8 flow direction implementation this produces null direction codes (typically encoded as 0 or −1 depending on the library). For the Freeman (1991) MFD weight formula used in Multiple Flow Direction Methods, the consequences are worse: the weight calculation divides by the sum of downslope tangent slopes, which is zero across a flat, producing a divide-by-zero that propagates NaNs through the entire accumulation solve.

The three-stage fix is:

1. Fill depressions — removes closed basins so every cell eventually drains to the DEM boundary.
2. Resolve flats — assigns a micro-gradient to each flat zone that is consistent with the direction water would have to flow to exit the zone.
3. Compute direction — runs the chosen routing algorithm on the corrected surface.

The diagram below illustrates how a single filled depression creates a flat plateau and how ResolveFlats imposes a gradient toward the outlet.

Core technique explanation

richdem.ResolveFlats implements the Barnes, Lehman & Mulla (2014) epsilon-gradient algorithm. Rather than adding a fixed perturbation — which can over-correct on very gentle terrain — it assigns the minimum gradient that propagates consistently toward each flat zone’s outlet. The magnitude is determined by the floating-point epsilon of the DEM’s dtype, so it survives dtype-preserving I/O without rounding to zero.

The algorithm works in two passes over each contiguous flat region:

1. Gradient toward higher ground — assigns increasing labels (one per cell) moving away from the flat’s inlet edges, so cells near the high-side boundary drain inward.
2. Gradient toward lower ground — assigns increasing labels moving away from the outlet edges, so cells near the outlet drain outward.

The two label sets are added together. The result is a monotonically increasing surface from inlet to outlet, guaranteed to carry D8 and MFD routing across without null directions.

FillDepressions must run first. It uses the Wang & Liu (2006) priority-flood algorithm (O(n log n)) to raise every closed basin to the elevation of its lowest spill point. This step creates new flat plateaus at the fill level — which is exactly what ResolveFlats is then applied to. For a deeper treatment of depression-filling strategies, see DEM pit filling algorithms.

Annotated code example

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

logger = logging.getLogger(__name__)


def remove_flat_artifacts(
    dem_path: str,
    output_dem_path: str,
    output_fd_path: str,
    fd_method: str = "D8",
    water_mask_path: str | None = None,
) -> tuple[str, str]:
    """
    Fill depressions, resolve flat areas, and write a clean flow direction raster.

    Parameters
    ----------
    dem_path : str
        Path to input DEM (single-band float32/float64 GeoTIFF, metric CRS).
    output_dem_path : str
        Where to write the resolved DEM.
    output_fd_path : str
        Where to write the flow direction raster.
    fd_method : str
        'D8' or 'Dinf' (D-Infinity). Passed to richdem.FlowDirection.
    water_mask_path : str | None
        Optional boolean raster; True cells are set to nodata before flat resolution
        to prevent artificial routing through open water.

    Returns
    -------
    tuple[str, str]
        Paths to the resolved DEM and flow direction raster.
    """
    dem_path = Path(dem_path)
    if not dem_path.exists():
        raise FileNotFoundError(f"DEM not found: {dem_path}")

    # ── 1. Read spatial metadata via rasterio ─────────────────────────────────
    with rasterio.open(dem_path) as src:
        profile = src.profile.copy()
        # Fall back to -9999 if the raster has no declared nodata.
        nodata = src.nodata if src.nodata is not None else -9999.0
        logger.info(
            "Loaded DEM %s: shape=%s CRS=%s nodata=%s",
            dem_path.name,
            (src.height, src.width),
            src.crs,
            nodata,
        )

    # ── 2. Load into richdem ───────────────────────────────────────────────────
    # rd.LoadGDAL reads the GDAL-compatible GeoTIFF and returns an rdarray
    # that carries the affine transform and CRS alongside the numpy data.
    dem_rd = rd.LoadGDAL(str(dem_path))
    dem_rd.no_data = nodata  # override to match rasterio-read nodata

    # ── 3. Apply permanent water mask (optional) ──────────────────────────────
    # Open-water cells are hydrologically flat but should not be routed through.
    # Set them to nodata so ResolveFlats leaves them untouched.
    if water_mask_path is not None:
        wm_path = Path(water_mask_path)
        if not wm_path.exists():
            raise FileNotFoundError(f"Water mask not found: {wm_path}")
        with rasterio.open(wm_path) as wm_src:
            water_mask = wm_src.read(1).astype(bool)
        dem_arr = np.array(dem_rd, dtype=np.float64)
        dem_arr[water_mask] = nodata
        dem_rd[:, :] = dem_arr  # write masked values back into the rdarray
        logger.info(
            "Applied water mask: %d cells set to nodata", int(water_mask.sum())
        )

    # ── 4. Fill depressions ────────────────────────────────────────────────────
    # Wang & Liu (2006) priority-flood guarantees O(n log n) and preserves
    # natural drainage divides; in_place=False returns a new rdarray.
    filled = rd.FillDepressions(dem_rd, in_place=False)
    n_filled = int(np.sum(np.array(filled) > np.array(dem_rd)))
    logger.info("FillDepressions raised %d cells above original elevation.", n_filled)

    # ── 5. Resolve flat areas ─────────────────────────────────────────────────
    # Barnes et al. (2014) epsilon-gradient: assigns the minimum gradient
    # consistent with FP precision across each contiguous flat zone.
    resolved = rd.ResolveFlats(filled, in_place=False)
    logger.info("ResolveFlats complete.")

    # ── 6. Compute flow direction ─────────────────────────────────────────────
    # Normalise user-supplied string so 'dinf', 'D-Inf', etc. all work.
    method_key = "Dinf" if fd_method.upper().replace("-", "") in ("DINF", "DINFINITY") else "D8"
    flow_dir = rd.FlowDirection(resolved, method=method_key)
    logger.info("FlowDirection computed using method=%s.", method_key)

    # ── 7. Validate: count null-direction cells ───────────────────────────────
    fd_arr = np.array(flow_dir)
    null_count = int(np.sum(fd_arr == flow_dir.no_data))
    valid_total = int(np.sum(np.array(dem_rd) != nodata))
    null_pct = 100.0 * null_count / max(valid_total, 1)
    logger.info(
        "Null-direction cells: %d / %d (%.3f%%)", null_count, valid_total, null_pct
    )
    if null_pct > 0.1:
        logger.warning(
            "More than 0.1%% null directions remain (%d cells). "
            "Check for lakes without a water mask or DEM boundary clipping.",
            null_count,
        )

    # ── 8. Write resolved DEM ─────────────────────────────────────────────────
    resolved_arr = np.array(resolved, dtype=np.float32)
    profile.update(dtype="float32", count=1, nodata=nodata)
    with rasterio.open(output_dem_path, "w", **profile) as dst:
        dst.write(resolved_arr, 1)
    logger.info("Resolved DEM written to %s.", output_dem_path)

    # ── 9. Write flow direction raster ────────────────────────────────────────
    fd_profile = profile.copy()
    fd_profile.update(dtype=str(fd_arr.dtype), nodata=flow_dir.no_data)
    with rasterio.open(output_fd_path, "w", **fd_profile) as dst:
        dst.write(fd_arr, 1)
    logger.info("Flow direction raster written to %s.", output_fd_path)

    return output_dem_path, output_fd_path

Parameter reference

Worked example: output interpretation

A correctly resolved D8 output on a 1000 × 1000 DEM should contain:

• Values 1–8 encoding the eight compass directions (richdem encodes 1=E, 2=NE … 8=SE in the standard D8 scheme).
• flow_dir.no_data only at DEM boundary cells and nodata input cells — never interior to the valid terrain.
• Zero residual null interior cells (null_pct < 0.01% in normal terrain).

Run these checks after calling the function above:

import numpy as np
import rasterio
import logging

logger = logging.getLogger(__name__)


def validate_flow_direction(fd_path: str, expected_nodata: float) -> None:
    """Log key statistics for a D8 flow direction raster and warn on anomalies."""
    with rasterio.open(fd_path) as src:
        fd = src.read(1)
        nodata = src.nodata if src.nodata is not None else expected_nodata

    valid = fd[fd != nodata]
    unique_vals = np.unique(valid)
    logger.info("D8 direction codes present: %s", sorted(unique_vals.tolist()))

    if not set(unique_vals).issubset({1, 2, 3, 4, 5, 6, 7, 8}):
        unexpected = set(unique_vals) - {1, 2, 3, 4, 5, 6, 7, 8}
        logger.warning("Unexpected direction codes: %s", unexpected)

    null_interior = int(np.sum(fd == 0))  # richdem uses 0 for unresolved
    if null_interior > 0:
        logger.warning("%d cells carry direction code 0 (unresolved).", null_interior)
    else:
        logger.info("No unresolved (0-coded) cells. Flat resolution successful.")

What anomalies indicate:

Gotchas and edge cases

• Order is non-negotiable. FillDepressions creates flat plateaus at the spill elevation; running ResolveFlats before FillDepressions misses those synthetically created flats. Always fill first.

• float32 can under-resolve large flats. The epsilon assigned by ResolveFlats must be representable in the DEM’s dtype. For flat zones spanning more cells than the exponent range of float32 allows, elevation differences collapse to zero after rounding. Process any flat region wider than roughly 5 km as float64 and downcast only after direction computation.

• Lake masking is not optional for wetland DEMs. Without a permanent water mask, ResolveFlats routes water across open-water surfaces in ways that produce hydrologically meaningless channel networks. Use JRC Global Surface Water or a National Wetlands Inventory layer aligned to your DEM.

• Boundary cells behave differently from interior cells. Cells at the DEM edge have fewer than eight neighbours. richdem assigns these cells a boundary-outlet direction automatically. If your DEM is tiled, add a buffer of at least 10 cells from an adjacent tile before running FillDepressions to prevent artificial edge drainage.

• Parallel tile processing introduces seam artifacts. If you process large DEMs in tiles, flat zones that straddle a tile boundary will be resolved independently and may produce conflicting micro-gradients. Use a minimum 500 m overlap and re-run ResolveFlats across the merged seam strip. The spatial resolution tradeoffs discussion covers tile-based pipeline design in more detail.

Frequently Asked Questions

Why does ResolveFlats need to run after FillDepressions rather than before?

FillDepressions creates new flat plateaus at the fill level. Running ResolveFlats first misses these synthetically created flat zones. The correct order is always fill first, then resolve.

Will the epsilon-gradients distort real topography?

No. The Barnes et al. (2014) algorithm assigns the smallest gradient consistent with floating-point precision — typically sub-millimetre increments on metre-scale DEMs. This is orders of magnitude below real topographic variation and does not affect slope-derived metrics outside the flat zone.

Should genuine water bodies be excluded before resolving flats?

Yes. Apply a permanent water mask and set lake cells to nodata before calling ResolveFlats. Otherwise, the algorithm routes water across open-water surfaces in ways that are hydrologically meaningless.

Related topics

• Multiple Flow Direction Methods — parent page; MFD proportion calculation that depends on this preprocessing step
• D8 Flow Direction Implementation — sibling routing approach that also requires flat-free terrain
• DEM Pit Filling Algorithms — deep treatment of depression-filling strategies and their hydrologic trade-offs
• Best Practices for Filling Sinks in High-Resolution LiDAR Data — LiDAR-specific edge cases for the fill step that precedes flat resolution
• D-Infinity Routing Patterns — alternative routing algorithm for divergent hillslopes that also requires resolved flats