Extracting drainage basins from digital elevation models
2026-02-27
Water quality managers in the Athabasca River basin need to answer a question that seems simple but requires a surprising amount of spatial computation: which parcels of land — which farms, which industrial sites, which municipalities — drain into this reach of river? The answer determines whose effluent license applies, whose operations require monitoring, and whose land practices affect the downstream fishery. Before digital elevation models and flow-routing algorithms, watershed boundaries were traced by hand on topographic maps by experienced technicians. Today a GIS can delineate hundreds of nested watersheds in a large basin automatically from a DEM, consistently and reproducibly, in minutes.
Watershed delineation is the automatic extraction of drainage basin boundaries from terrain data. The key insight is that every cell in a DEM can be assigned a flow direction to the steepest downslope neighbour; a complete map of flow directions is a directed graph; and the watershed of any outlet point is simply all the cells that eventually drain to it by following the graph. This is the same logic as the hydrological flow model in Series 1, now applied at the scale of complete basins and implemented with algorithms (D8 flow direction, breach filling of sinks, flow accumulation) that handle the practical complications of real terrain data. This model derives those algorithms, implements the watershed extraction procedure, and discusses the sensitivity of results to DEM resolution and sink-filling choices.
If it rains on this mountain, where does the water eventually flow?
Watershed (also called catchment or drainage basin): Area that drains to a common outlet.
Applications: - Flood prediction: How much area contributes flow to this point? - Water quality: Which land uses drain into this stream? - Reservoir planning: What’s the contributing area to a dam site? - Erosion modelling: Where does sediment originate and accumulate? - Ecological assessment: Upstream impacts on downstream habitats
The mathematical question: Given a DEM, how do we automatically extract: 1. Flow direction from each cell 2. Flow accumulation (contributing area) 3. Stream network (where flow accumulates) 4. Watershed boundaries for any pour point
Water flows downhill along path of steepest descent.
For each cell: Determine which neighbor receives flow.
D8 method (8-direction): - Examine all 8 neighbors - Flow to neighbor with steepest descent - Ties: Use deterministic rule (e.g., choose first in priority order)
Flow direction encoding:
32 64 128
16 [X] 1
8 4 2Powers of 2 → single value encodes direction.
Example: Flow to east (right) → direction = 1
For each cell: Count how many upstream cells drain through it.
Algorithm: 1. Start with cells having no upstream contributors (ridges) 2. Process in topological order (high to low elevation) 3. Each cell passes its accumulated flow to its downstream neighbor
Result: Raster where value = number of contributing cells.
High accumulation = stream channels (many cells drain here)
Threshold accumulation to identify streams:
\text{stream} = (\text{accumulation} > T)
Typical threshold: T = 100-1000 cells (depends on DEM resolution)
Result: Binary raster of stream network.
Pour point: Outlet location where we want to define watershed.
Watershed: All cells that drain to (flow through) pour point.
Algorithm: Trace upslope from pour point, marking all contributing cells.
For cell at (i, j) with elevation z_{i,j}:
Compute slope to each neighbor:
s_n = \frac{z_{i,j} - z_n}{d_n}
Where: - z_n = elevation of neighbor - d_n = distance to neighbor (1 for cardinal, \sqrt{2} for diagonal)
Flow direction:
\text{dir}_{i,j} = \arg\max_n(s_n)
(Direction of maximum slope)
Edge cases:
Flat area (s_n = 0 for all neighbors): - Mark as undefined or use flat resolution algorithm
Depression (all neighbors higher): - Mark as sink or fill depression
Depressions (local minima) disrupt flow routing.
Fill algorithm:
1. Initialize output DEM = input DEM
2. Set boundary cells to input elevation
3. Set interior cells to very high value
4. Iterate:
- For each interior cell:
- New elevation = min(current, max(neighbors) + epsilon)
5. Repeat until convergenceResult: DEM with depressions filled to spill elevation.
Epsilon (small value like 0.01m) ensures flow out of filled areas.
Pseudocode:
function flow_accumulation(flow_direction):
# Initialize
accumulation = ones(rows, cols) # Each cell contributes itself
# Sort cells by elevation (high to low)
cells = sort_by_elevation_descending(DEM)
# Process in topological order
for each cell in cells:
downstream = get_downstream_cell(cell, flow_direction)
if downstream exists:
accumulation[downstream] += accumulation[cell]
return accumulationComplexity: O(n) after sorting O(n \log n)
Key: Topological ordering ensures upstream cells processed before downstream.
Upslope area from pour point:
function extract_watershed(flow_direction, pour_point):
watershed = empty_set()
to_process = {pour_point}
while to_process not empty:
cell = to_process.pop()
watershed.add(cell)
# Find all neighbors that flow into this cell
for each neighbor of cell:
if flow_direction[neighbor] points to cell:
to_process.add(neighbor)
return watershedResult: Set of all cells contributing flow to pour point.
Problem: Compute flow direction and accumulation for this DEM.
Elevation (meters):
j=0 j=1 j=2 j=3
i=0 50 48 46 44
i=1 52 45 43 42
i=2 54 47 44 40
i=3 56 49 46 38Use D8 method (8 neighbors).
Step 1: Flow Direction
Cell (0,1) elevation 48:
Neighbors and slopes: - (0,0): (48-50)/1 = -2 (uphill, ignore) - (0,2): (48-46)/1 = 2 ✓ - (1,0): (48-52)/√2 = -2.83 (uphill) - (1,1): (48-45)/1 = 3 ✓ - (1,2): (48-43)/√2 = 3.54 ✓ Maximum
Flow direction: Southeast diagonal (code = 2)
Cell (1,1) elevation 45:
Neighbors: - (0,0): -7/√2 (uphill) - (0,1): -3 (uphill) - (0,2): 2 - (1,0): -7 (uphill) - (1,2): 2 - (2,0): -9/√2 (uphill) - (2,1): -2 (uphill) - (2,2): 1
Maximum slope: (1,2) or (0,2) both = 2. Choose (1,2) by convention.
Flow direction: East (code = 1)
Continue for all cells…
Flow direction map:
j=0 j=1 j=2 j=3
i=0 1 2 2 4
i=1 1 1 2 4
i=2 1 1 2 4
i=3 1 1 1 4(1=E, 2=SE, 4=S)
Step 2: Flow Accumulation
Process cells in descending elevation order:
(3,0) elevation 56: Accumulation = 1, flows to (3,1) - Update (3,1): 1 + 1 = 2
(2,0) elevation 54: Accumulation = 1, flows to (2,1) - Update (2,1): 1 + 1 = 2
(1,0) elevation 52: Accumulation = 1, flows to (1,1) - Update (1,1): 1 + 1 = 2
(0,0) elevation 50: Accumulation = 1, flows to (0,1) - Update (0,1): 1 + 1 = 2
(3,1) elevation 49: Accumulation = 2, flows to (3,2) - Update (3,2): 1 + 2 = 3
Continue…
Final accumulation:
j=0 j=1 j=2 j=3
i=0 1 2 4 8
i=1 1 2 5 10
i=2 1 2 6 13
i=3 1 2 3 16Cell (3,3): Accumulates flow from all 16 cells (entire grid drains here).
This is the watershed outlet.
Below is an interactive watershed delineation tool.
<label>
Terrain type:
<select id="terrain-type">
<option value="valley" selected>Simple Valley</option>
<option value="ridge">Ridge and Valleys</option>
<option value="complex">Complex Terrain</option>
</select>
</label>
<label>
Stream threshold:
<input type="range" id="stream-threshold" min="10" max="200" step="10" value="50">
<span id="stream-threshold-val">50</span> cells
</label>
<label>
Show flow direction:
<input type="checkbox" id="show-flow-dir">
</label>
<label>
Show streams:
<input type="checkbox" id="show-streams" checked>
</label>
<div class="watershed-info">
<p><strong>Click any cell to extract its watershed</strong></p>
<p><strong>Watershed area:</strong> <span id="watershed-area">--</span> cells</p>
<p><strong>Max accumulation:</strong> <span id="max-accum">--</span> cells</p>
</div><canvas id="watershed-canvas" width="700" height="350" style="border: 1px solid #ddd; cursor: crosshair;"></canvas>Try this: - Click right panel to set pour point (red dot) and extract watershed (green overlay) - Simple valley: Water flows down center channel - Ridge and valleys: Divides drainage into separate basins - Complex terrain: Multiple drainage patterns - Stream threshold: Higher = only major streams shown (blue) - Show flow direction: Red arrows show where each cell drains - Left: Elevation (lighter = higher) - Right: Flow accumulation (darker blue = more contributing area)
Key insight: Watershed boundaries follow ridges—all rain inside the green area flows to the red pour point!
Contributing area = watershed size upstream of point.
Larger watershed → more runoff → higher flood potential
Example: Dam site with 500 km² watershed - Expected peak flow during storm - Reservoir sizing for flood control
Pollutant sources within watershed affect downstream quality.
Non-point source pollution: - Agricultural runoff (fertilizers, pesticides) - Urban stormwater (oil, sediment) - Septic systems
Delineate watershed → identify all potential sources → prioritize management.
Strahler stream order:
Higher order = larger river
Application: Stream classification for habitat assessment
Runoff coefficient (C):
Q = C \times P \times A
Where: - Q = runoff volume - P = precipitation - A = watershed area - C = 0.1 (forest) to 0.9 (pavement)
Reservoir inflow = accumulated runoff from watershed
D8 fails where no downslope neighbor exists.
Solutions:
1. Impose flow direction: - Route to nearest definite flow (edge of flat)
2. D-infinity algorithm: - Allow fractional flow to multiple neighbors - More realistic but more complex
3. Epsilon slope: - Add tiny random elevation noise to break ties
Natural depressions: - Lakes (real features) - DEM artifacts (errors)
Consequences: - Flow accumulation terminates - Watershed extraction incomplete
Solutions:
Fill depressions: - Raise elevation to spill point - Loses lake information
Breach depressions: - Cut channel through lowest barrier - Preserves topography better
Hybrid: - Fill small artifacts, preserve large lakes
Coarse DEM misses small drainage features: - Ditches - Culverts - Small channels
Result: Incorrectly routed flow
Solution: Use finest available DEM (1m LiDAR ideal)
Watershed extends beyond DEM boundary:
Flow direction undefined at edges.
Solution: - Use larger DEM extent - Or explicitly model boundary conditions
D8 limitation: Flow restricted to 8 directions (multiples of 45°).
D-infinity: Flow in any direction.
Method:
Flow partitioning:
If steepest descent at angle \theta between neighbors A and B:
f_A = \frac{\theta_B - \theta}{\theta_B - \theta_A}
f_B = 1 - f_A
Result: More realistic flow divergence
Cost: Higher computational complexity
Topological sort: Linear ordering of directed acyclic graph (DAG) such that for every edge u \to v, u comes before v in the ordering.
Flow network = DAG where edges point downstream.
Topological order = process cells from high to low elevation.
Ensures: When processing cell, all upstream contributors already processed.
L = empty list (topological order)
S = set of nodes with no incoming edges
while S not empty:
remove node n from S
add n to L
for each node m with edge n → m:
remove edge n → m
if m has no other incoming edges:
add m to SFor flow: Start with ridge cells (no upstream), process downhill.
Simpler for DEMs: Sort by elevation instead of graph structure.
Works because: Higher elevation → upstream in flow network.