Watershed Modeling Python for hydrology automation

Breaching Road Embankments and Culverts with WhiteboxTools

Breaching is the conditioning edit that removes an obstruction instead of drowning the ground behind it, and BreachDepressionsLeastCost is the most practical implementation available from Python. It finds, for each depression, the cheapest path out — cheapest meaning the least total elevation that has to be cut — and carves that path rather than raising the depression to its spill point. This guide covers its parameters and the ways it goes wrong, as part of the hydro-conditioning for culverts and road crossings topic within hydrology data preparation and DEM processing.

Prerequisites

  • whitebox installed and able to download its binary on first run (pip install whitebox, then one call to warm the cache).
  • A DEM in a projected CRS with square cells — all distances are specified in cells and converted mentally to metres.
  • Knowledge of the tallest embankment in the area of interest, which sets max_cost.

Core Technique: Least-Cost Breaching

For each depression the tool identifies the pour point, then searches outward for a cell lower than the depression’s minimum. Every candidate path has a cost equal to the total elevation that would have to be removed along it. The cheapest path within the distance limit wins, and the tool carves it.

Three Candidate Paths, One Winner A depression sits behind a road embankment. Path A crosses at the culvert line, cutting 1.4 metres over 18 metres of embankment for a cost of 12. Path B crosses at a higher point for a cost of 41. Path C runs around the end of the embankment for a cost of 9 but a distance of 140 cells, beyond the search radius. depression road A — cost 12, dist 18 cells ✓ chosen B — cost 41, cuts a higher part of the crest C — cost 9, but 140 cells away rejected: beyond the dist limit Cost is total elevation removed along the path, not path length — a long shallow cut can beat a short deep one, which is why dist matters. plan view — depression, embankment and three candidate exits

The behaviour that surprises people is path C. A long path around the end of an embankment can have a lower cost than a short path through it, because cost counts elevation, not distance. The dist parameter is what stops the tool from preferring it — and setting dist too large is how a breach ends up running two hundred metres along a ditch line instead of through the culvert.

Annotated Code Example

python 
import logging
import os

import numpy as np
import rasterio
from whitebox import WhiteboxTools

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


def breach_then_fill(
    dem_path: str,
    out_dir: str,
    dist_cells: int = 40,
    max_cost_m: float = 6.0,
    min_dist: bool = True,
    flat_increment: float = 0.001,
) -> dict:
    """
    Breach what can be breached, then fill the residue.

    Parameters
    ----------
    dem_path       : Input DEM (projected CRS, square cells).
    out_dir        : Working directory for intermediate and final rasters.
    dist_cells     : Maximum breach search radius, in cells.
    max_cost_m     : Maximum total elevation a single path may cut through.
    min_dist       : Prefer the shortest path among those within max_cost.
    flat_increment : Epsilon gradient applied to flats by the fill pass.

    Returns
    -------
    Dict of output paths and the depression counts at each stage.
    """
    os.makedirs(out_dir, exist_ok=True)
    wbt = WhiteboxTools()
    wbt.set_working_dir(os.path.abspath(out_dir))
    wbt.verbose = False           # the tool is extremely chatty by default

    dem_abs = os.path.abspath(dem_path)
    breached = "dem_breached.tif"
    conditioned = "dem_conditioned.tif"
    sinks_before = "sinks_before.tif"
    sinks_after = "sinks_after.tif"

    # --- Count depressions before, so the effect is measurable rather than
    # assumed. Sink() labels each closed depression with a unique id. ---
    wbt.sink(dem_abs, sinks_before)
    n_before = _count_labels(os.path.join(out_dir, sinks_before))
    log.info("Depressions before conditioning: %d", n_before)

    # --- Breach. min_dist=True makes the tool prefer the shortest of the
    # acceptable paths, which keeps a breach on the culvert line rather than
    # letting it wander along a ditch to a marginally cheaper exit. ---
    log.info("Breaching with dist=%d cells, max_cost=%.2f m", dist_cells, max_cost_m)
    wbt.breach_depressions_least_cost(
        dem=dem_abs,
        output=breached,
        dist=dist_cells,
        max_cost=max_cost_m,
        min_dist=min_dist,
        flat_increment=flat_increment,
        fill=False,               # fill separately so each stage is auditable
    )

    wbt.sink(os.path.join(out_dir, breached), sinks_after)
    n_after = _count_labels(os.path.join(out_dir, sinks_after))
    log.info("Depressions after breaching: %d (%.1f %% removed)",
             n_after, 100.0 * (n_before - n_after) / max(1, n_before))

    # --- Fill the residue. Anything the breach could not reach within its
    # limits is either a genuine closed basin or a structure taller than
    # max_cost, and filling is the right treatment for both. ---
    wbt.fill_depressions(
        dem=os.path.join(out_dir, breached),
        output=conditioned,
        fix_flats=True,
        flat_increment=flat_increment,
    )
    log.info("Filled the residual depressions → %s", conditioned)

    # --- Report how much terrain actually moved. A conditioning run that
    # changed millions of cells has flooded something. ---
    changed, max_cut, max_raise = _difference_stats(
        dem_abs, os.path.join(out_dir, conditioned)
    )
    log.info("Cells modified: %d; deepest cut %.2f m; highest fill %.2f m",
             changed, max_cut, max_raise)

    return {
        "breached": os.path.join(out_dir, breached),
        "conditioned": os.path.join(out_dir, conditioned),
        "depressions_before": n_before,
        "depressions_after_breach": n_after,
        "cells_modified": changed,
    }


def _count_labels(path: str) -> int:
    with rasterio.open(path) as src:
        arr = src.read(1)
        nodata = src.nodata
    valid = arr[arr != nodata] if nodata is not None else arr
    return int(np.unique(valid[valid > 0]).size)


def _difference_stats(before_path: str, after_path: str):
    with rasterio.open(before_path) as a, rasterio.open(after_path) as b:
        x = a.read(1).astype("float64")
        y = b.read(1).astype("float64")
        nodata = a.nodata
    m = np.isfinite(x) & np.isfinite(y)
    if nodata is not None:
        m &= x != nodata
    d = np.where(m, y - x, 0.0)
    return int((np.abs(d) > 1e-6).sum()), float(-d.min()), float(d.max())


# --- Example usage ---
# breach_then_fill("basin_1m.tif", "./conditioning", dist_cells=40, max_cost_m=6.0)

Parameter Reference

Parameter Typical value Effect
dist 25–60 cells Must exceed the widest obstruction in cells. Larger values slow the tool sharply and allow implausible paths
max_cost 3–8 m Separates embankments from dams. Depressions needing a deeper cut are left for the fill pass
min_dist True Prefers the shortest acceptable path, keeping the breach on the crossing rather than around it
flat_increment 0.001 m Epsilon gradient on flats, so the breached channel is not a flat surface
fill False Keep the stages separate so each is measurable; run fill_depressions explicitly afterwards
Where max_cost Stops Helping and Starts Draining Reservoirs Two curves against max_cost from 1 to 30 metres. Depressions removed rises steeply to about 92 percent by 6 metres then flattens. Cells modified stays low until about 10 metres then climbs sharply as the tool begins carving through dam crests. usable band 1 m 4 m 10 m 18 m 30 m max_cost — maximum elevation a single breach may cut 0 50 % 100 % depressions removed cells modified — reservoirs being drained Nothing is gained past about 6 m, and past 10 m the tool starts cutting through dam crests.

Worked Example: Output Interpretation

A 1 m LiDAR DEM over a 60 km² rural basin with dist=40, max_cost=6.0:

Stage Depressions Cells modified Deepest cut
Original DEM 18 420 — —
After breaching 1 490 41 300 5.8 m
After fill 0 218 000 —

Read those numbers as a set. Breaching removed 92 % of depressions by modifying 41 000 cells — a tiny fraction of a 60-million-cell grid, which is what a targeted edit should look like. The fill pass then modified far more cells, but almost all of that is the epsilon gradient applied across flats, not real flooding. A deepest cut of 5.8 m sits just under max_cost, which means at least one structure was at the limit and worth inspecting manually.

If the breaching stage had modified two million cells, the cause would be max_cost set high enough to carve through dam crests, and the reservoirs would be gone.

The two parameters interact, and the interaction is what decides whether a run targets structures or reshapes the landscape. Reading them as a grid makes the safe corner obvious.

Where the Two Parameters Are Safe Together A two by three grid of dist against max cost. Small dist with small max cost misses wide embankments. Moderate dist with moderate max cost is the safe combination. Large dist with large max cost both wanders along ditches and cuts through dam crests. misses wide embankments the safe corner dist 25–60, cost 3–8 m breaches wander along ditch lines deep cuts, short reach reservoirs start to drain both failures at once dist < 20 dist 25–60 dist > 100 search radius, cells cost 3–8 m cost > 12 m max_cost

Gotchas and Edge Cases

  • dist in cells, embankments in metres. On a 1 m DEM a 25 m road needs dist above 25; on a 10 m DEM the same road needs dist above 3. Converting in the wrong direction is the most common configuration error.
  • max_cost above the dam height. Reservoirs drain silently. Compare depression volume before and after, split by whether a mapped waterbody sits inside.
  • min_dist=False with a large dist. Breaches wander along ditches to marginally cheaper exits, producing channels that do not correspond to any structure on the ground.
  • Skipping the fill pass. Breaching does not guarantee a depressionless surface. The flow router will produce sinks wherever a depression survived.
  • flat_increment left at zero. The breached channel becomes a flat, and flow direction across it is arbitrary — see removing flat area artifacts from flow direction grids.
  • Relative paths. WhiteboxTools resolves outputs against its working directory but inputs against the current process directory. Passing absolute paths for inputs and bare names for outputs, as above, avoids the most common source of “file not found” on a tool that has already run.