⚙️ How-To Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 6 of 20 Updated

How to Delineate a Watershed from a Pour Point in Python

Problem statement

You have a DEM and an outlet — a gauge, a bridge, a sampling point — and need the area that drains to it, as a polygon you can measure, map and overlay. The steps are standard: condition the DEM, route flow, snap the outlet to the channel, trace upstream, polygonise. What is not standard is how good the result is, and that depends mostly on the DEM and the outlet rather than on the code.

Measured for Esopus Creek above the USGS gauge at Coldbrook, New York — published drainage area 497.3 km² — against the USGS NLDI basin polygon:

  • With the 3DEP 10 m DEM, the watershed was 492.44 km², 0.97% under the published area, with an intersection-over-union of 0.985 against the NLDI polygon.
  • With the 30 m Copernicus GLO-30 surface model, it was 450.39 km², 9.4% short, IoU 0.897 and 51.3 km² of disagreement.
  • Conditioning took 18.7 s at 10 m and 0.22 s at 30 m; flow routing, snapping, tracing and polygonising together took 0.77 s and 0.17 s.
  • The 10 m polygon came out in 3 parts, and its geodesic area was 492.80 km² against 492.44 km² measured on the UTM grid.

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")                  # a conditioned, projected DEM
fdir = grid.flowdir(grid.read_raster("esopus_breached.tif"))
acc = grid.accumulation(fdir)
x, y = grid.snap_to_mask(acc > 10000, (560427.2, 4651640.7))   # outlet moved onto a channel > 1 km2
catchment = grid.catchment(x=x, y=y, fdir=fdir, xytype="coordinate")
grid.clip_to(catchment)
polygons = [shape(g) for g, v in grid.polygonize(grid.view(catchment, dtype=np.uint8)) if v == 1]
watershed = gpd.GeoDataFrame(geometry=polygons, crs=grid.crs.srs).dissolve()
print(f"{watershed.area.iloc[0] / 1e6:.2f} km2")
492.02 km2

Condition the DEM first, snap the outlet to a channel of the expected size, and compare the area with a published figure.

Flow for delineating a watershed: condition the DEM, route flow, snap the outlet, trace the catchment, polygonise and validate, with the time each step took on the 10 metre Esopus DEM.
Conditioning dominated the run time; the delineation itself took seconds.

Step-by-step solution

1. Choose and prepare the DEM

Use a bare-earth DEM in a projected CRS with metre cells, covering the whole catchment with margin. The 3DEP 1/3 arc-second DEM, reprojected to UTM zone 18N at 10 m, covered the basin in 2,406 × 3,470 cells; GLO-30 at 30 m in 802 × 1,157. Check that nodata is tagged: untagged nodata borders become cliffs that attract flow.

2. Condition it

Least-cost breaching with filling as a fallback — WhiteboxTools' breach_depressions_least_cost with a search of 1 km (100 cells at 10 m, 33 at 30 m) — took 18.68 s and 0.22 s. See breaching or filling a DEM for the options and their side effects.

3. Route flow

pysheds' flowdir and accumulation on the conditioned 10 m DEM took 0.42 s together. Accumulation is what the outlet snaps to and what the watershed's size can be read from before tracing it.

4. Snap the outlet

The gauge's coordinates lay beside the derived channel. The highest accumulation within 150 m was 160 m away diagonally on the 10 m grid and 163 m on the 30 m grid, on cells draining 492.44 km² and 450.39 km². For outlets near confluences, snap to the nearest channel of the expected size instead — see pour points explained.

5. Trace the catchment

grid.catchment walks flow directions upstream from the outlet and returns a boolean grid: 4,924,422 cells at 10 m, matching the accumulation at the outlet exactly. Snapping and tracing took 0.18 s.

6. Polygonise and dissolve

grid.polygonize on the clipped catchment returned the raster's outline as polygons; dissolved, the 10 m watershed had 3 parts — a main polygon and cells joined only at corners — with a total area of 492.44 km² on the UTM grid. Polygonising took 0.17 s. Remove slivers smaller than a few cells if you need a single polygon, and report that you did.

7. Measure area correctly

Area measured in UTM metres was 492.44 km²; the geodesic area of the same polygon on the WGS84 ellipsoid was 492.80 km². The difference, 0.07%, is the UTM scale factor at this location. For published figures use geodesic or equal-area measurements; see catchment areas on geographic DEMs.

8. Validate against a reference

The 10 m watershed was 0.97% under the published 497.28 km² and overlapped the NLDI polygon with an IoU of 0.9853, differing by 7.31 km². The 30 m watershed was 9.43% under, IoU 0.8972, differing by 51.26 km². The two DEM watersheds disagreed with each other by 47.35 km² (IoU 0.9044): the source, not the software, made the difference.

Bar chart of watershed area and error for the 3DEP 10 metre DEM and the GLO-30 30 metre surface model against the published drainage area.
The same code on two DEMs: under 1% error at 10 m bare earth, over 9% at 30 m surface.

Code examples

Example 1 — delineate and polygonise, with timings 

import os
import time

import rasterio
import whitebox
from pyproj import Geod, Transformer


def delineate(raw_dem, lon, lat, radius_m=150, workdir="work"):
    os.makedirs(workdir, exist_ok=True)
    timings = {}
    start = time.perf_counter()
    wbt = whitebox.WhiteboxTools()
    wbt.set_verbose_mode(False)
    with rasterio.open(raw_dem) as src:
        cell = abs(src.res[0])
    conditioned = os.path.abspath(os.path.join(workdir, os.path.basename(raw_dem).replace(".tif", "_cond.tif")))
    wbt.breach_depressions_least_cost(os.path.abspath(raw_dem), conditioned, dist=int(1000 / cell), fill=True)
    timings["condition"] = time.perf_counter() - start

    start = time.perf_counter()
    grid = Grid.from_raster(conditioned)
    fdir = grid.flowdir(grid.read_raster(conditioned))
    acc = np.asarray(grid.accumulation(fdir))
    timings["route"] = time.perf_counter() - start

    start = time.perf_counter()
    x, y = Transformer.from_crs("EPSG:4269", grid.crs.srs, always_xy=True).transform(lon, lat)
    col, row = ~grid.affine * (x, y)
    k = int(np.ceil(radius_m / cell))
    window = acc[int(row) - k:int(row) + k + 1, int(col) - k:int(col) + k + 1]
    dr, dc = np.unravel_index(window.argmax(), window.shape)
    ox, oy = grid.affine * (int(col) - k + dc + 0.5, int(row) - k + dr + 0.5)
    catchment = grid.catchment(x=ox, y=oy, fdir=fdir, xytype="coordinate", snap="center")
    timings["snap + trace"] = time.perf_counter() - start

    start = time.perf_counter()
    grid.clip_to(catchment)
    parts = [shape(g) for g, v in grid.polygonize(grid.view(catchment, dtype=np.uint8)) if v == 1]
    watershed = gpd.GeoDataFrame(geometry=parts, crs=grid.crs.srs).dissolve()
    timings["polygonise"] = time.perf_counter() - start

    geodesic = abs(Geod(ellps="WGS84").geometry_area_perimeter(watershed.to_crs(4326).geometry.iloc[0])[0]) / 1e6
    print(f"{os.path.basename(raw_dem)}: {len(parts)} part(s), {watershed.area.iloc[0] / 1e6:.2f} km2 projected, "
          f"{geodesic:.2f} km2 geodesic; outlet moved {np.hypot(ox - x, oy - y):.0f} m")
    print("  " + ", ".join(f"{step} {seconds:.2f} s" for step, seconds in timings.items()))
    return watershed


w10 = delineate("esopus_3dep13_utm_basin.tif", -74.2701944, 42.0144722)
w30 = delineate("esopus_glo30_utm_basin.tif", -74.2701944, 42.0144722)
esopus_3dep13_utm_basin.tif: 3 part(s), 492.44 km2 projected, 492.80 km2 geodesic; outlet moved 160 m
  condition 18.06 s, route 0.44 s, snap + trace 0.18 s, polygonise 0.18 s
esopus_glo30_utm_basin.tif: 1 part(s), 450.39 km2 projected, 450.72 km2 geodesic; outlet moved 163 m
  condition 0.23 s, route 0.05 s, snap + trace 0.04 s, polygonise 0.06 s

Example 2 — compare with a reference polygon and a published area 

def compare(watershed, reference, published_km2, label):
    a, b = watershed.geometry.iloc[0], reference.to_crs(watershed.crs).geometry.iloc[0]
    iou = a.intersection(b).area / a.union(b).area
    area = watershed.area.iloc[0] / 1e6
    print(f"{label}: {area:.2f} km2 ({100 * (area / published_km2 - 1):+.2f}% vs published), "
          f"IoU {iou:.4f}, symmetric difference {a.symmetric_difference(b).area / 1e6:.2f} km2")


nldi = gpd.read_file("nldi_basin_01362500.geojson")
compare(w10, nldi, 497.28, "3DEP 10 m")
compare(w30, nldi, 497.28, "GLO-30")
compare(w30, w10, w10.area.iloc[0] / 1e6, "GLO-30 against 3DEP")
3DEP 10 m: 492.44 km2 (-0.97% vs published), IoU 0.9853, symmetric difference 7.31 km2
GLO-30: 450.39 km2 (-9.43% vs published), IoU 0.8972, symmetric difference 51.26 km2
GLO-30 against 3DEP: 450.39 km2 (-8.54% vs published), IoU 0.9044, symmetric difference 47.35 km2

In the last line the "published" figure is the 3DEP watershed itself.

Example 3 — save the watershed with its provenance 

w10.assign(source="USGS 3DEP 1/3 arc-second, UTM 18N 10 m", conditioning="least-cost breach 1 km + fill",
           outlet="USGS 01362500, max accumulation within 150 m", area_km2=w10.area / 1e6).to_file("esopus_watershed.gpkg")
print(gpd.read_file("esopus_watershed.gpkg").iloc[0][["source", "conditioning", "outlet", "area_km2"]])
source                USGS 3DEP 1/3 arc-second, UTM 18N 10 m
conditioning                   least-cost breach 1 km + fill
outlet          USGS 01362500, max accumulation within 150 m
area_km2                                            492.4422
Name: 0, dtype: object

Keeping the DEM source, conditioning and snapping rule with the polygon makes the area reproducible and explains why it differs from someone else's.

Explanation

Why the DEM matters more than the tool

Every step after conditioning is deterministic: the same flow directions give the same catchment in any software. What changes the answer is where the DEM puts the drainage divides and channels. A 30 m surface model smooths narrow ridges, includes tree canopy and blurs valley floors; along a divide 100 km long, a shift of one or two cells moves square kilometres of area from one side to the other.

Why the polygon has several parts

A D8 catchment is a set of cells connected by flow, not by shared edges. Two cells can belong to the catchment while touching only at a corner, and polygonising treats corner-touching cells as separate polygons. Small extra parts are an artefact of the raster representation, not separate watersheds.

Why projected and geodesic areas differ

UTM preserves shapes locally but scales distances by up to 0.1% within a zone. Area measured on the grid inherits that scale squared. Geodesic area on the ellipsoid does not, which is why it was 0.36 km² larger here. Both are far more accurate than the DEM itself; the point is to state which one you report.

Why a good watershed still misses the published area

The published drainage area is itself an estimate from maps and older DEMs, and the NLDI polygon from NHDPlus data. The 10 m watershed matched both to about 1%, which is the level at which DEM resolution, conditioning and reference data all contribute. Differences above a few per cent point at a specific problem worth finding.

Bar chart of the time spent conditioning, routing, snapping and tracing, and polygonising for the 10 metre and 30 metre DEMs.
At 10 m the whole delineation after conditioning took under a second.

Edge cases or notes

  • Catchments larger than the DEM are truncated; check that no watershed cell touches the edge.
  • Reservoirs and lakes make flow paths across them arbitrary; mask or burn if the outlet is below one.
  • Transboundary basins may need DEMs from several sources; mosaic and reproject them to one grid first.
  • Very large basins at 10 m can exhaust memory in pysheds; process at coarser resolution or in WhiteboxTools, which works from files.
  • Several outlets are better handled together; see splitting a catchment into sub-basins.
  • Geographic DEMs need per-row cell areas; see catchment areas on a latitude–longitude DEM.
  • Simplifying the polygon for display changes its area; measure before simplifying.

FAQ

How do I delineate a watershed in Python?

Condition the DEM, compute flow direction and accumulation, snap the outlet to the channel, call pysheds' catchment or WhiteboxTools' watershed, then polygonise the result. For Esopus Creek at 10 m it matched the published area within 1%.

Which DEM should I use for watershed delineation?

A bare-earth DEM at the finest practical resolution. The 3DEP 10 m DEM was 0.97% off the published area; the 30 m GLO-30 surface model 9.43%.

How long does watershed delineation take?

For a 500 km² basin at 10 m, about 19 s to condition the DEM and under 1 s for routing, snapping, tracing and polygonising.

Why is my watershed polygon made of several parts?

Cells joined only at corners polygonise as separate parts. The Esopus watershed had 3 parts; dissolve them, and drop tiny slivers if a single polygon is needed.

How should I measure the watershed area?

Geodesically or in an equal-area projection for reporting. The UTM area was 492.44 km² and the geodesic area 492.80 km².

How do I check a delineated watershed?

Compare its area with a published drainage area and its boundary with a reference polygon. IoU against the USGS NLDI basin was 0.985 at 10 m.