🧭 Concept Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 1 of 20 Updated

DEM Hydrology Explained: From Elevation to Where Water Goes

Problem statement

A digital elevation model says how high every cell is. It does not say where rain goes — that has to be derived, in a chain of steps where each one depends on the last: remove the depressions that would trap water, give every cell a flow direction, count how many cells drain through each one, pick an outlet on the resulting channel, and trace everything upstream of it. Skip or mishandle one step and the chain still produces a map, only of the wrong watershed.

Measured on the Esopus Creek catchment above the USGS gauge at Coldbrook, New York — published drainage area 497.3 km² — with the USGS 3DEP 10 m elevation model:

  • On the raw DEM, no cell drained more than 11.95 km². 5,768 pits and 45,736 flat cells stopped flow long before it reached the gauge.
  • Only 0.96% of the cells sat in depressions, and the median depression was 0.15 m deep. Conditioning changed less than 1% of the DEM and made the whole river connect.
  • The conditioned DEM gave a 492.4 km² watershed, 1.0% under the published area, overlapping the USGS reference polygon with an intersection-over-union of 0.985.
  • A gauge coordinate used as the outlet without snapping drained one cell. Moved 160 m onto the channel, it drained the whole basin.

Quick answer

The standard pipeline, one step per line, with WhiteboxTools:

import whitebox

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.fill_depressions("dem.tif", "filled.tif", fix_flats=True)        # 1. condition
wbt.d8_pointer("filled.tif", "d8.tif")                                # 2. flow direction
wbt.d8_flow_accumulation("d8.tif", "acc.tif", out_type="cells", pntr=True)   # 3. flow accumulation
wbt.snap_pour_points("gauge.shp", "acc.tif", "outlet.shp", snap_dist=150)    # 4. outlet on the channel
wbt.watershed("d8.tif", "outlet.shp", "watershed.tif")                # 5. everything upstream

Use a projected DEM with metre cells, condition it before anything else, and check the watershed area against a published value before trusting it.

Flow from a raw DEM through depression conditioning, D8 flow direction, flow accumulation and outlet snapping to a delineated watershed, with the measured result of each step for Esopus Creek.
Five steps, each consuming the previous one's output; the watershed is only as good as the weakest.

Step-by-step solution

1. Start from a projected, bare-earth DEM

Flow routing compares neighbouring elevations and measures distances and areas in cell units, so cells should be square and in metres. The 3DEP data were reprojected from geographic coordinates to UTM zone 18N at 10 m. A bare-earth model (a DTM) is the right source: surface models include trees and buildings that act as dams. Resolution matters too: the same workflow on the 30 m Copernicus GLO-30 DSM produced a 450.7 km² watershed, 9.4% short — see how DEM resolution changes a drainage network.

2. Remove depressions

Real landscapes have few closed hollows; elevation models have many, from measurement noise, bridges, culverts and interpolation. In the Esopus basin, 80,057 cells — 0.96% — were in depressions, most of them shallow: the median depth was 0.148 m, the deepest 10.13 m. Every depression is a local sink where flow stops. Filling raises each one to its spill point; breaching cuts a channel through its rim. Filling changed 80,209 cells; least-cost breaching changed 52,499. Both gave watersheds within 0.4 km² of each other — see breaching or filling a DEM.

3. Assign a flow direction to every cell

D8, the most common method, sends each cell's flow to the steepest of its eight neighbours. On the raw DEM, pysheds could not assign a direction to 45,736 flat cells and 5,768 pits. On the conditioned DEM those numbers fell to 3 flats and 194 pits, all at the edge of the data. D8 is simple and makes a clean single-thread network; multi-direction methods spread flow over hillslopes more realistically — see flow direction explained.

4. Accumulate flow

Flow accumulation counts, for every cell, how many cells drain through it. Multiplied by the cell area it is the upstream drainage area. On the raw DEM the maximum was 11.95 km², because every pit ended a flow path; on the conditioned DEM it was 513.5 km² at the edge of the data and 492.4 km² at the gauge. The values span six orders of magnitude: half of all cells received flow from 10 cells or fewer.

5. Put the outlet on the channel

A gauge's recorded coordinates are rarely exactly on the one-cell-wide channel the DEM produces. At Coldbrook, the gauge's cell had an accumulation of 1 — it drained only itself. The nearest cell with more than 1 km² of drainage was 40 m away. Snapping the outlet to the highest accumulation within 150 m moved it onto Esopus Creek and gave the full watershed — see pour points explained.

6. Trace the watershed

Walking the flow directions upstream from the outlet marks every cell that drains to it. The 3DEP watershed held 4,924,422 cells: 492.44 km², 0.97% under the published 497.28 km², with 7.31 km² of symmetric difference against the USGS NLDI basin polygon. Converted to a polygon, its geodesic area was 492.80 km².

7. Derive streams from accumulation

A stream network is the set of cells whose accumulation exceeds a threshold. A 1 km² threshold gave 257 segments and 332.9 km of streams; 0.1 km² gave 2,397 segments and 865.6 km. The threshold is a choice, not a property of the DEM — see extracting a stream network.

8. Validate against something independent

A published drainage area, a mapped basin boundary, or a reference stream network each catch different errors. A watershed that looks plausible can still be 74% too small, as one pysheds conditioning sequence produced here — see fixing a watershed that comes out as a few pixels.

Bar chart of the largest drainage area reachable on the raw Esopus DEM against the conditioned DEM, and the watershed area at the gauge against the published value.
Changing under 1% of the cells turned 12 km² of disconnected drainage into a 492 km² river.

Code examples

Example 1 — delineate a watershed from a DEM and a gauge coordinate 

import os

import geopandas as gpd
import rasterio
import whitebox
from pyproj import Transformer
from shapely.geometry import Point


def delineate(dem, lon, lat, workdir="work", snap_m=150):
    """Fill, D8, accumulate, snap and trace; return the watershed area in km2."""
    os.makedirs(workdir, exist_ok=True)
    wbt = whitebox.WhiteboxTools()
    wbt.set_verbose_mode(False)
    wbt.set_working_dir(os.path.abspath(workdir))
    with rasterio.open(dem) as src:
        crs, cell = src.crs, src.res[0]
    x, y = Transformer.from_crs("EPSG:4269", crs, always_xy=True).transform(lon, lat)
    gpd.GeoDataFrame(geometry=[Point(x, y)], crs=crs).to_file(os.path.join(workdir, "gauge.shp"))
    wbt.fill_depressions(os.path.abspath(dem), "filled.tif", fix_flats=True)
    wbt.d8_pointer("filled.tif", "d8.tif")
    wbt.d8_flow_accumulation("d8.tif", "acc.tif", out_type="cells", pntr=True)
    wbt.snap_pour_points("gauge.shp", "acc.tif", "outlet.shp", snap_dist=snap_m)
    wbt.watershed("d8.tif", "outlet.shp", "watershed.tif")
    with rasterio.open(os.path.join(workdir, "watershed.tif")) as src:
        cells = int((src.read(1) == 1).sum())
    return cells * cell * cell / 1e6


area = delineate("esopus_3dep13_utm_basin.tif", -74.2701944, 42.0144722)
print(f"watershed {area:.2f} km2, {100 * (area / 497.28 - 1):+.2f}% against the published 497.28 km2")
watershed 492.42 km2, -0.98% against the published 497.28 km2

The gauge coordinates are USGS NAD83 values, so the transformer starts from EPSG:4269.

Example 2 — what conditioning changes 

import numpy as np

if not hasattr(np, "in1d"):          # pysheds 0.5 calls np.in1d, removed in NumPy 2.4
    np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
from pysheds.grid import Grid


def routing_summary(path):
    grid = Grid.from_raster(path)
    dem = grid.read_raster(path)
    fdir = grid.flowdir(dem)
    acc = grid.accumulation(fdir)
    cell_km2 = abs(grid.affine.a * grid.affine.e) / 1e6
    print(f"{os.path.basename(path)}: flats {int((fdir == -1).sum()):,}, pits {int((fdir == -2).sum()):,}, "
          f"largest drainage area {float(acc.max()) * cell_km2:.2f} km2")


routing_summary("esopus_3dep13_utm_basin.tif")
routing_summary("work/filled.tif")
esopus_3dep13_utm_basin.tif: flats 45,736, pits 5,768, largest drainage area 11.95 km2
filled.tif: flats 3, pits 194, largest drainage area 513.49 km2

pysheds marks cells with no downslope neighbour as −1 (flat) or −2 (pit) in its flow-direction grid, which makes the effect of conditioning easy to count.

Example 3 — how much of the DEM sits in depressions 

def depression_report(dem, workdir="work"):
    wbt = whitebox.WhiteboxTools()
    wbt.set_verbose_mode(False)
    wbt.set_working_dir(os.path.abspath(workdir))
    wbt.depth_in_sink(os.path.abspath(dem), "sink_depth.tif", zero_background=True)
    with rasterio.open(dem) as src:
        valid = src.read(1) != src.nodata
    with rasterio.open(os.path.join(workdir, "sink_depth.tif")) as src:
        depth = src.read(1)
    sinks = valid & (depth > 0)
    print(f"cells in depressions {sinks.sum():,} ({sinks.sum() / valid.sum():.2%}); "
          f"median depth {np.median(depth[sinks]):.3f} m, deepest {depth[valid].max():.2f} m")


depression_report("esopus_3dep13_utm_basin.tif")
cells in depressions 80,057 (0.96%); median depth 0.148 m, deepest 10.13 m

Explanation

Why raw elevation cannot route water

Flow routing needs every cell to have a lower neighbour until the flow reaches the edge of the data or a real sink. A DEM built from lidar or radar has millions of tiny irregularities, and any cell lower than all eight neighbours ends every flow path that reaches it. With 5,768 such pits in the basin, the longest connected flow path drained 11.95 km² — a stream the size of a small tributary.

Why a small change has a large effect

Depressions are small and shallow but they sit exactly where water collects: in valley bottoms, behind road embankments and bridges, along the channel. Raising 1% of cells by a median of 15 cm reconnected the valley floor, and because accumulation is cumulative, reconnecting a single cell on the main stem adds everything upstream of it.

Why the outlet must be snapped

A derived channel is one cell wide, and the gauge is a point measured independently of the DEM. Unless it happens to fall in exactly the right 10 m cell, its accumulation belongs to a hillslope beside the river. Snapping to high accumulation nearby corrects that, at the risk of jumping to a larger neighbouring river if the search radius is too generous.

Why the answer still differs from the published area

The published drainage area comes from a different method and source. The DEM channel may leave the true river at braided reaches, reservoirs or culverts, and a 10 m grid cannot represent the drainage divide more precisely than a cell. A 1% difference is normal; 10% or more points at the DEM source, the conditioning or the outlet.

Bar chart of flow accumulation percentiles inside the Esopus basin, from a median of 10 cells to 35 km2 at the 99.9th percentile, on a logarithmic scale.
Accumulation is almost all hillslope: half the cells drain 10 cells or fewer, a tiny fraction carry the river.

Edge cases or notes

  • pysheds 0.5 and NumPy 2.4+ fail with AttributeError: module 'numpy' has no attribute 'in1d'; the one-line shim in Example 2 restores it.
  • Geographic DEMs give wrong areas unless each cell's area is computed from its latitude; see catchment areas on a latitude–longitude DEM.
  • Reservoirs and lakes are large flat areas; flow crosses them along arbitrary paths after flat resolution.
  • Karst and closed basins have real depressions that should not be filled.
  • Roads and bridges create artificial dams in the DEM; see streams that stop at roads.
  • Edges of the DEM cut off upstream area if the data do not cover the whole catchment.
  • Surface models (DSMs) include vegetation and buildings; prefer a bare-earth DTM.

FAQ

What is DEM hydrology?

Deriving where water flows from an elevation model: conditioning the DEM, computing flow directions and accumulation, and using them to delineate watersheds and stream networks.

Why do I need to fill or breach a DEM before flow routing?

Depressions stop flow. On the raw Esopus DEM no cell drained more than 11.95 km²; after conditioning under 1% of cells, the watershed at the gauge was 492.4 km².

How accurate is a DEM-derived watershed?

With a 10 m bare-earth DEM and a snapped outlet, the Esopus watershed was 1.0% under the published area. The 30 m GLO-30 surface model gave a watershed 9.4% too small.

Why does my watershed contain only one cell?

The outlet is not on the derived channel. The Coldbrook gauge's own cell drained nothing else; snapping it 160 m to the highest accumulation nearby gave the whole basin.

Which Python libraries do DEM hydrology?

WhiteboxTools and pysheds are the common choices; both condition DEMs, compute D8 flow direction and accumulation and delineate watersheds. RichDEM and TauDEM are alternatives.

What cell size should the DEM have?

As fine as the landscape needs and in metres. A 10 m DEM resolved the Esopus basin to within 1%; at 30 m, narrow valleys and divides were lost.