⚙️ How-To Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 15 of 20 Updated

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 70.5% of mapped streams of order 3 and above lay within 30 m of a derived channel. Lowering the mapped lines by 10 m raised that to 87.4%; WhiteboxTools' fill_burn raised it to 96.6%.
  • The catchment barely changed: 490.46 km² unburned, 491.48 km² burned, against 497.3 km² published; overlap with the USGS basin polygon rose from IoU 0.9794 to 0.9840.
  • A 2 m burn did little; 50 m aligned more than 10 m but gave the same catchment: 70.5%, 74.0%, 87.4% and 96.1% of major mapped streams were found with no burn and with 2, 10 and 50 m, while the catchment was 491.48 km² at 10 m and 491.50 km² at 50 m.
  • On the 10 m 3DEP DEM burning was nearly redundant: 94.0% of order-3-and-above streams were already within 30 m, 97.5% after fill_burn.

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 the mapped network, not just the catchment area.

Bar chart of the share of mapped streams of order 3 and above within 30 metres of a derived channel on GLO-30 and 3DEP, unburned, with a 10 metre burn and with WhiteboxTools fill_burn.
Burning fixed alignment on the 30 m surface model and changed little on the 10 m bare-earth DEM.

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, 97.1% 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. Burning it lowered 56,214 GLO-30 cells by up to 32.6 m along a route no surface water takes. It happened not to move the gauge's catchment, but 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 almost as scattered as before: 86.9% of channel cells within 30 m of a mapped line and 74.0% of major mapped streams found. 10 m brought 96.1% and 87.4%; 50 m reached 100% and 96.1%, as good as fill_burn, without changing the catchment. 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, fills depressions and routes consistently in one step, and took 0.4 s on GLO-30 and 3.7 s on 3DEP. It changed 68,472 GLO-30 cells and gave the best alignment of any method: 96.6% of mapped streams of order 3 and above within 30 m of a derived channel.

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, burning moved the Coldbrook catchment by 1 km² and raised IoU with the USGS basin polygon from 0.9794 to 0.9840; alignment is what changed.

Table of burn depths on GLO-30 — none, 2, 10 and 50 metres and fill_burn — with catchment area, IoU with the USGS basin and the share of derived channel cells within 30 metres of mapped lines.
Past 10 m, deeper burning improved alignment but left the catchment unchanged.

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_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 490.46 km2 (-1.4%), IoU with the USGS basin 0.9794
burn 2 m   catchment 490.49 km2 (-1.4%), IoU with the USGS basin 0.9795
burn 10 m  catchment 491.48 km2 (-1.2%), IoU with the USGS basin 0.9840
burn 50 m  catchment 491.50 km2 (-1.2%), IoU with the USGS basin 0.9840

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.4 s, 68,472 cells changed; catchment 491.48 km2, IoU 0.9840

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 84.2%; mapped within 30 m of derived: all orders 34.9%, order 3+ 70.5%
burn 2 m   derived cells within 30 m of mapped 86.9%; mapped within 30 m of derived: all orders 37.1%, order 3+ 74.0%
burn 10 m  derived cells within 30 m of mapped 96.1%; mapped within 30 m of derived: all orders 42.8%, order 3+ 87.4%
burn 50 m  derived cells within 30 m of mapped 100.0%; mapped within 30 m of derived: all orders 45.0%, order 3+ 96.1%
fill_burn  derived cells within 30 m of mapped 100.0%; mapped within 30 m of derived: all orders 45.4%, order 3+ 96.6%

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, most of which lie far from the channels being burned. Burning moves where water runs inside the catchment, which changes which cells count as channels, but rarely which side of a ridge a cell drains to. On GLO-30 the gauge catchment moved by 0.2%; alignment on order-3-and-above streams rose by 26 percentage points.

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, 34.9% of mapped length was within 30 m of a derived channel unburned and 45.4% 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 burns and then fills and routes so the burned network drains continuously downstream; on GLO-30 it matched the alignment of a 50 m trench, 96.6% against 96.1% of major mapped streams.

Two panels comparing GLO-30 and 3DEP 10 m: unburned alignment, alignment after fill_burn, catchment area and cells changed by burning.
Burning is a correction for a coarse or noisy DEM; a good DEM needs little of it.

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.

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 70.5% to 87.4%.

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 87.4% of major mapped streams and 50 m 96.1%, with the same catchment. Deeper burns override the DEM more, so use them only where the mapped network is reliable.

Does stream burning change watershed area?

Usually only slightly. The Esopus catchment on GLO-30 went from 490.46 to 491.48 km²; burning mainly changes where channels run inside the catchment.

Should I include pipelines when burning streams?

No. The Shandaken Tunnel in NHDPlus HR would have lowered 56,214 GLO-30 cells by up to 32.6 m along an underground route.

Is stream burning needed for a 10 m DEM?

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

What does WhiteboxTools FillBurn do?

It burns a vector stream network into a DEM and fills the result so the network drains continuously. It took 0.4 s on the 30 m DEM and aligned 96.6% of major mapped streams.