Pour Points Explained: Why Watershed Delineation Needs a Snapped Outlet
On this page
Problem statement
A watershed is defined by its outlet: every cell that drains to one chosen cell. In practice the outlet comes from somewhere else — a gauge's coordinates, a bridge, a sampling site — and a DEM-derived channel is a line one cell wide that rarely passes exactly through it. Delineate from the raw coordinate and the "watershed" is a scrap of riverbank. Move the point too far and it lands on a bigger river, and the watershed is someone else's.
Measured on the conditioned 10 m DEM of Esopus Creek, New York, at the USGS gauge at Coldbrook (published 497.3 km²) and at 87 other gauges inside the basin with published drainage areas:
- The Coldbrook gauge's own cell drained 1 cell. The nearest cell draining more than 1 km² was 40 m away.
- Snapping to the highest accumulation within 30 m found a 0.39 km² rill; within 50 m, Esopus Creek at 491.9 km² — 1.1% under the published area.
- Across the 87 gauges, the highest accumulation within 500 m was more than 50% wrong at 28 of them, because the search jumped into larger rivers.
- Snapping to the nearest cell draining at least 0.5 km², within 150 m, got 79 of 87 within 5%, with a median error of 1.1%.
Quick answer
Snap to the nearest channel cell, not to the biggest one nearby, and check the result against a known area:
import numpy as np
def snap_to_channel(acc, transform, x, y, min_cells, radius_m):
"""Nearest cell with accumulation >= min_cells within radius_m of (x, y); returns row, col."""
col, row = ~transform * (x, y)
row, col = int(row), int(col)
k = int(np.ceil(radius_m / abs(transform.a)))
r0, c0 = max(row - k, 0), max(col - k, 0)
window = acc[r0:row + k + 1, c0:col + k + 1]
rr, cc = np.mgrid[r0:r0 + window.shape[0], c0:c0 + window.shape[1]]
distance = np.hypot(rr - row, cc - col) * abs(transform.a)
distance = np.where((window >= min_cells) & (distance <= radius_m), distance, np.inf)
i = np.unravel_index(distance.argmin(), distance.shape)
return (rr[i], cc[i]) if np.isfinite(distance[i]) else (row, col)Choose min_cells from the smallest stream the outlet could be on, and radius_m from the positional accuracy of the coordinates.
What the outlet is, and two ways to snap it
1. Understand what the outlet is
Watershed tools trace flow directions upstream from a cell. The pour point is that cell. Its accumulation is the watershed's size in cells, so the accumulation grid itself is the fastest check: a pour point on the right channel has an accumulation close to the expected drainage area.
2. Expect coordinates to miss the channel
Gauges are surveyed at the bank or the bridge; DEM channels follow the lowest cells, which may be one or several cells away, and may even be on the wrong side of a mid-channel bar. At Coldbrook, the gauge's cell had an accumulation of 1. Cells draining more than 0.01 km² started 10 m away — the neighbouring cell — and cells draining more than 1 km² 40 m away.
3. Snapping to maximum accumulation: fast and risky
Moving the point to the highest accumulation within a radius is the classic approach, and the default of WhiteboxTools' snap_pour_points. At Coldbrook:
radius moved watershed error
none 0 m 0.000 km2 -100.0%
30 m 22 m 0.390 km2 -99.9%
50 m 50 m 491.934 km2 -1.1%
150 m 149 m 492.354 km2 -1.0%
1000 m 998 m 494.057 km2 -0.6%
5000 m 2830 m 513.540 km2 +3.3%
From 50 m to 1,000 m the answer barely changed, because Esopus Creek was the largest channel in every search window. At 5 km the search reached the lowest cell of the DEM.
4. Snapping to the nearest channel: safer where rivers meet
Searching for the nearest cell whose accumulation exceeds a threshold keeps the outlet on the closest stream of at least that size. With thresholds of 0.01 and 0.1 km², the nearest qualifying cell at Coldbrook was a 0.389 km² side channel 10 m away; with 1 km² and above, Esopus Creek 40 m away. The threshold encodes what counts as "the river".
5. Test the rule on many outlets, not one
Coldbrook is on a large river with no bigger one nearby, so almost any rule works there. Across 87 gauges with published areas inside the basin:
rule median error within 5% off by >50%
no snapping 100.0% 23/87 64/87
max accumulation within 50 m 1.1% 72/87 10/87
max accumulation within 150 m 1.4% 63/87 15/87
max accumulation within 500 m 8.7% 36/87 28/87
nearest cell >= 0.5 km2 within 150 m 1.1% 79/87 4/87
nearest cell >= 0.5 km2 within 500 m 1.1% 79/87 3/87
best cell within 500 m (knowing the answer) 0.1% 84/87 0/87
The last line is an upper bound: the cell whose area best matched the published figure. Nearest-channel snapping came close to it.
6. Watch for jumps into a larger river
With maximum-accumulation snapping, Fox Hollow tributary at Allaben (10.31 km²) became 164.11 km² at 150 m — the Esopus main stem next to its mouth — and Little Beaver Kill at Beechford (42.73 km²) became 491.90 km² at 500 m. Gauges near confluences are exactly where outlets are most often placed and where this rule fails.
7. Know what your tool's snapping does
WhiteboxTools' snap_pour_points (maximum accumulation) moved Coldbrook 22.9 m and found the 0.39 km² rill with snap_dist=50, and Esopus Creek, 78.1 m away, with 150. jenson_snap_pour_points, which snaps to the nearest cell of a stream raster, did not move at all with 50 m — no stream cell was within reach — and found Esopus Creek, 42.5 m away, at 150 m. pysheds' snap_to_mask(acc > 100) found the rill; acc > 10000 the creek.
8. Check the snapped point's cell convention
pysheds' catchment() takes coordinates and by default snaps them to a cell corner. Given the exact centre of the cell chosen by maximum-accumulation snapping, snap="corner" resolved to the diagonal neighbour and returned a catchment under 0.001 km²; snap="center" returned the full 492.35 km². At the cell chosen by nearest-channel snapping both modes happened to agree, which is what makes the trap hard to spot. Coordinates computed as cell centres need snap="center" — see fixing a watershed that comes out as a few pixels.
Scoring snapping rules against published areas
The data used on this page
The setup builds every input this page reads, with the same code as DEM hydrology explained: 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); the 10 m DEM conditioned by least-cost breaching in WhiteboxTools (esopus_breached.tif).
import json
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
if not os.path.exists("esopus_breached.tif"): # least-cost breaching, 20 cells, one thread so it repeats
import whitebox
wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_max_procs(1)
wbt.breach_depressions_least_cost(os.path.abspath("esopus_3dep13_utm_basin.tif"), os.path.abspath("esopus_breached.tif"),
dist=20, fill=True)
if not os.path.exists("gauges_bbox.json"): # USGS stream gauges in the basin box, as the API returns them
Path("gauges_bbox.json").write_text(get(f"https://api.waterdata.usgs.gov/ogcapi/v0/collections/monitoring-locations/items?f=json"
f"&bbox={W},{S},{E},{N}&site_type_code=ST&limit=500").text)
for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm_basin.tif', 'esopus_breached.tif', 'gauges_bbox.json']:
if f.endswith(".json"):
print(f"{f}: {len(json.load(open(f))['features']):,} features")
elif 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
esopus_breached.tif: 3,471 x 2,408 cells of 10 m
gauges_bbox.json: 500 features
Example 1 — the two snapping rules at one gauge
import numpy as np
if not hasattr(np, "in1d"):
np.in1d = lambda a, b, **kw: np.isin(np.ravel(a), b, **kw)
from pyproj import Transformer
from pysheds.grid import Grid
grid = Grid.from_raster("esopus_breached.tif")
fdir = grid.flowdir(grid.read_raster("esopus_breached.tif"))
acc = np.asarray(grid.accumulation(fdir))
T, cell_km2 = grid.affine, abs(grid.affine.a * grid.affine.e) / 1e6
def snap_to_max(acc, transform, x, y, radius_m):
col, row = ~transform * (x, y)
row, col = int(row), int(col)
k = int(np.ceil(radius_m / abs(transform.a)))
r0, c0 = max(row - k, 0), max(col - k, 0)
window = acc[r0:row + k + 1, c0:col + k + 1]
rr, cc = np.mgrid[r0:r0 + window.shape[0], c0:c0 + window.shape[1]]
inside = np.hypot(rr - row, cc - col) * abs(transform.a) <= radius_m
i = np.unravel_index(np.where(inside, window, -1).argmax(), window.shape)
return rr[i], cc[i]
x, y = Transformer.from_crs("EPSG:4269", grid.crs.srs, always_xy=True).transform(-74.2701944, 42.0144722)
for radius in (30, 50, 150, 1000):
r, c = snap_to_max(acc, T, x, y, radius)
print(f"max within {radius:>4} m: {acc[r, c] * cell_km2:8.3f} km2")
for km2 in (0.1, 1.0):
r, c = snap_to_channel(acc, T, x, y, km2 / cell_km2, 150)
print(f"nearest >= {km2} km2: {acc[r, c] * cell_km2:8.3f} km2")max within 30 m: 0.390 km2
max within 50 m: 491.935 km2
max within 150 m: 492.354 km2
max within 1000 m: 494.057 km2
nearest >= 0.1 km2: 0.389 km2
nearest >= 1.0 km2: 491.931 km2
snap_to_channel is the function from the quick answer.
Example 2 — score both rules against published drainage areas
import json
import geopandas as gpd
from shapely.geometry import Point
MI2 = 2.589988110336
features = json.load(open("gauges_bbox.json"))["features"]
gauges = gpd.GeoDataFrame(
[{"name": f["properties"]["monitoring_location_name"], "km2": f["properties"]["drainage_area"] * MI2,
"geometry": Point(f["geometry"]["coordinates"])}
for f in features if (f["properties"].get("drainage_area") or 0) >= 1.0],
crs=4269).to_crs(grid.crs.srs)
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(grid.crs.srs).geometry.iloc[0].buffer(200)
gauges = gauges[gauges.within(basin)].reset_index(drop=True)
rules = {
"max within 150 m": lambda gx, gy: snap_to_max(acc, T, gx, gy, 150),
"max within 500 m": lambda gx, gy: snap_to_max(acc, T, gx, gy, 500),
"nearest >= 0.5 km2, 150 m": lambda gx, gy: snap_to_channel(acc, T, gx, gy, 0.5 / cell_km2, 150),
}
for name, rule in rules.items():
errors = np.array([abs(acc[rule(p.x, p.y)] * cell_km2 / km2 - 1) for p, km2 in zip(gauges.geometry, gauges.km2)])
print(f"{name:28} {len(errors)} gauges: median error {np.median(errors):6.1%}, within 5% {int((errors <= 0.05).sum())}, "
f"off by >50% {int((errors > 0.5).sum())}")max within 150 m 87 gauges: median error 1.4%, within 5% 63, off by >50% 15
max within 500 m 87 gauges: median error 8.7%, within 5% 36, off by >50% 28
nearest >= 0.5 km2, 150 m 87 gauges: median error 1.1%, within 5% 79, off by >50% 4
Example 3 — the cell-centre trap in pysheds
for label, (r, c) in [("max within 150 m", snap_to_max(acc, T, x, y, 150)),
("nearest >= 1 km2", snap_to_channel(acc, T, x, y, 1.0 / cell_km2, 150))]:
cx, cy = T * (c + 0.5, r + 0.5) # the centre of the snapped cell
for mode in ("corner", "center"):
catchment = np.asarray(grid.catchment(x=cx, y=cy, fdir=fdir, xytype="coordinate", snap=mode))
print(f"{label}, cell ({r}, {c}), snap={mode!r}: {catchment.astype(bool).sum() * cell_km2:.3f} km2")max within 150 m, cell (2138, 2163), snap='corner': 0.000 km2
max within 150 m, cell (2138, 2163), snap='center': 492.354 km2
nearest >= 1 km2, cell (2124, 2164), snap='corner': 491.931 km2
nearest >= 1 km2, cell (2124, 2164), snap='center': 491.931 km2
Why one cell changes everything
Why a single cell makes such a difference
Accumulation jumps across the width of a valley: the channel cell carries the whole upstream network, its neighbours only their own few metres of slope. There is no gentle gradient to land on approximately. The outlet is either on the channel or not.
Why maximum accumulation fails near confluences
A search window near the mouth of a tributary contains cells of both the tributary and the main river. The main river's accumulation is always larger, so the maximum rule chooses it whenever it is in range. The larger the radius, the more tributary gauges it captures: 15 of 87 were off by more than half at 150 m, 28 at 500 m.
Why nearest-channel snapping is more robust
It moves the point the shortest distance that puts it on a stream of meaningful size. Unless the tributary is smaller than the threshold, or the main river is nearer than the tributary's channel, it stays on the right stream. That made it better than every maximum-accumulation radius on the 87 gauges and much less sensitive to the radius.
Why validation data matter
A snapped outlet that finds a plausible river produces a plausible-looking watershed. Only a known drainage area, a reference boundary or a visual check against the mapped river network tells you whether it is the intended one. With published areas, the gap to the best achievable cell — 0.1% median error — also shows how much of the remaining error comes from the DEM rather than the outlet.
Reservoirs, braided channels and swapped coordinates
- Gauges on reservoirs or lakes sit on flats where the derived channel is arbitrary.
- Braided channels can split flow between parallel DEM channels; the outlet captures only one.
- Coordinates in a different CRS than the DEM put the point far outside the grid; pysheds raised a "pour point out of bounds"
ValueError. - Swapped x and y produce the same error or, worse, a point inside the grid somewhere else.
- Coarser DEMs move channels further from true positions and need larger radii.
- Multiple outlets that snap to the same cell produce identical watersheds; check for duplicates.
- Stream burning reduces snapping distance by forcing channels onto mapped positions.
Internal links
- How to delineate a watershed from a pour point in Python — the delineation itself
- Fixing a watershed that comes out as a few pixels — outlet and cell-convention failures
- How to calculate flow accumulation from a DEM in Python — the grid outlets snap to
- How to split a catchment into sub-basins at many outlets — snapping many gauges at once
- How to extract a stream network from a DEM — the threshold behind nearest-channel snapping
- DEM hydrology explained: from elevation to where water goes — the whole pipeline
- How to burn a known river network into a DEM — moving channels onto mapped rivers
- Projected vs geographic CRS: what actually changes when you reproject — keeping outlets and DEM in one CRS
FAQ
What is a pour point?
The outlet cell of a watershed: the cell whose upstream area the watershed is. Delineation traces flow directions upstream from it.
Why do I need to snap a pour point?
Because the derived channel is one cell wide and the outlet's coordinates rarely fall on it. The Coldbrook gauge's own cell drained only itself; the channel was 40 m away.
What snap distance should I use?
Enough to reach the channel from your coordinates, and no more. Maximum-accumulation snapping at 500 m was more than 50% wrong at 28 of 87 gauges; nearest-channel snapping within 150 m at 4.
What is the difference between snap_pour_points and jenson_snap_pour_points?
WhiteboxTools' snap_pour_points moves to the highest accumulation within the distance; jenson_snap_pour_points moves to the nearest cell of a stream raster. The second avoids jumping into larger rivers.
How do I know the snapped outlet is right?
Compare the accumulation at the snapped cell with a known drainage area. Across 87 gauges, nearest-channel snapping put 79 within 5% of their published areas.
Why does pysheds return a tiny catchment at my snapped point?
The default snap="corner" can resolve cell-centre coordinates to a neighbouring cell. Pass snap="center" when your coordinates are cell centres.