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

Pour Points Explained: Why Watershed Delineation Needs a Snapped Outlet

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.40 km² rill; within 50 m, Esopus Creek at 492.0 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.

Bar chart of the Coldbrook watershed area when the outlet is snapped to the highest accumulation within radii from 0 to 5,000 metres.
Too small a radius finds the bank or a rill; too large a radius finds a different river.

Step-by-step solution

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 14 m away, 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          30 m      0.396 km2    -99.9%
50 m          50 m    492.023 km2     -1.1%
150 m        149 m    492.442 km2     -1.0%
1000 m      1000 m    494.070 km2     -0.6%
5000 m      2811 m    513.484 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.386 km² side channel 14 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%        22/87       64/87
max accumulation within 50 m                     1.3%        72/87       10/87
max accumulation within 150 m                    1.5%        61/87       15/87
max accumulation within 500 m                    8.6%        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.2%        85/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.45 km² at 150 m — the Esopus main stem next to its mouth — and Little Beaver Kill at Beechford (42.73 km²) became 491.99 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 26.5 m and found the 0.387 km² rill with snap_dist=50, and Esopus Creek 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 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.44 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.

Bar chart of the number of 87 gauges whose snapped watershed was within 5 percent of the published area, for no snapping, maximum accumulation within 50, 150 and 500 metres, and nearest channel within 150 metres.
The rule that looked safe at one big-river gauge was the worst across many.

Code examples

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.396 km2
max within   50 m:  492.023 km2
max within  150 m:  492.442 km2
max within 1000 m:  494.070 km2
nearest >= 0.1 km2:      0.386 km2
nearest >= 1.0 km2:    492.020 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.5%, within 5% 61, off by >50% 15
max within 500 m             87 gauges: median error   8.6%, 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.442 km2
nearest >= 1 km2, cell (2124, 2164), snap='corner': 492.020 km2
nearest >= 1 km2, cell (2124, 2164), snap='center': 492.020 km2

Explanation

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.2% median error — also shows how much of the remaining error comes from the DEM rather than the outlet.

Two panels showing a gauge on a tributary near its confluence: maximum-accumulation snapping jumps to the main river, nearest-channel snapping stays on the tributary.
Near a confluence the largest nearby accumulation belongs to the other river.

Edge cases or notes

  • 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.

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.