How to Burn a Known River Network into a DEM

Problem statement

Channels derived from a DEM do not always sit where the rivers are. A 30 m surface model puts tree canopy on valley floors, smooths narrow gorges and lets flow cut across meander necks; the derived stream wanders tens of metres from the mapped one, and some tributaries join the wrong river. Stream burning lowers the DEM along a mapped network so that flow routing follows it. It is a blunt tool: it fixes alignment where the mapped network is right and imposes errors where it is not.

Measured on the Esopus Creek basin, New York, burning the USGS NHDPlus High Resolution network (805.7 km of streams, artificial paths and connectors in the basin) into the 30 m Copernicus GLO-30 DEM and the 10 m 3DEP DEM, with streams defined at 1 km²:

• On GLO-30, only 76.4% of mapped streams of order 3 and above lay within 30 m of a derived channel, and the catchment was 448.62 km², 9.8% short. Lowering the mapped lines by 10 m raised alignment to 91.8% and put the catchment back at 491.73 km² (IoU 0.985 with the USGS basin).
• A 2 m burn did nothing for the catchment; 50 m aligned more than 10 m but gave a slightly worse one: 448.83 km² at 2 m, 491.73 km² at 10 m, 489.84 km² at 50 m.
• WhiteboxTools' fill_burn aligned channels best of all, 96.0%, and broke the catchment: 70.01 km² at the gauge, because 3.6 km upstream the main river crossed into a neighbouring burned line.
• On the 10 m 3DEP DEM burning was nearly redundant: 93.5% of order-3-and-above streams were already within 30 m, 97.2% after fill_burn, and the catchment stayed at 492.24 km² whatever was burned.

Quick answer

import numpy as np
import rasterio
from rasterio.features import rasterize

with rasterio.open("dem.tif") as src:
    profile, dem = src.profile, src.read(1)
lines = rasterize(streams.to_crs(profile["crs"]).geometry, out_shape=dem.shape, transform=profile["transform"]).astype(bool)
burned = np.where(lines & (dem != profile["nodata"]), dem - 10, dem).astype(profile["dtype"])
with rasterio.open("dem_burned.tif", "w", **profile) as dst:
    dst.write(burned, 1)

Then condition and route the burned DEM as usual. Burn only surface flow paths — leave out pipelines and tunnels — and check the result against both the mapped network and the catchment area.

no burn    catchment 448.62 km2 (-9.8%), IoU with the USGS basin 0.8965
burn 2 m   catchment 448.83 km2 (-9.7%), IoU with the USGS basin 0.8969
burn 10 m  catchment 491.73 km2 (-1.1%), IoU with the USGS basin 0.9846
burn 50 m  catchment 489.84 km2 (-1.5%), IoU with the USGS basin 0.9808

Step-by-step solution

1. Decide whether you need to burn

Burn when derived channels must match a mapped network: to attach model results to mapped reaches, to compute distance to a known river, or where the DEM is too coarse or too noisy to find channels. Measure the alignment first (Example 3). On the 10 m 3DEP DEM, 96.9% of derived channel cells were within 30 m of a mapped line without any burning.

2. Choose and clean the network

Use the lines water actually follows on the surface. NHDPlus HR in the basin included 8.54 km of pipeline, part of the Shandaken Tunnel that carries water between reservoirs underground. Passing it to fill_burn with the rest changed 30,015 more GLO-30 cells, sinking a tunnel's route by the same 10,012 m as every other burned line, along a path no surface water takes. It creates channels across ridges. Keep stream/river, artificial path and connector types; drop pipelines.

3. Match the network to the DEM's CRS and extent

Reproject the lines to the DEM's CRS and clip them to its extent. Lines that end exactly at the edge are fine; lines outside it are ignored by rasterisation.

4. Rasterise the lines

rasterize marks every cell a line passes through; the default marks only cells whose centres are touched by the line's rendering, all_touched=True marks more. A one-cell-wide channel is what D8 routing needs.

5. Lower the lines, not the whole valley

Subtract a fixed depth. On GLO-30, 2 m left derived channels scattered: 89.1% of channel cells within 30 m of a mapped line and 79.8% of major mapped streams found, and the catchment still 9.7% short. 10 m brought 98.0% and 91.8% and recovered the catchment; 50 m reached 100% and 96.4%, but the catchment fell back to 489.84 km². The deeper the burn, the more the mapped network overrides the DEM, and very deep burns create artificial canyons that distort slope and HAND; go deep only where the network is trustworthy.

6. Or use a tool that burns and fills together

WhiteboxTools' fill_burn burns the network and fills the result in one step, and took 0.3 s on GLO-30 and 19.0 s on 3DEP. It lowers the burned lines by about 10,000 m: on GLO-30 it changed 65,985 cells, lowering one by 10,012 m. Every mapped line becomes a trench so deep that nothing else matters, and where two lines touch in the raster the trenches join. Alignment was the best of any method, 96.0% of major mapped streams — but on GLO-30 the Esopus main stem's flow crossed into an adjacent burned artificial path about 3.6 km above the Coldbrook gauge and left the DEM at its east edge, so the gauge's catchment came out as 70.01 km². On 3DEP the same tool left the catchment at 492.24 km². Check the catchment after any burn, not just the alignment.

7. Condition and route the burned DEM

A burned DEM still has pits and flats. Fill or breach it, then compute flow directions and accumulation as for any DEM; see conditioning a DEM.

8. Check alignment and area

Compare two ways: how much of the derived network sits on mapped lines, and how much of the mapped network is found by derived channels. Check the catchment at a gauge too. On GLO-30, a 10 m burn moved the Coldbrook catchment from 448.62 to 491.73 km² and IoU with the USGS basin polygon from 0.8965 to 0.9846; fill_burn moved it to 70.01 km².

Code examples

Example 1 — burn at several depths and route each DEM

import os

import geopandas as gpd
import numpy as np
import rasterio
import whitebox
from pyproj import Transformer
from rasterio.features import rasterize
from shapely.geometry import Point

wbt = whitebox.WhiteboxTools()
wbt.set_max_procs(1)                      # one thread: parallel runs give slightly different results
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("."))
gx, gy = Transformer.from_crs("EPSG:4269", "EPSG:32618", always_xy=True).transform(-74.2701944, 42.0144722)
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(32618)
flowlines = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(32618)
network = flowlines[flowlines.ftype.isin([460, 558, 334])]        # stream/river, artificial path, connector

with rasterio.open("esopus_glo30_utm.tif") as src:
    profile, dem, transform = src.profile, src.read(1), src.transform
cell = transform.a
in_basin = rasterize(basin.geometry, out_shape=dem.shape, transform=transform).astype(bool)


def evaluate(dem_path, tag):
    wbt.fill_depressions(dem_path, f"{tag}_filled.tif", fix_flats=True)
    wbt.d8_pointer(f"{tag}_filled.tif", f"{tag}_d8.tif")
    wbt.d8_flow_accumulation(f"{tag}_d8.tif", f"{tag}_acc.tif", out_type="cells", pntr=True)
    with rasterio.open(f"{tag}_acc.tif") as src:
        acc = src.read(1)
    col, row = ~transform * (gx, gy)
    r, c, k = int(row), int(col), int(150 / cell)
    window = acc[r - k:r + k + 1, c - k:c + k + 1]
    dr, dc = np.unravel_index(window.argmax(), window.shape)
    gpd.GeoDataFrame(geometry=[Point(*(transform * (c - k + dc + 0.5, r - k + dr + 0.5)))], crs=32618).to_file(f"{tag}_outlet.shp")
    wbt.watershed(f"{tag}_d8.tif", f"{tag}_outlet.shp", f"{tag}_ws.tif")
    with rasterio.open(f"{tag}_ws.tif") as src:
        ws = src.read(1)
        ws = (ws != src.nodata) & (ws > 0)
    streams = (acc * cell * cell / 1e6 >= 1.0) & in_basin
    return streams, window.max() * cell * cell / 1e6, (ws & in_basin).sum() / (ws | in_basin).sum()


lines = rasterize(network.geometry, out_shape=dem.shape, transform=transform).astype(bool) & (dem != profile["nodata"])
results = {"no burn": evaluate("esopus_glo30_utm.tif", "raw")}
for depth in (2, 10, 50):
    with rasterio.open(f"burn_{depth}m.tif", "w", **profile) as dst:
        dst.write(np.where(lines, dem - depth, dem).astype(profile["dtype"]), 1)
    results[f"burn {depth} m"] = evaluate(f"burn_{depth}m.tif", f"burn{depth}")
for name, (_, km2, iou) in results.items():
    print(f"{name:10} catchment {km2:.2f} km2 ({km2 / 497.3 - 1:+.1%}), IoU with the USGS basin {iou:.4f}")

no burn    catchment 448.62 km2 (-9.8%), IoU with the USGS basin 0.8965
burn 2 m   catchment 448.83 km2 (-9.7%), IoU with the USGS basin 0.8969
burn 10 m  catchment 491.73 km2 (-1.1%), IoU with the USGS basin 0.9846
burn 50 m  catchment 489.84 km2 (-1.5%), IoU with the USGS basin 0.9808

Example 2 — WhiteboxTools fill_burn

import time

network[["geometry"]].to_file("nhd_network.shp")
start = time.perf_counter()
wbt.fill_burn("esopus_glo30_utm.tif", "nhd_network.shp", "fill_burn.tif")
seconds = time.perf_counter() - start
with rasterio.open("fill_burn.tif") as src:
    changed = (np.abs(src.read(1) - dem) > 1e-3) & (dem != profile["nodata"])
results["fill_burn"] = evaluate("fill_burn.tif", "fillburn")
_, km2, iou = results["fill_burn"]
print(f"fill_burn {seconds:.1f} s, {changed.sum():,} cells changed; catchment {km2:.2f} km2, IoU {iou:.4f}")

fill_burn 0.3 s, 65,985 cells changed; catchment 70.01 km2, IoU 0.0072

Example 3 — how well derived channels find the mapped network

from scipy.ndimage import distance_transform_edt

mapped = gpd.clip(flowlines[flowlines.ftype.isin([460, 558])], basin.geometry.make_valid())
samples = {}
for lowest in (1, 3):
    points = [part.interpolate(d) for g in mapped[mapped.streamorde >= lowest].geometry
              for part in getattr(g, "geoms", [g]) for d in np.arange(0, part.length, 10.0)]
    samples[lowest] = (((np.array([p.y for p in points]) - transform.f) / -cell).astype(int),
                       ((np.array([p.x for p in points]) - transform.c) / cell).astype(int))
mapped_cells = rasterize(mapped.geometry, out_shape=dem.shape, transform=transform, all_touched=True).astype(bool)
to_mapped = distance_transform_edt(~mapped_cells) * cell

for name, (streams, _, _) in results.items():
    to_derived = distance_transform_edt(~streams) * cell
    print(f"{name:10} derived cells within 30 m of mapped {np.mean(to_mapped[streams] <= 30):.1%}; mapped within 30 m of derived: "
          f"all orders {np.mean(to_derived[samples[1]] <= 30):.1%}, order 3+ {np.mean(to_derived[samples[3]] <= 30):.1%}")

no burn    derived cells within 30 m of mapped 86.2%; mapped within 30 m of derived: all orders 37.0%, order 3+ 76.4%
burn 2 m   derived cells within 30 m of mapped 89.1%; mapped within 30 m of derived: all orders 38.9%, order 3+ 79.8%
burn 10 m  derived cells within 30 m of mapped 98.0%; mapped within 30 m of derived: all orders 44.1%, order 3+ 91.8%
burn 50 m  derived cells within 30 m of mapped 100.0%; mapped within 30 m of derived: all orders 45.2%, order 3+ 96.4%
fill_burn  derived cells within 30 m of mapped 100.0%; mapped within 30 m of derived: all orders 45.5%, order 3+ 96.0%

The mapped network is sampled every 10 m, and each sample is checked against the distance to the nearest derived channel cell.

Explanation

Why derived channels miss mapped rivers

GLO-30 is a surface model: over forest its elevations include canopy, and a 30 m cell averages across narrow valley floors. The lowest cells in a valley may be on one bank, and flow can cut corners across meanders. The derived channel is then consistently offset from the mapped centreline — a median of 30 m, one cell, on GLO-30 — and, where a ridge is smoothed away, joins the wrong valley.

Why burning changes alignment more than area

Catchment area is set by divides, and most divides lie far from the channels being burned. But where a surface model has smoothed a low divide away, a whole tributary can drain the wrong way, and then a burn that forces the mapped network moves the divide back. That is what happened on GLO-30: unburned, the catchment was 43 km² short; a 10 m burn restored it to within 1.1% of the published area. A burn can move a catchment in the other direction just as easily — fill_burn's joined trenches sent the main river past the gauge — so area is the check that catches both.

Why low-order recall stays low

The comparison counts mapped streams of every order, and NHDPlus HR maps many first-order streams draining less than 1 km². Channels defined at 1 km² cannot find them whatever the DEM: across all orders, 37.0% of mapped length was within 30 m of a derived channel unburned and 45.5% after fill_burn. The order-3-and-above figure is the fair measure of burning.

Why fill_burn beats a simple subtraction

A fixed-depth burn creates a trench whose bottom still follows the DEM's noisy profile, and conditioning afterwards has to fill the trench's own pits. fill_burn lowers the lines so far — about 10 km — that the trench's own profile no longer matters, and fills what is left; on GLO-30 it matched the alignment of a 50 m trench, 96.0% against 96.4% of major mapped streams. The cost of that depth is that trenches which touch become one channel, whichever way the river really goes.

Edge cases or notes

• Outdated networks impose old channels where rivers have moved; check the network's date against the DEM's.
• Pipelines, tunnels and canals are not surface flow paths; exclude them unless you are modelling managed transfers.
• Braided or looping lines create multiple burned paths; D8 can follow only one.
• Network gaps at culverts or confluences leave burned segments that do not connect; snap and dissolve lines first.
• Deep burns distort slope, HAND and anything else computed from the burned DEM; keep the unburned DEM for those.
• Walls — raising the DEM along mapped divides — are the complement of burning where catchment boundaries are known.
• Resolution mismatch: a network mapped at 1:24,000 is finer than a 30 m DEM can represent; burning forces it anyway.

Internal links

• How to extract a stream network from a DEM — measuring derived against mapped streams
• How to condition a DEM for hydrology: breaching or filling — what to run after burning
• How DEM resolution and source change a drainage network — why GLO-30 needs burning more
• How to delineate a watershed from a pour point in Python — catchments from burned DEMs
• Fixing streams that stop at roads, bridges and embankments — burning at crossings only
• How to calculate height above nearest drainage (HAND) in Python — a product distorted by deep burns
• Stream order explained: Strahler, Shreve and what they measure — the orders used in the comparison
• How to rasterize a vector layer in Python — turning lines into cells

FAQ

What is stream burning in a DEM?

Lowering the DEM's elevations along a mapped river network so that flow routing follows the mapped channels. Lowering NHDPlus HR lines by 10 m in GLO-30 raised aligned major streams from 76.4% to 91.8%.

How deep should I burn streams into a DEM?

Deep enough to beat the relief across the valley floor. On GLO-30, 2 m was too little; 10 m found 91.8% of major mapped streams and 50 m 96.4%, and 10 m gave the better catchment. Deeper burns override the DEM more, so use them only where the mapped network is reliable.

Does stream burning change watershed area?

It can, in either direction. On GLO-30 a 10 m burn took the Esopus catchment from 448.62 to 491.73 km² by restoring a divide the surface model had lost, and fill_burn took it to 70.01 km² by joining two burned lines above the gauge. On the 10 m 3DEP DEM no burn changed it at all.

Should I include pipelines when burning streams?

No. The Shandaken Tunnel in NHDPlus HR would have added 30,015 burned GLO-30 cells along an underground route.

Is stream burning needed for a 10 m DEM?

Often not. On the 10 m 3DEP DEM, 93.5% of major mapped streams were already within 30 m of a derived channel; fill_burn raised it to 97.2%.

What does WhiteboxTools FillBurn do?

It burns a vector stream network into a DEM, about 10,000 m deep, and fills the result so the network drains continuously. It took 0.3 s on the 30 m DEM and aligned 96.0% of major mapped streams — but joined touching lines, which cost the Coldbrook gauge most of its catchment. Check the catchment afterwards.