Burning Stream Networks into DEMs with the AGREE Method

When a grid genuinely cannot resolve a channel — a 30 m global DEM over a narrow incised valley, or a low-relief coastal plain where the channel is shallower than the vertical accuracy — no amount of depression filling recovers the drainage. Stream burning forces the issue by editing the elevation surface so the mapped network becomes the steepest path. This guide implements the AGREE method, which is the version worth using; it belongs to the hydro-conditioning for culverts and road crossings topic within hydrology data preparation and DEM processing, and it should be reached for only after breaching has been ruled out.

Prerequisites

Beyond a projected DEM and rasterio, this technique needs:

• A vector stream network you are prepared to treat as ground truth, because the extracted network will reproduce it.
• scipy.ndimage for the distance transform that drives the smooth buffer.
• A decision, made in advance, about what you will validate against — the burned network cannot validate itself.

Core Technique: Two Drops, Not One

The name AGREE comes from the original ArcInfo implementation, and its insight is that a burn needs two components with different geometries.

The sharp drop applies only to the cells the flowline passes through. It guarantees that once flow reaches the channel it stays there, because no neighbouring cell can be lower.

The smooth drop applies to a buffer around the channel and decays linearly with distance from it. It is what makes the surrounding hillslope drain into the channel rather than alongside it. Without it, the burn produces a trench with vertical walls, and the cells beside that trench route parallel to it — the classic artefact of naive burning.

Annotated Code Example

The implementation is short because scipy.ndimage.distance_transform_edt does the geometric work. Each line that is not obvious is commented.

import logging

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize
from scipy.ndimage import distance_transform_edt

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


def agree_burn(
    dem_path: str,
    streams_path: str,
    out_path: str,
    smooth_drop_m: float = 1.5,
    sharp_drop_m: float = 8.0,
    buffer_cells: int = 4,
) -> dict:
    """
    Apply an AGREE stream burn to a DEM.

    Parameters
    ----------
    dem_path      : Input DEM, projected CRS, square cells.
    streams_path  : Vector stream network, any GDAL-readable format.
    out_path      : Destination for the burned DEM.
    smooth_drop_m : Total drop applied at the channel by the smooth buffer.
    sharp_drop_m  : Additional drop applied to the channel cells alone.
    buffer_cells  : Width of the smooth buffer, measured in cells either side.

    Returns
    -------
    Summary dict with the cell counts touched by each component.
    """
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype(np.float32)
        profile = src.profile.copy()
        transform, crs = src.transform, src.crs
        nodata = src.nodata if src.nodata is not None else -9999.0

    streams = gpd.read_file(streams_path).to_crs(crs)

    # --- Rasterize the network. all_touched keeps a diagonal flowline
    # connected; without it a 45-degree reach becomes a dotted line. ---
    channel = rasterize(
        ((geom, 1) for geom in streams.geometry if geom is not None),
        out_shape=dem.shape, transform=transform, fill=0,
        dtype=np.uint8, all_touched=True,
    ).astype(bool)
    log.info("Rasterized %d stream features to %d channel cells",
             len(streams), int(channel.sum()))

    valid = (dem != nodata) & np.isfinite(dem)

    # --- Distance, in cells, from every cell to the nearest channel cell.
    # The transform runs on the INVERSE mask: distance to the nearest True. ---
    dist_cells = distance_transform_edt(~channel)

    # --- Smooth component: full drop at the channel, tapering linearly to
    # zero at the buffer edge. Cells beyond the buffer are untouched. ---
    ramp = np.clip(1.0 - dist_cells / float(buffer_cells), 0.0, 1.0)
    smooth = ramp * smooth_drop_m
    n_buffer = int((ramp > 0).sum())

    burned = dem.copy()
    burned[valid] -= smooth[valid].astype(np.float32)

    # --- Sharp component: an extra drop on the channel cells only, applied
    # after the smooth pass so the two are additive at the centreline. ---
    burned[channel & valid] -= np.float32(sharp_drop_m)

    total_drop = float(smooth_drop_m + sharp_drop_m)
    log.info("Smooth component touched %d cells (%.2f %% of the grid)",
             n_buffer, 100.0 * n_buffer / valid.sum())
    log.info("Channel cells dropped by %.2f m in total", total_drop)

    profile.update(dtype="float32", nodata=nodata, compress="lzw", tiled=True)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(burned, 1)
    log.info("Wrote burned DEM: %s", out_path)

    return {
        "channel_cells": int(channel.sum()),
        "buffer_cells_touched": n_buffer,
        "total_channel_drop_m": total_drop,
        "output": out_path,
    }


# --- Example usage ---
# agree_burn(
#     dem_path="basin_30m.tif",
#     streams_path="nhd_flowlines.gpkg",
#     out_path="basin_30m_burned.tif",
#     smooth_drop_m=1.5, sharp_drop_m=8.0, buffer_cells=4,
# )

Parameter Reference

The relationship between the two drops is the parameter choice that matters. A large sharp drop with no smooth buffer is the trench artefact. A large smooth buffer with a small sharp drop produces a broad shallow valley that competes with real terrain a few cells away.

Worked Example: Reading the Output

On a 30 m DEM over a 900 km² coastal-plain basin with a 1.5 m smooth drop, an 8 m sharp drop and a 4-cell buffer:

The agreement figure of 0.97 is not evidence of anything: it is what burning guarantees. The number that carries information is the delineated area against the published gauge area, which moved from 15 % low to under 1 % — and that comparison is independent of the burned network.

The 6.8 % of cells with materially changed slope is the price. Any analysis on the burned surface that depends on slope, elevation or depth is affected in that fraction of the basin, which is why a burned DEM should be kept as a separate product from the conditioned DEM used for terrain analysis.

Keep the burned surface as a separate product. The moment it is reused for terrain analysis, the buffer’s smooth drop shows up as a valley that is wider and shallower than the real one.

Gotchas and Edge Cases

• Validating against the burn source. The agreement is constructed, not observed. Use gauge drainage areas, an independent channel map or field survey.
• all_touched=False. Leaves a dotted channel on diagonal reaches. The extracted network then fragments exactly where the burn was supposed to help.
• Burning across a divide. A mapped flowline with a digitising error that crosses a ridge will burn a channel through it, capturing a neighbouring basin. Inspect the burned surface for new channels crossing contours at right angles.
• Burning before breaching. A burn passes straight through a road embankment because the flowline does, which happens to be the desired outcome — but it does so by lowering the road, and any subsequent analysis sees a road that is metres below its real elevation.
• Reusing the burned DEM downstream. Slope, curvature, wetness index and any depth calculation are all distorted near the channel. Keep the burned surface as an input to routing only.
• Distance transform units. distance_transform_edt returns distance in cells unless a sampling argument is given. A buffer specified in metres and applied in cells is off by the cell size, which on a 30 m grid is a factor of thirty.

Related Topics

• Hydro-Conditioning for Culverts and Road Crossings — the parent topic, and the breaching approach that should be preferred where the DEM can resolve the channel
• Breaching Road Embankments and Culverts with WhiteboxTools — the local, auditable alternative to a global burn
• Validating Hydro-Conditioning Against NHD Flowlines — how to check a conditioned surface without the circularity a burn introduces
• DEM Pit Filling Algorithms — the fill pass that follows any burn