11.5 Watersheds and Stream Networks
Modelling water flow on a DEM — filling, flow direction, accumulation, and watershed delineation.
Key takeaways
- Hydrology on a DEM follows a standard pipeline: fill depressions → flow direction → flow accumulation → streams + watersheds.
- D8 and D-infinity are the two common flow algorithms.
- Watershed delineation is the foundation of most hydrological GIS.
Introduction
Water flows downhill. A DEM encodes enough to model where. This lesson walks through the classic hydrology pipeline — a workflow that underlies flood modelling, erosion, pollution tracking, and countless environmental applications.
The five-step pipeline
- Fill depressions — remove pits that would trap water.
- Compute flow direction — which neighbour each cell drains to.
- Compute flow accumulation — upstream contributing area per cell.
- Extract stream network — cells above a flow-accumulation threshold.
- Delineate watersheds — group cells by common outlet.
Fill depressions
A pit is a cell lower than all its neighbours. In reality these are often real features (lakes), but for flow modelling they trap water. Standard practice: raise pit cells to the lowest neighbour's elevation before computing flow.
1from pysheds.grid import Grid
2grid = Grid.from_raster('dem.tif')
3dem = grid.read_raster('dem.tif')
4flooded = grid.fill_depressions(dem)
5inflated = grid.resolve_flats(flooded)Or with GRASS: r.fill.dir.
Flow direction
For each cell, which way does water flow?
D8
The simplest: water goes to the lowest of the 8 neighbours. Result: one of 8 discrete directions.
Advantages: simple, fast. Disadvantages: forces water into one neighbour; real flow can fan out.
D-infinity
Water is apportioned between the two downslope neighbours that form a triangle. Output is a continuous direction angle.
Advantages: more realistic dispersion on hillslopes. Disadvantages: slightly more complex; some tools still D8-only.
fdir = grid.flowdir(inflated, routing='d8') # or 'dinf'Flow accumulation
How many upstream cells flow into each cell? Accumulation is calculated by summing contributions following flow direction.
facc = grid.accumulation(fdir, routing='d8')Cells with high accumulation are channels; low-accumulation cells are on hillslopes.
Extract streams
Threshold the accumulation raster:
streams = (facc > 1000).astype('uint8')Choose the threshold based on drainage density and the smallest channel you care to capture. For a 30 m DEM, 1 000 cells = ~1 km² contributing area — a modest stream. For fine rasters, thousands of cells.
Vectorise to a stream network:
gdal_polygonize.py streams.tif streams.gpkgOr use PySheds extract_river_network for proper stream ordering.
Stream ordering
Assign each stream segment a Strahler order:
- Order 1 — headwaters.
- Order 2 — junction of two order-1.
- Order 3 — junction of two order-2.
- ... up to major rivers.
Useful for classifying rivers by size in cartography or analysis.
Delineate watersheds
Given an outlet point, the watershed is the set of cells that drain to it:
1x_outlet, y_outlet = -73.96, 40.78
2x_snap, y_snap = grid.snap_to_mask(facc > 100, (x_outlet, y_outlet))
3watershed = grid.catchment(x=x_snap, y=y_snap,
4 fdir=fdir, dirmap=dirmap,
5 xytype='coordinate')Repeat for multiple outlets (a gauge network) to produce all sub-watersheds.
Fill / flow / accumulation in GDAL
GRASS r.watershed wraps the whole pipeline. QGIS has GRASS tools built in.
Limitations and realism
- Flat areas — algorithms must resolve which direction water flows; imposed gradients introduce artefacts.
- Real hydrology — streams are not just topographic channels; culverts, drainage tiles, and engineered structures reroute water.
- Urban environments — DEMs from LiDAR include buildings and sidewalks; hydrologic models need terrain-only DEMs with storm-drain networks separately.
- Scale — 30 m DEMs miss small channels and drainage density.
Modern alternatives
- SCALGO Live — commercial flood modelling with DEM processing at continental scale.
- MERIT-HYDRO — global 3-arcsec hydrography with stream widths.
- HydroSHEDS — global drainage network derived from SRTM.
- 2D hydrodynamic models (HEC-RAS, LISFLOOD-FP) — for actual flood inundation.
Worked example — watershed for a gauge
Given a USGS stream gauge at (−83.12, 36.48):
1import pysheds.grid as pg
2grid = pg.Grid.from_raster('dem.tif')
3dem = grid.read_raster('dem.tif')
4flooded = grid.fill_depressions(dem)
5inflated = grid.resolve_flats(flooded)
6fdir = grid.flowdir(inflated)
7facc = grid.accumulation(fdir)Export catch as a raster or vectorise for display.
Self-check exercises
1. Why fill depressions before computing flow direction?
A pit — a cell lower than all its neighbours — would trap water. Water in reality often passes through such cells (buffered by ponds, soaked into soil, flowed over during storms). Filling raises pits to the lowest neighbour's elevation, ensuring continuous flow paths from headwaters to outlets. Without filling, flow accumulation stops at the pits.
2. D8 vs D-infinity — when does it matter?
D8 forces all water into one neighbour; D-infinity apportions between two. On hillslopes where water really does diverge, D-infinity is more realistic. On channels (high accumulation cells), D8 is usually fine since most water is in a well-defined line. D-infinity matters most for small-scale hillslope processes (erosion, hillslope hydrology); D8 is adequate for watershed-scale modelling.
3. Your extracted stream network includes straight lines crossing a road. What's happening?
The road embankment acts as a topographic dam; the DEM thinks water flows over it. In reality a culvert carries water underneath. Fix by explicitly lowering DEM cells at known culvert locations (a "stream burning" technique), or by representing culverts as vector features and integrating them into a hydrologic model that supports them.
Summary
- Pipeline: fill → flow direction → accumulation → streams → watersheds.
- D8 is simple; D-infinity is more realistic on hillslopes.
- Stream thresholds and ordering are cartographic and analytical choices.
- DEM hydrology is a starting point; real modelling adds culverts, soils, and dynamics.
Further reading
- Wilson & Gallant — Terrain Analysis: Principles and Applications.
- pysheds documentation.
- HydroSHEDS project materials.
- Lehner, B. & Grill, G. — Global river hydrography and network routing.
