Deriving Topographic Wetness Index for Flood Models

Compute topographic wetness index in Python: fill sinks, accumulate flow, combine upslope area with slope, and avoid the divide-by-zero that ruins flat-terrain flood features.

Topographic wetness index is ln(a / tan(β)), where a is the upslope contributing area per unit contour width and β is the local slope. High values mark places that collect water from a large area and cannot shed it quickly — valley bottoms, footslopes, hollows. It is one of the strongest single predictors available for flood susceptibility, soil moisture and wetland mapping, and it costs one DEM and about forty lines of code.

The index builds directly on the gradients described in computing slope and aspect from a DEM with NumPy, and slots into the wider stack described in DEM and Terrain Derivative Features.

From elevation to wetness index in four operations A raw digital elevation model is depression filled, then flow directions are assigned per cell, then flow is accumulated downslope to give contributing area, and finally the contributing area divided by the tangent of slope is passed through a natural logarithm to produce the topographic wetness index. Each stage depends on the one before it being correct 1. fill sinks no pits left for flow to die in 2. direction D8 or multi-flow per cell 3. accumulate contributing area per cell 4. combine ln(a / tan β) slope floored TWI Skip step 1 and every cell upstream of a pit reports too little contributing area. Skip the slope floor in step 4 and flat terrain returns infinity. Both failures produce a raster that opens fine and trains a worse model.

Why This Fails in Geospatial ML Pipelines

The index has a pole. As slope approaches zero the denominator approaches zero and the ratio diverges, so any DEM containing genuinely flat ground — a floodplain, a lake surface, a reclaimed polder — produces inf after the division and inf after the logarithm. NumPy emits a warning at most. Downstream, inf propagates through scaling, poisons a StandardScaler, and turns an entire feature column into nan for every sample, at which point most estimators either crash or silently drop the column.

The second failure is subtler and more damaging: unfilled depressions. Real DEMs contain thousands of one-cell pits from interpolation noise, plus genuine closed basins. A D8 algorithm routes flow to the steepest downslope neighbour; a pit has none, so flow stops. Every cell draining through that pit reports a contributing area truncated at the pit. The resulting wetness map looks entirely reasonable — it has the right texture and range — but the values along the drainage network, exactly where flood risk concentrates, are systematically too low.

Third, resolution changes what the feature means. Contributing area is measured in cells, and converting to square metres requires the cell size. On a 30 m DEM, a 10,000 m² contributing area is 11 cells; on a 2 m lidar DEM it is 2,500. Train on one and infer on the other and the model sees a feature shifted by two orders of magnitude — the training-serving skew warned about throughout DEM and Terrain Derivative Features.

Core Principles

  • Fill or breach depressions first. Choose deliberately: filling erases genuine basins, breaching cuts artificial channels through them. Record which you used.
  • Floor the slope, never the output. Clamping tan β to a small positive value keeps the physics interpretable; clipping the final index hides the problem.
  • Convert accumulation to real units. Multiply cell counts by cell area, and divide by cell width to get area per unit contour length.
  • Prefer multiple-flow-direction for hillslopes. D8’s single-neighbour routing produces stripes on smooth slopes that a model will happily memorise.
  • Assert finiteness before writing. One assert np.isfinite(twi).all() catches every pole and every void in one line.
  • Pin the DEM resolution in configuration. The feature is only comparable across runs at a fixed cell size.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
import rasterio
from scipy import ndimage

logger = logging.getLogger(__name__)

# Offsets and distances for the eight neighbours, in (row, col) order.
NEIGHBOURS = [(-1, 0), (-1, 1), (0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1)]
MIN_TAN_BETA = np.tan(np.radians(0.1))    # slope floor: 0.1 degrees


def fill_depressions(dem: np.ndarray, max_iter: int = 200) -> np.ndarray:
    """Iterative priority-free depression fill.

    Raises each cell to the minimum of its neighbours' maxima until stable. Simple,
    deterministic, and adequate for the small pits that dominate real DEMs.
    """
    filled = dem.copy()
    for i in range(max_iter):
        neighbour_max = ndimage.maximum_filter(filled, size=3, mode="nearest")
        candidate = np.minimum(neighbour_max, np.nanmax(filled))
        raised = np.where(filled < candidate, np.minimum(candidate, filled + 0.5), filled)
        if np.allclose(raised, filled, equal_nan=True):
            logger.info("depression fill converged after %d passes", i)
            return raised
        filled = raised
    logger.warning("depression fill hit max_iter=%d — large basins may remain", max_iter)
    return filled


def flow_accumulation_mfd(dem: np.ndarray, cell_size: float, p: float = 1.1) -> np.ndarray:
    """Multiple-flow-direction accumulation, in cell counts.

    Flow leaves each cell towards every downslope neighbour, weighted by
    (drop / distance) ** p. Cells are processed from highest to lowest so a cell's
    own inflow is complete before it distributes.
    """
    rows, cols = dem.shape
    acc = np.ones_like(dem, dtype="float64")
    order = np.argsort(dem, axis=None)[::-1]          # highest first

    for flat_idx in order:
        r, c = divmod(int(flat_idx), cols)
        z = dem[r, c]
        if not np.isfinite(z):
            continue

        drops, targets = [], []
        for dr, dc in NEIGHBOURS:
            rr, cc = r + dr, c + dc
            if not (0 <= rr < rows and 0 <= cc < cols):
                continue
            zn = dem[rr, cc]
            if not np.isfinite(zn) or zn >= z:
                continue
            dist = cell_size * (np.hypot(dr, dc))
            drops.append(((z - zn) / dist) ** p)
            targets.append((rr, cc))

        if not drops:
            continue                                   # outlet or filled pit
        weights = np.asarray(drops) / np.sum(drops)
        for (rr, cc), w in zip(targets, weights):
            acc[rr, cc] += acc[r, c] * w
    return acc


def topographic_wetness_index(dem: np.ndarray, cell_size: float,
                              slope_deg: np.ndarray) -> np.ndarray:
    """ln(a / tan beta) with a slope floor and real-world units."""
    filled = fill_depressions(dem)
    acc_cells = flow_accumulation_mfd(filled, cell_size)

    # Specific catchment area: contributing area per unit contour width.
    a = (acc_cells * cell_size ** 2) / cell_size

    tan_beta = np.maximum(np.tan(np.radians(slope_deg)), MIN_TAN_BETA)
    twi = np.log(a / tan_beta).astype("float32")

    if not np.isfinite(twi[np.isfinite(dem)]).all():
        raise ValueError("TWI contains non-finite values over valid terrain")
    return twi

Step-by-Step Walkthrough

Step 1 — Load a metric DEM and its slope

from math import isclose

with rasterio.open("dem_utm33n.tif") as src:
    assert src.crs.is_projected, "reproject to a metric CRS first"
    dem = src.read(1, masked=True).filled(np.nan).astype("float32")
    cell = abs(src.transform.a)
    profile = src.profile.copy()

slope_deg, _ = slope_aspect(dem, cell)          # from the slope and aspect guide

Step 2 — Fill and inspect what changed

filled = fill_depressions(dem)
raised = filled - dem
print(f"{(raised > 0.01).mean():.2%} of cells raised, max {np.nanmax(raised):.2f} m")

If more than a few percent of cells are raised, or the maximum lift exceeds a few metres, the DEM has genuine closed basins. Decide explicitly whether to keep them — an endorheic basin filled flat will report a huge contributing area draining nowhere.

Step 3 — Accumulate and sanity-check the range

acc = flow_accumulation_mfd(filled, cell)
print(f"accumulation: min {acc.min():.1f}, median {np.median(acc):.1f}, max {acc.max():.0f}")
assert acc.min() >= 1.0, "every cell contributes at least itself"

The maximum should be close to the number of cells draining through the catchment outlet. A maximum barely above the median means flow is terminating early — go back to the fill step.

Step 4 — Combine, write, and record the parameters

twi = topographic_wetness_index(dem, cell, slope_deg)
print(f"TWI range {np.nanmin(twi):.2f} to {np.nanmax(twi):.2f}")

profile.update(count=1, dtype="float32", nodata=np.nan, compress="deflate", tiled=True)
with rasterio.open("twi.tif", "w", **profile) as dst:
    dst.write(twi, 1)
    dst.set_band_description(1, "twi")
    dst.update_tags(slope_floor_deg="0.1", flow_algorithm="MFD p=1.1",
                    cell_size_m=str(cell))

Writing the parameters into the file tags is what lets a colleague — or you, in six months — tell whether two wetness rasters are comparable.

An unfilled pit truncates every catchment upstream of it Two rows of cells are compared. In the unfilled row a one-cell depression sits in the middle of the drainage path, so accumulation stops there and downstream cells report only their own local inflow. In the filled row flow passes through and accumulation grows steadily to the outlet. Accumulation along one drainage path unfilled 1 3 8 pit 1 2 3 catchment restarts at 1 filled 1 3 8 14 21 29 38 accumulation grows to the outlet The unfilled map has the right texture and the wrong values exactly where flooding happens.

Verification

import numpy as np


def test_twi_is_finite_on_a_perfect_plane():
    """A plane has constant slope; TWI must be finite everywhere."""
    cell = 10.0
    plane = (np.arange(60, dtype="float32") * cell * 0.05)[None, :].repeat(60, 0)
    slope, _ = slope_aspect(plane, cell)
    twi = topographic_wetness_index(plane, cell, slope)
    assert np.isfinite(twi).all()


def test_flat_terrain_is_capped_not_infinite():
    """Perfectly flat ground must hit the slope floor rather than diverge."""
    flat = np.full((40, 40), 25.0, dtype="float32")
    slope = np.zeros_like(flat)
    twi = topographic_wetness_index(flat, 10.0, slope)
    assert np.isfinite(twi).all()
    assert twi.max() < np.log(40 * 40 * 10.0 / MIN_TAN_BETA) + 1


def test_valley_scores_higher_than_ridge():
    """The whole point of the feature: hollows must beat crests."""
    yy, xx = np.mgrid[0:80, 0:80].astype("float32")
    valley = 100.0 + 0.02 * (xx - 40) ** 2 - 0.4 * yy      # a V running north-south
    slope, _ = slope_aspect(valley, 10.0)
    twi = topographic_wetness_index(valley, 10.0, slope)
    assert twi[40, 40] > twi[40, 5], "valley floor should be wetter than the valley side"

For a visual check, overlay the top decile of the index on a stream network. The high-index cells should trace the channels and the footslopes beside them; if they trace ridge lines instead, the dz_dy sign convention is inverted upstream in the slope computation.

The slope floor is what keeps the index finite A curve plots the wetness index against slope. As slope decreases towards zero the curve rises steeply towards infinity. A vertical line marks the slope floor at one tenth of a degree, and to the left of it the curve is replaced by a flat capped segment, bounding the maximum value the feature can take. ln(a / tan β) has a pole at β = 0 slope β (degrees) → TWI floor 0.1° capped without the floor this region is inf steep ground → low index Record the floor with the raster: it sets the feature’s maximum, so two runs with different floors are not comparable.

FAQ

Why does the index blow up on flat terrain?

Because it divides by tan β, which approaches zero as the ground flattens. Clamp tan β to a small positive floor before dividing — 0.1° is a common choice — and write that floor into the output metadata, since it determines the maximum value the feature can reach.

Do I have to fill sinks before computing flow accumulation?

For single-flow-direction routing, yes: an unfilled pit has no downslope neighbour, so flow terminates and every upstream cell under-reports its contributing area. Fill, or breach if the depressions are real features you want to preserve, and log which choice you made.

Should I use D8 or a multiple-flow-direction algorithm?

D8 is correct in incised channels but produces stripy, artificial-looking accumulation on hillslopes. Multiple-flow-direction spreads flow proportionally to slope and gives a smoother surface that usually predicts better — which is why the implementation above uses it.

Part of: DEM and Terrain Derivative Features Part of: Spatial Feature Engineering for Machine Learning