Watershed Modeling Python for hydrology automation

Breaching vs Filling: Choosing a Depression Treatment

Every depression in a DEM is one of three things: a data artefact, an artificial blockage, or a real closed basin. The three call for different treatments, and the usual practice of applying one algorithm to all of them is why so many conditioned surfaces are wrong in ways nobody notices. This guide is the decision, as part of the DEM pit filling algorithms topic within hydrology data preparation and DEM processing.

Prerequisites

  • A DEM with its depressions labelled, so their count, area and volume are known rather than assumed.
  • Some knowledge of the landscape: is it glaciated, karst, arid, heavily roaded, or none of those.
  • A stated purpose for the conditioned surface, because a DEM conditioned for routing and one conditioned for terrain analysis are not the same product.

Core Technique: What Each Treatment Does to the Surface

Filling raises terrain until every depression spills. Breaching lowers terrain along a path until every depression drains. Both make the surface routable, and they distribute the damage very differently.

Three Treatments of One Depression A depression in a valley. Filling raises the surface across the whole depression to the spill elevation. Breaching cuts a narrow channel from the depression's low point through the confining rim. Leaving it in place preserves both the depression and its storage volume. fill terrain raised across the whole depression millions of cells changed storage destroyed breach one channel cut through the confining rim tens of thousands of cells storage mostly preserved keep and model storage no cells changed; the basin routes only spill the only correct answer where the depressions are real

Parameter Reference: The Decision Table

Landscape / cause Best treatment Why
Sensor noise pits, sub-metre LiDAR Fill, small epsilon The depressions are artefacts; there is nothing to preserve
Road and rail embankments Breach at the crossings The blockage is artificial and the crossing location is known
Sinkholes on karst Keep and model The depression is the hydrology; routing through it is fiction
Prairie potholes Keep, with a fill-spill model Contributing area varies with wetness; a static answer is wrong
Playas and closed arid basins Keep They genuinely do not contribute to the outlet
Quarries, borrow pits Keep Real closed basins, even though artificial
Vegetation-induced pits in canopy DEMs Fill Artefacts of surface, not ground
Unknown cause, regional 10–30 m DEM Fill Most depressions at that resolution are artefacts

The row that most often gets misapplied is the third and fourth: in a landscape with genuine closed basins, filling inflates the delineated contributing area at the outlet by everything those basins hold, which on a pothole landscape can be tens of percent.

Annotated Code Example

The function below classifies depressions before treating them, and applies a different rule to each class.

python 
import logging

import numpy as np
import rasterio
import richdem as rd
from scipy import ndimage

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


def classify_depressions(
    dem: np.ndarray,
    cell_size: float,
    nodata: float = -9999.0,
    artefact_max_cells: int = 12,
    artefact_max_depth_m: float = 0.35,
) -> dict:
    """
    Label every depression and sort it into artefact, blockage or real basin.

    The two thresholds separate sensor noise from everything else. Anything
    larger or deeper than a noise pit needs a decision rather than a default.
    """
    grid = rd.rdarray(dem.copy().astype(np.float32), no_data=nodata)
    rd.FillDepressions(grid, epsilon=False, in_place=True)
    raised = np.asarray(grid) - dem

    labels, n = ndimage.label(raised > 1e-6)
    log.info("Found %d depressions", n)
    if n == 0:
        return {"n": 0, "artefact": [], "candidate": []}

    sizes = ndimage.sum(np.ones_like(labels), labels, index=range(1, n + 1))
    depths = ndimage.maximum(raised, labels, index=range(1, n + 1))
    volumes = ndimage.sum(raised, labels, index=range(1, n + 1)) * cell_size ** 2

    artefact, candidate = [], []
    for i in range(n):
        rec = {"label": i + 1, "cells": int(sizes[i]),
               "depth_m": float(depths[i]), "volume_m3": float(volumes[i])}
        if rec["cells"] <= artefact_max_cells and rec["depth_m"] <= artefact_max_depth_m:
            artefact.append(rec)
        else:
            candidate.append(rec)

    art_vol = sum(r["volume_m3"] for r in artefact)
    cand_vol = sum(r["volume_m3"] for r in candidate)
    log.info("%d artefact pits (%.0f m³ total) and %d larger depressions "
             "(%.0f m³) needing a decision",
             len(artefact), art_vol, len(candidate), cand_vol)
    if candidate:
        biggest = max(candidate, key=lambda r: r["volume_m3"])
        log.info("Largest depression: %d cells, %.2f m deep, %.0f m³ — "
                 "inspect this one before choosing a treatment",
                 biggest["cells"], biggest["depth_m"], biggest["volume_m3"])
    return {"n": n, "artefact": artefact, "candidate": candidate,
            "labels": labels, "raised": raised}


def treatment_impact(dem_path: str, cell_size: float | None = None) -> dict:
    """
    Report what each treatment would cost, before choosing one.

    Runs both a fill and a breach and measures the cells modified, the storage
    volume removed and the change in the maximum accumulation — the three
    numbers the decision actually turns on.
    """
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype(np.float32)
        nodata = src.nodata if src.nodata is not None else -9999.0
        cell_size = cell_size or abs(src.transform.a)

    base = rd.rdarray(dem.copy(), no_data=nodata)
    acc_before = float(np.asarray(rd.FlowAccumulation(
        rd.FillDepressions(base, epsilon=True), method="D8")).max())

    filled = rd.rdarray(dem.copy(), no_data=nodata)
    rd.FillDepressions(filled, epsilon=True, in_place=True)
    fill_changed = int((np.asarray(filled) - dem > 1e-6).sum())
    fill_volume = float((np.asarray(filled) - dem).sum() * cell_size ** 2)

    breached = rd.rdarray(dem.copy(), no_data=nodata)
    rd.BreachDepressions(breached, in_place=True)
    breach_changed = int((np.abs(np.asarray(breached) - dem) > 1e-6).sum())
    breach_volume = float(np.abs(np.minimum(
        np.asarray(breached) - dem, 0.0)).sum() * cell_size ** 2)

    log.info("Fill:   %d cells changed, %.0f m³ of storage removed",
             fill_changed, fill_volume)
    log.info("Breach: %d cells changed, %.0f m³ of terrain excavated",
             breach_changed, breach_volume)
    log.info("Breaching touches %.1f%% as many cells as filling",
             100.0 * breach_changed / max(1, fill_changed))
    return {
        "fill_cells": fill_changed, "fill_volume_m3": fill_volume,
        "breach_cells": breach_changed, "breach_volume_m3": breach_volume,
        "max_accumulation": acc_before,
    }


# --- Example usage ---
# impact = treatment_impact("basin_1m.tif")

Worked Example: The Same DEM, Three Ways

A 1 m LiDAR DEM over a 40 km² basin with mixed agriculture, roads and a few farm ponds:

Treatment Cells changed Storage removed Delineated area at the outlet Max accumulation
Fill everything 4 210 000 1.82 M m³ 40.1 km² 40.1 M cells
Breach everything 68 400 0.11 M m³ 39.8 km² 39.8 M cells
Breach roads, keep ponds 41 200 0.02 M m³ 38.6 km² 38.6 M cells

The third row is the defensible one and it is 3.7 % smaller than the first. That difference is the farm ponds: filling them routes their catchments to the outlet, which is what happens in a large storm and not what happens the rest of the time. Which row is right depends on whether the model is a design flood or a water balance — and the point is that the choice is visible and stated, rather than made by a default.

The Size of the Edit Each Treatment Makes Filling everything modifies 4.2 million cells and removes 1.82 million cubic metres of storage. Breaching everything modifies 68 thousand cells and removes 0.11 million. Breaching roads while keeping ponds modifies 41 thousand cells and removes 0.02 million. fill everything 4.2 M breach everything 68 400 breach roads, keep ponds 41 200 DEM cells modified (log scale) A hundredfold difference in how much of the original surface survives. 40 km² basin, 1 m LiDAR

The three treatments also differ in what they leave available downstream. A surface conditioned for routing has had its storage removed, so it cannot answer a storage question afterwards.

What Each Treatment Leaves You Able to Ask A matrix of three treatments against four downstream questions. Filling supports routing and delineation but not storage volume or inundation extent. Breaching supports all four with minor caveats. Keeping the depressions supports storage questions and needs an explicit spill model for routing. treatment routing delineation storage volume inundation fill everything yes yes no no breach structures yes yes mostly mostly keep and model with a spill model yes yes yes Keep the unconditioned DEM: it is the only product that can answer the last two columns exactly.

Gotchas and Edge Cases

  • One treatment for all depressions. The commonest error. Classify first; the classification takes seconds.
  • Filling a pothole landscape. Inflates the contributing area by everything the potholes hold, which on some prairie basins exceeds a third of the nominal area.
  • Breaching a dam. Drains a reservoir that exists. Guard with a maximum-cost limit and a check against mapped waterbodies.
  • Conditioning once for every purpose. A surface conditioned for routing has had its storage removed; using it for a depression-storage analysis returns zero by construction. Keep both products.
  • Epsilon left off. A filled depression becomes a perfectly flat surface, and the router assigns arbitrary directions across it.
  • Not measuring the impact. Both treatments produce a routable surface, and only the cell counts and area changes reveal how much they differ. Measure before choosing.