⚙️ How-To Hydrology & watersheds Terrain & elevation Geometry & topology Hydrology & Watersheds · step 7 of 20 Updated

How to Extract a Stream Network from a DEM

Problem statement

A DEM stream network is every cell whose flow accumulation exceeds a threshold, turned into lines. The code is short. The threshold is not a detail: it decides how many streams exist, how long the network is, how far up each valley the channels start, and what Strahler order the main river ends up with. There is no single correct value, only values that match a reference network better or worse.

Measured on the conditioned 10 m DEM of the 492 km² Esopus Creek catchment in New York, against the USGS NHDPlus High Resolution flowlines (805.7 km of mapped streams in the basin):

  • Thresholds from 0.05 to 10 km² gave networks from 1,260 km to 114 km long — 5,095 segments down to 27.
  • A 0.1 km² threshold matched the mapped network's density best: 1.76 km of stream per km² against 1.64 for NHDPlus HR, with 85% of the derived network within 50 m of a mapped stream and 89% of the mapped network within 50 m of a derived one.
  • At 1 km², 99% of the derived network lay on mapped streams but it covered only 45% of them.
  • The main river's Strahler order at the outlet went from 7 to 4 across the same range of thresholds.

Quick answer 

import numpy as np

if not hasattr(np, "in1d"):
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
import geopandas as gpd
from pysheds.grid import Grid
from shapely.geometry import shape

grid = Grid.from_raster("esopus_breached.tif")
fdir = grid.flowdir(grid.read_raster("esopus_breached.tif"))
acc = grid.accumulation(fdir)

threshold_km2 = 0.1
cells = threshold_km2 * 1e6 / abs(grid.affine.a * grid.affine.e)
network = grid.extract_river_network(fdir, acc > cells)
streams = gpd.GeoDataFrame(geometry=[shape(f["geometry"]) for f in network["features"]], crs=grid.crs.srs)
print(len(streams), "segments,", round(streams.length.sum() / 1000, 1), "km")

On the whole DEM, before clipping to the catchment, that printed 4049 segments, 1458.7 km. Express the threshold as an area, not a cell count, so it means the same thing at any resolution. Then compare the result with a mapped network before choosing.

Bar chart of total derived stream length in the Esopus basin for accumulation thresholds from 0.05 to 10 km2, with the NHDPlus High Resolution length for comparison.
The threshold sets the network's length over a factor of eleven; the mapped network sits near 0.1 km².

Step-by-step solution

1. Route flow on a conditioned DEM

Stream extraction needs flow direction and accumulation from a depression-free DEM; see flow accumulation. Clip to the catchment of interest so lengths and densities refer to one basin: here the catchment above the Coldbrook gauge, 492.02 km².

2. Choose the threshold as an area

A threshold of 1,000 cells is 0.1 km² on a 10 m grid and 0.9 km² on a 30 m grid. Convert from an area to cells for each DEM. The area is the channel-initiation area: how much upslope land it takes, in this landscape and this DEM, to start a stream.

3. Sweep thresholds and measure the network

On the Esopus DEM:

threshold   segments   length     density (km/km2)   max Strahler   derived within 50 m of mapped
0.05 km2      5,095   1,259.9 km   2.561                7              71%
0.1 km2       2,393     865.0 km   1.758                7              85%
0.25 km2        946     583.2 km   1.185                6              95%
0.5 km2         501     443.6 km   0.902                6              98%
1 km2           257     332.8 km   0.676                5              99%
2 km2           139     242.8 km   0.494                5              98%
5 km2            59     162.5 km   0.330                4              98%
10 km2           27     114.2 km   0.232                4              99%

Extraction and stream ordering together took 0.5–6.3 s per threshold, apart from the first, which spent 15.8 s mostly compiling pysheds' routines.

4. Compare with a mapped network both ways

Two numbers matter. Precision: how much of the derived network lies near a mapped stream. Recall: how much of the mapped network lies near a derived one. At 0.1 km², 85% of derived length was within 50 m of NHDPlus HR and 89% of NHDPlus HR length within 50 m of the derived network. At 1 km² the figures were 99% and 45%: every derived stream was real, and more than half the real ones were missing. Pick the threshold where both are acceptable (Example 2).

5. Match the reference's level of detail

Reference networks differ as much as thresholds do. NHDPlus HR mapped 805.7 km of flowlines in the basin, a density of 1.64 km/km²; the medium-resolution NHDPlus V2 network mapped 148.1 km, 0.30 km/km². A threshold tuned to one would be badly wrong against the other. State which reference you matched.

6. Remove what the reference counts differently

NHDPlus HR in the basin included 128.8 km of artificial paths — flowlines through lakes and wide rivers — and pipelines that carry water between reservoirs underground. Artificial paths are real flow routes; pipelines are not surface streams and should be excluded before comparing.

7. Check where channels start, not just how long they are

A threshold that matches total length can still start channels in the wrong places: on steep, convergent slopes streams begin with smaller areas than on gentle, divergent ones. Map the derived channel heads over the reference; if they are systematically too far up or down valleys, a single threshold may not suit the whole basin.

8. Record the threshold with the network

Every downstream product — stream order, stream density, distance to stream, HAND — depends on it. Keep it in the output's attributes or file name (Example 3).

Bar chart of the share of derived stream length within 50 m of mapped streams and the share of mapped stream length within 50 m of derived streams, for thresholds of 0.1, 0.25, 0.5 and 1 km2.
Raising the threshold makes the derived network more correct and less complete.

Code examples

Example 1 — sweep thresholds 

import time

cell_km2 = abs(grid.affine.a * grid.affine.e) / 1e6
x, y = grid.snap_to_mask(acc > 1e6 / cell_km2 / 1e6 * 1.0, (560427.2, 4651640.7))
catchment = grid.catchment(x=x, y=y, fdir=fdir, xytype="coordinate")
grid.clip_to(catchment)
fdir_c, acc_c = grid.view(fdir), grid.view(acc)
basin_km2 = float(np.asarray(catchment).astype(bool).sum()) * cell_km2

networks = {}
for threshold_km2 in (0.05, 0.1, 0.25, 0.5, 1, 2, 5, 10):
    start = time.perf_counter()
    network = grid.extract_river_network(fdir_c, acc_c > threshold_km2 / cell_km2)
    lines = gpd.GeoDataFrame(geometry=[shape(f["geometry"]) for f in network["features"]], crs=grid.crs.srs)
    order = np.asarray(grid.stream_order(fdir_c, acc_c > threshold_km2 / cell_km2))
    networks[threshold_km2] = lines
    length_km = lines.length.sum() / 1000
    print(f"{threshold_km2:>5} km2: {len(lines):>5} segments, {length_km:7.1f} km, {length_km / basin_km2:.3f} km/km2, "
          f"max Strahler {order.max()}, {time.perf_counter() - start:.2f} s")
 0.05 km2:  5095 segments,  1259.9 km, 2.561 km/km2, max Strahler 7, 15.80 s
  0.1 km2:  2393 segments,   865.0 km, 1.758 km/km2, max Strahler 7, 6.28 s
 0.25 km2:   946 segments,   583.2 km, 1.185 km/km2, max Strahler 6, 1.50 s
  0.5 km2:   501 segments,   443.6 km, 0.902 km/km2, max Strahler 6, 2.43 s
    1 km2:   257 segments,   332.8 km, 0.676 km/km2, max Strahler 5, 1.91 s
    2 km2:   139 segments,   242.8 km, 0.494 km/km2, max Strahler 5, 1.03 s
    5 km2:    59 segments,   162.5 km, 0.330 km/km2, max Strahler 4, 0.66 s
   10 km2:    27 segments,   114.2 km, 0.232 km/km2, max Strahler 4, 0.46 s

The outlet is snapped to a cell draining more than 1 km² so the catchment is the whole basin; see pour points explained.

Example 2 — precision and recall against a mapped network 

mapped = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(grid.crs.srs)
mapped = mapped[mapped.ftype.isin([460, 558, 334])]          # stream/river, artificial path, connector; no pipelines
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(grid.crs.srs)
mapped = gpd.clip(mapped, basin)
mapped_zone = mapped.geometry.union_all().buffer(50)

for threshold_km2 in (0.1, 0.25, 0.5, 1):
    derived = networks[threshold_km2]
    precision = derived.geometry.intersection(mapped_zone).length.sum() / derived.length.sum()
    recall = mapped.geometry.intersection(derived.geometry.union_all().buffer(50)).length.sum() / mapped.length.sum()
    print(f"{threshold_km2:>4} km2: derived near mapped {precision:.0%}, mapped near derived {recall:.0%}")
print(f"mapped network: {len(mapped)} flowlines, {mapped.length.sum() / 1000:.1f} km")
 0.1 km2: derived near mapped 85%, mapped near derived 89%
0.25 km2: derived near mapped 95%, mapped near derived 74%
 0.5 km2: derived near mapped 98%, mapped near derived 59%
   1 km2: derived near mapped 99%, mapped near derived 45%
mapped network: 1948 flowlines, 805.7 km

NHDPlus HR ftype codes are 460 for stream/river, 558 for artificial path, 334 for connector and 428 for pipeline.

Example 3 — save the chosen network with its threshold 

chosen = 0.1
out = networks[chosen].assign(threshold_km2=chosen, length_m=lambda d: d.length.round(1))
out.to_file("esopus_streams.gpkg", layer=f"streams_{chosen}km2")
print(out.describe().loc[["count", "mean", "max"], ["length_m"]])
          length_m
count  2393.000000
mean    361.476766
max    2588.500000

Explanation

Why there is no natural threshold

In the landscape, channels start where concentrated flow erodes a defined bed, which depends on slope, soil, vegetation and climate. In a DEM, every cell has some accumulation and the transition from hillslope to channel is gradual. The threshold is a model of channel initiation, calibrated by comparison, not a property that can be read off the grid.

Why precision and recall pull apart

A low threshold creates channels up every hollow, including many that are not mapped streams, so precision falls. A high threshold keeps only large rivers, which are all mapped, but omits the small tributaries that make up most of any network's length, so recall falls. First-order streams were 443.6 km of NHDPlus HR's 805.7 km in the basin; they are exactly what high thresholds drop.

Why Strahler order depends on the threshold

Strahler order increases by one only where two streams of equal order meet. Adding small first-order streams at a low threshold creates more junctions of equal order all the way down the network, so the outlet's order rises. The Esopus at Coldbrook was order 7 at 0.05 km² and order 4 at 10 km²; NHDPlus HR gives it order 5. Orders from different thresholds or sources are not comparable; see stream order explained.

Why derived channels drift from mapped ones

A 10 m DEM places the channel in the lowest cells, which on a wide valley floor may be tens of metres from the mapped centreline, and conditioning moves channels through breaches. The 50 m tolerance absorbs most of this; where it does not, burning the mapped network into the DEM forces agreement.

Bar chart of the maximum Strahler order in the Esopus basin for accumulation thresholds from 0.05 to 10 km2.
The river did not change; the count of small tributaries did, and the order followed.

Edge cases or notes

  • Flat valley floors produce straight parallel channels; see fixing streams on flats.
  • Roads and bridges can end channels abruptly; see streams that stop at roads.
  • Coarser DEMs need larger thresholds and miss narrow valleys; see how DEM resolution changes a drainage network.
  • Variable thresholds based on slope and area (such as an area–slope product) represent channel heads better in mixed terrain.
  • Segment counts depend on how the tool splits lines at junctions; compare lengths rather than counts across tools.
  • Lakes and reservoirs appear as straight lines across flats; clip or flag them.
  • Clip before extracting to avoid computing networks outside the catchment and paying for them.

FAQ

What threshold should I use to extract streams from a DEM?

One calibrated against a mapped network for your area. On the 10 m Esopus DEM, 0.1 km² matched NHDPlus High Resolution best, with 85% of derived and 89% of mapped length agreeing within 50 m.

Why does my DEM stream network have too many streams?

The threshold is too low for the reference you compare with. At 0.05 km² the Esopus network was 1,260 km long, against 806 km mapped.

How do I extract a stream network with pysheds?

Compute flow direction and accumulation, then call grid.extract_river_network(fdir, acc > threshold_cells), which returns GeoJSON line features you can load into GeoPandas.

Should the threshold be in cells or in area?

In area. The same number of cells is nine times the area on a 30 m grid as on a 10 m grid; convert an area threshold to cells for each DEM.

Why does the Strahler order of my river change?

Because it depends on how many small streams the threshold creates. The Esopus outlet was order 7 at 0.05 km² and order 4 at 10 km².

How do I check a derived stream network?

Buffer each network and measure how much of the derived length lies near mapped streams and how much of the mapped length lies near derived ones. Both numbers are needed.