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 13.03 km². 5,781 pits and 46,539 flat cells stopped flow long before it reached the gauge.
• Only 1.00% of the cells sat in depressions, and the median depression was 1.27 m deep. Conditioning changed about 1% of the DEM and made the whole river connect.
• The conditioned DEM gave a 492.2 km² watershed, 1.0% under the published area, overlapping the USGS reference polygon with an intersection-over-union of 0.985.
• How far you clip the DEM decided the answer. With the DEM cut 1 km around the basin, filling gave 492.2 km². Cut 1.5 km around it, the edge ran round the Ashokan Reservoir instead of through it, filling flooded the valley 61.7 m deep at the gauge, and the watershed came out 448.3 km².
• A gauge coordinate used as the outlet without snapping drained one cell. Moved 78 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_max_procs(1)                      # one thread: parallel runs give slightly different results
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.

Bare earth, no depressions, flow direction, accumulation, outlet

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 448.2 km² watershed, 9.9% 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, 83,652 cells — 1.00% — were in depressions: the median depth was 1.27 m, the deepest 20.47 m. Filling raises each one to its spill point; breaching cuts a channel through its rim. Filling changed 84,086 cells; least-cost breaching with a 20-cell search changed 70,263. Both gave watersheds within 0.3 km² of each other — on this DEM, clipped this way — 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 46,539 flat cells and 5,781 pits. On the conditioned DEM those numbers fell to no flats and 195 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 the cells draining through each cell; times the cell area, it is the upstream drainage area. On the raw DEM the maximum was 13.03 km², because every pit ended a flow path; on the conditioned DEM it was 513.4 km² at the edge of the data and 492.2 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,922,072 cells: 492.21 km², 1.02% under the published 497.28 km², with 7.25 km² of symmetric difference against the USGS NLDI basin polygon. As a polygon its geodesic area was 492.56 km² (measured here).

7. Derive streams from accumulation

A stream network is the set of cells whose accumulation exceeds a threshold. In extracting a stream network, a 1 km² threshold gave 263 segments and 331.2 km of streams; 0.1 km² gave 2,360 segments and 865.1 km. The threshold is a choice, not a property of the DEM.

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 8% too small: one pass of pysheds' fill-and-resolve sequence gave 459.0 km² here, and only a second pass reached 491.0 km² — see fixing a watershed that comes out as a few pixels.

From DEM and gauge to a delineated watershed

The data used on this page

The setup builds every input this page reads: the basin of Esopus Creek above Coldbrook, New York (USGS 01362500), from the USGS NLDI; USGS 3DEP 1/3 arc-second elevation, resampled bilinearly to UTM 18N at 10 m over the basin plus 1 km (esopus_3dep13_utm_basin.tif).

import os
import time
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import requests
from rasterio.merge import merge
from rasterio.warp import Resampling, reproject, transform_bounds


def get(url, **kw):
    """GET with patient retries: the USGS services answer 504 when busy, and recover."""
    for attempt in range(6):
        try:
            r = requests.get(url, timeout=600, **kw)
            r.raise_for_status()
            return r
        except requests.RequestException:
            if attempt == 5:
                raise
            time.sleep(30 * (attempt + 1))


# basins from the USGS NLDI
for site in ('01362500',):
    if not os.path.exists(f"nldi_basin_{site}.geojson"):
        Path(f"nldi_basin_{site}.geojson").write_text(get(f"https://api.water.usgs.gov/nldi/linked-data/nwissite/USGS-{site}/basin?f=json").text)
basin = gpd.read_file("nldi_basin_01362500.geojson")     # Esopus Creek at Coldbrook
W, S, E, N = basin.total_bounds + np.array([-0.03, -0.03, 0.03, 0.03])

TILES = {"3dep13": ([f"/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/{t}/USGS_13_{t}.tif"
                      for t in ("n43w075", "n42w075")], 10),
         "glo30": ([f"/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_{t}_00_W075_00_DEM/"
                    f"Copernicus_DSM_COG_10_{t}_00_W075_00_DEM.tif" for t in ("N42", "N41")], 30)}


def dem_box(name, buffer_m, out):
    """UTM 18N grid of whole cells around the basin plus buffer_m, bilinear from the source tiles, no nodata."""
    if os.path.exists(out):
        return
    paths, res = TILES[name]
    srcs = [rasterio.open(p) for p in paths]
    mosaic, tf = merge(srcs, bounds=(W - 0.1, S - 0.1, E + 0.1, N + 0.1))
    bx = basin.to_crs(32618).total_bounds
    x0, y0 = np.floor((bx[0] - buffer_m) / res) * res, np.floor((bx[1] - buffer_m) / res) * res
    x1, y1 = np.ceil((bx[2] + buffer_m) / res) * res, np.ceil((bx[3] + buffer_m) / res) * res
    w, h = int((x1 - x0) / res), int((y1 - y0) / res)
    grid_tf = rasterio.transform.from_origin(x0, y1, res, res)
    grid = np.full((h, w), -9999, "float32")
    reproject(mosaic[0], grid, src_transform=tf, src_crs=srcs[0].crs, dst_transform=grid_tf, dst_crs="EPSG:32618",
              resampling=Resampling.bilinear, src_nodata=srcs[0].nodata, dst_nodata=-9999)
    with rasterio.open(out, "w", driver="GTiff", height=h, width=w, count=1, dtype="float32", crs="EPSG:32618",
                       transform=grid_tf, nodata=-9999, compress="deflate") as dst:
        dst.write(grid, 1)


dem_box("3dep13", 1000, "esopus_3dep13_utm_basin.tif")   # 3DEP 1/3 arc-second at 10 m, basin + 1 km

for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm_basin.tif']:
    if f.endswith(".tif"):
        with rasterio.open(f) as r:
            print(f"{f}: {r.width:,} x {r.height:,} cells of {r.res[0]:.6g} {'m' if r.crs.is_projected else 'degrees'}")
    elif f.endswith(".zip"):
        print(f"{f}: {os.path.getsize(f) / 1e6:.1f} MB")
    else:
        print(f"{f}: {len(gpd.read_file(f)):,} features")

nldi_basin_01362500.geojson: 1 features
esopus_3dep13_utm_basin.tif: 3,471 x 2,408 cells of 10 m

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_max_procs(1)                      # one thread: parallel runs give slightly different results
    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.21 km2, -1.02% 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 46,539, pits 5,781, largest drainage area 13.03 km2
filled.tif: flats 0, pits 195, largest drainage area 513.42 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_max_procs(1)                      # one thread: parallel runs give slightly different results
    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 83,652 (1.00%); median depth 1.273 m, deepest 20.47 m

Example 4 — build the DEM these examples use

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.merge import merge
from rasterio.transform import from_origin
from rasterio.warp import Resampling, reproject

basin = gpd.read_file("nldi_basin_01362500.geojson")            # USGS NLDI basin for the Coldbrook gauge
w, s, e, n = basin.total_bounds + np.array([-0.03, -0.03, 0.03, 0.03])
tiles = [rasterio.open(f"/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/{t}/USGS_13_{t}.tif")
         for t in ("n43w075", "n42w075")]
mosaic, mosaic_transform = merge(tiles, bounds=(w - 0.1, s - 0.1, e + 0.1, n + 0.1))

x0, y0, x1, y1 = basin.to_crs(32618).total_bounds + np.array([-1000, -1000, 1000, 1000])   # the basin plus 1 km
x0, y0, x1, y1 = np.floor(x0 / 10) * 10, np.floor(y0 / 10) * 10, np.ceil(x1 / 10) * 10, np.ceil(y1 / 10) * 10
width, height = int((x1 - x0) / 10), int((y1 - y0) / 10)
dem = np.full((height, width), -9999, "float32")
reproject(mosaic[0], dem, src_transform=mosaic_transform, src_crs=tiles[0].crs, dst_transform=from_origin(x0, y1, 10, 10),
          dst_crs="EPSG:32618", resampling=Resampling.bilinear, src_nodata=tiles[0].nodata, dst_nodata=-9999)
with rasterio.open("esopus_3dep13_utm_basin.tif", "w", driver="GTiff", height=height, width=width, count=1, dtype="float32",
                   crs="EPSG:32618", transform=from_origin(x0, y1, 10, 10), nodata=-9999, compress="deflate") as dst:
    dst.write(dem, 1)
print(f"{width} x {height} cells, {int((dem == -9999).sum())} without data")

3471 x 2408 cells, 0 without data

The 1 km margin matters more than it looks: see the note on clipping under Edge cases.

Why raw elevation cannot route water

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,781 such pits in the basin, the longest connected flow path drained 13.03 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 1.28 m 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.

Library versions, reservoirs, karst and DEM edges

• 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.
• Multithreaded WhiteboxTools runs differ: two runs of the same conditioning here differed in 10,471 cells, by up to 9.8 m. wbt.set_max_procs(1) makes them repeat.

Internal links

• Flow direction explained: D8, D-infinity and MFD compared — the routing step in detail
• How to condition a DEM for hydrology: breaching or filling — removing depressions
• How to calculate flow accumulation from a DEM in Python — the accumulation step
• Pour points explained: why watershed delineation needs a snapped outlet — choosing the outlet
• How to delineate a watershed from a pour point in Python — the full delineation
• How to extract a stream network from a DEM — thresholds and channels
• Digital elevation models explained: DEM, DSM and DTM — choosing the source
• How to fill sinks and derive flow direction from a DEM — a shorter introduction

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 13.03 km²; after conditioning about 1% of cells, the watershed at the gauge was 492.2 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.9% 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 78 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.