🔧 Fix Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 16 of 20 Updated

Fixing a Watershed That Comes Out as a Few Pixels

Problem statement

You pass the coordinates of a river gauge to a watershed tool and get back a speck: one cell, a handful of cells, a polygon a few metres across. Or the call fails with an "out of bounds" error. The river is real, the DEM is fine, and the outlet is exactly where the gauge is. That is the problem — the gauge is on the riverbank as the DEM sees it, not in the modelled channel.

Reproduced at the Esopus Creek gauge at Coldbrook, New York (USGS 01362500, published area 497.3 km²), on a conditioned 10 m 3DEP DEM:

  • The gauge's own coordinates gave a catchment of 1 cell, 0.0001 km²: that cell received flow from nothing but itself.
  • Snapping to cells draining more than 0.001 or 0.01 km² moved the outlet 34.8 m onto a gully, still 0.01 km²; snapping to cells above 0.1 km² moved it 44.8 m onto Esopus Creek and gave 491.24 km².
  • WhiteboxTools' snap_pour_points with a snap distance of 0, 10, 30 or 50 m gave 0.00 km²; 100 m moved the outlet 65.6 m and gave 492.42 km².
  • Longitude and latitude passed to a UTM grid, or x and y swapped, raised ValueError: Pour point ... is out of bounds rather than a small watershed.

Quick answer 

x, y = grid.snap_to_mask(acc > 0.1 / cell_km2, (x, y))     # onto a channel draining at least 0.1 km2
catchment = grid.catchment(x=x, y=y, fdir=fdir, xytype="coordinate")

Check the accumulated area at the outlet before tracing: if it is a few cells, the outlet is not on the channel. Snap to a channel of a plausible size — ideally one close to the gauge's published drainage area — and make sure the coordinates are in the DEM's CRS and in x, y order.

Triage of tiny watershed results: an outlet on the bank, snapping to too small a channel, too short a snap distance, coordinates in the wrong CRS or order, and a flow direction raster with the wrong direction codes.
Tiny watersheds almost always come from the outlet, not the DEM.

Step-by-step solution

1. Look at the accumulation at the outlet cell

Before tracing anything, read the flow accumulation at the outlet. At the Coldbrook gauge's coordinates it was 1 cell. A river gauge on a 497 km² catchment should sit on a cell draining around 5 million 10 m cells. Any value orders of magnitude below the expected area means the outlet is off the channel.

2. Find the channel near the outlet

Within 150 m of the gauge, the largest accumulation — 491.27 km² — was 158 m away down the creek, and the nearest cell of Esopus Creek itself was 44.8 m away (Example 2). Being that far off is normal: gauge coordinates mark a gauge house or a bank, often recorded to a few seconds of arc, while a DEM-derived channel can sit a cell or two from the real one.

3. Snap to a channel of the right size

grid.snap_to_mask(acc > threshold, (x, y)) moves the outlet to the nearest cell above the threshold. The threshold decides which channel: 10 or 100 cells (0.001 or 0.01 km²) found a gully 34.8 m away with a 0.01 km² catchment; 1,000 cells (0.1 km²) or more found the creek, 44.8 m away, and 491.24 km². Choose a threshold that small side channels near the gauge cannot reach — a fraction of the expected area works well.

4. Give snapping tools enough distance

WhiteboxTools' snap_pour_points moves the outlet to the highest accumulation within snap_dist. At 0 and 10 m it moved the point 2.3 m to its cell centre; at 30 m and 50 m it found small flow lines 15.6 m and 26.5 m away. All four gave 0.00 km². At 100 m it reached the creek 65.6 m away: 492.42 km². At 300 m, 160.2 m away: 492.45 km². Too large a distance risks jumping to a bigger river nearby; see pour points explained.

5. Check the CRS and the axis order

Longitude and latitude passed to a grid in UTM metres fell far outside its bounds, and pysheds raised ValueError: Pour point (-74.2701944, 42.0144722) is out of bounds for dataset with bbox (536302.54, 4645423.40, 576242.54, 4676833.40). Northing and easting swapped did the same. Transform the gauge into the DEM's CRS with always_xy=True.

6. Check row and column against coordinates

Passing a column and row as if they were coordinates also raised the out-of-bounds error. With xytype="index" they were accepted — and gave 1 cell, because the unsnapped cell was still off the channel. Snapping has to happen whichever form the outlet takes.

7. Check the flow direction codes

Direction rasters from different tools use different codes. A WhiteboxTools pointer read by pysheds with its default ESRI codes gave 492.99 km², and with WhiteboxTools' codes, dirmap=(128, 1, 2, 4, 8, 16, 32, 64), 492.41 km². A mismatch does not always shrink the watershed, which makes it easy to miss; set dirmap explicitly whenever the pointer came from another tool.

8. Compare with the published area

After snapping, compare the catchment with the gauge's published drainage area: 491.24 km² against 497.3 km², −1.2%, is a match. A catchment within a factor of two but not closer points to conditioning or a different outlet channel; a factor of 100 or more points back to the outlet.

Bar chart of catchment area at the Coldbrook gauge for the unsnapped outlet, pysheds snapping thresholds and WhiteboxTools snap distances.
Areas jump from almost nothing to the full catchment once snapping reaches the creek; there is no middle ground.

Code examples

Example 1 — check an outlet before tracing 

import numpy as np

if not hasattr(np, "in1d"):                                     # pysheds 0.5 on NumPy 2.4+
    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_3dep13_utm.tif")
dem = grid.read_raster("esopus_3dep13_utm.tif")
surface = grid.resolve_flats(grid.fill_depressions(grid.fill_pits(dem)))
surface = grid.resolve_flats(grid.fill_depressions(surface))       # a second pass drains what the first left
fdir = grid.flowdir(surface)
acc = grid.accumulation(fdir)
cell_m = abs(grid.affine.a)
cell_km2 = cell_m * abs(grid.affine.e) / 1e6


def check_outlet(x, y, radius_m=150):
    xmin, ymin, xmax, ymax = grid.bbox
    if not (xmin <= x <= xmax and ymin <= y <= ymax):
        hints = []
        if abs(x) <= 180 and abs(y) <= 90:
            hints.append("looks like longitude/latitude")
        if xmin <= y <= xmax and ymin <= x <= ymax:
            hints.append("x and y look swapped")
        return f"({x:.1f}, {y:.1f}) is outside the grid: " + ("; ".join(hints) or "wrong CRS or wrong area")
    col, row = ~grid.affine * (x, y)
    r, c, k = int(row), int(col), int(radius_m / cell_m)
    window = np.asarray(acc)[r - k:r + k + 1, c - k:c + k + 1] * cell_km2
    i = np.unravel_index(window.argmax(), window.shape)
    away = np.hypot(i[0] - k, i[1] - k) * cell_m
    return (f"cell drains {acc[r, c] * cell_km2:.4f} km2; largest within {radius_m} m drains "
            f"{window.max():.2f} km2, {away:.0f} m away")


x, y = Transformer.from_crs("EPSG:4269", grid.crs.srs, always_xy=True).transform(-74.2701944, 42.0144722)
for label, point in (("gauge in UTM", (x, y)), ("gauge in degrees", (-74.2701944, 42.0144722)), ("x and y swapped", (y, x))):
    print(f"{label:17} {check_outlet(*point)}")
gauge in UTM      cell drains 0.0001 km2; largest within 150 m drains 491.27 km2, 158 m away
gauge in degrees  (-74.3, 42.0) is outside the grid: looks like longitude/latitude
x and y swapped   (4651640.7, 560427.2) is outside the grid: x and y look swapped

Example 2 — snap to channels of different sizes 

def catchment_km2(x, y):
    return np.asarray(grid.catchment(x=x, y=y, fdir=fdir, xytype="coordinate")).sum() * cell_km2


print(f"unsnapped: {catchment_km2(x, y):.4f} km2")
for threshold_km2 in (0.001, 0.01, 0.1, 10):
    sx, sy = grid.snap_to_mask(acc > threshold_km2 / cell_km2, (x, y))
    print(f"snap to channels > {threshold_km2:>5} km2: moved {np.hypot(sx - x, sy - y):5.1f} m, "
          f"catchment {catchment_km2(sx, sy):7.2f} km2")
unsnapped: 0.0001 km2
snap to channels > 0.001 km2: moved  34.8 m, catchment    0.01 km2
snap to channels >  0.01 km2: moved  34.8 m, catchment    0.01 km2
snap to channels >   0.1 km2: moved  44.8 m, catchment  491.24 km2
snap to channels >    10 km2: moved  44.8 m, catchment  491.24 km2

Example 3 — snap distances and direction codes in WhiteboxTools 

import os

import geopandas as gpd
import rasterio
import whitebox
from pysheds.sview import Raster
from shapely.geometry import Point

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("."))
wbt.fill_depressions("esopus_3dep13_utm.tif", "filled.tif", fix_flats=True)
wbt.d8_pointer("filled.tif", "d8.tif")
wbt.d8_flow_accumulation("d8.tif", "acc.tif", pntr=True)
gpd.GeoDataFrame(geometry=[Point(x, y)], crs=32618).to_file("gauge.shp")

for snap_dist in (0, 30, 50, 100, 300):
    wbt.snap_pour_points("gauge.shp", "acc.tif", "snapped.shp", snap_dist=snap_dist)
    wbt.watershed("d8.tif", "snapped.shp", "watershed.tif")
    with rasterio.open("watershed.tif") as src:
        ws = src.read(1)
        km2 = ((ws != src.nodata) & (ws > 0)).sum() * cell_km2
    moved = gpd.read_file("snapped.shp").geometry.iloc[0].distance(Point(x, y))
    print(f"snap_dist {snap_dist:>3} m: moved {moved:5.1f} m, watershed {km2:7.2f} km2")

with rasterio.open("d8.tif") as src:
    whitebox_codes = Raster(src.read(1).astype("int64"), fdir.viewfinder)
sx, sy = grid.snap_to_mask(acc > 10 / cell_km2, (x, y))
for label, dirmap in (("pysheds default codes", (64, 128, 1, 2, 4, 8, 16, 32)),
                      ("WhiteboxTools codes", (128, 1, 2, 4, 8, 16, 32, 64))):
    traced = grid.catchment(x=sx, y=sy, fdir=whitebox_codes, dirmap=dirmap, xytype="coordinate")
    print(f"{label:22} {np.asarray(traced).sum() * cell_km2:.2f} km2")
snap_dist   0 m: moved   2.3 m, watershed    0.00 km2
snap_dist  30 m: moved  15.6 m, watershed    0.00 km2
snap_dist  50 m: moved  26.5 m, watershed    0.00 km2
snap_dist 100 m: moved  65.6 m, watershed  492.42 km2
snap_dist 300 m: moved 160.2 m, watershed  492.45 km2
pysheds default codes  492.99 km2
WhiteboxTools codes    492.41 km2

Explanation

Why the gauge cell drains nothing

A D8 flow direction grid sends each cell to one neighbour, and flow accumulates only along the chains that result. Channels are one cell wide. A gauge a few tens of metres from the modelled channel lands on a hillside cell whose upstream area is itself and perhaps a few neighbours. Tracing upstream from it returns exactly that.

Why snapping thresholds produce jumps, not gradual changes

Accumulation along a hillside is tiny and along a river is huge, with little in between near a confluence. A threshold either reaches only side gullies or reaches the main channel; the catchment size jumps accordingly, from 0.01 km² to 491.24 km² between thresholds of 100 and 1,000 cells.

Why wrong coordinates raise errors in pysheds but not always elsewhere

pysheds checks that the pour point lies inside the grid's bounding box. Some tools instead clamp to the nearest edge cell or return an empty raster, which looks like a tiny watershed. A check of the outlet against the raster bounds, as in Example 1, makes all tools behave the same.

Why direction codes need care

ESRI, WhiteboxTools, TauDEM and GRASS number the eight directions differently. A pointer read with the wrong codes rotates every flow direction by 45° or more. At the Coldbrook gauge the rotated field still traced an area within 0.2% of the real one, so the size of the result is no proof that the codes were right; check the codes themselves.

Table of D8 direction codes for north, north-east, east and the other directions in ESRI and pysheds, WhiteboxTools, and TauDEM conventions.
Same eight directions, three numbering schemes; a pointer is meaningless without knowing which.

Edge cases or notes

  • Gauges beside confluences can snap to the wrong river with a large threshold; compare with the published drainage area.
  • Gauges upstream of a dam or lake may snap across the water surface to an unexpected channel.
  • A flat, unconditioned DEM can produce a small watershed even at the right cell; condition first.
  • Coordinates in NAD27 can be tens of metres off in the United States; transform from the stated datum.
  • Outlets on NoData cells trace nothing; check the DEM value at the outlet.
  • Very large thresholds can snap a tributary gauge onto the main stem beside it.
  • Geographic DEMs need snapping distances in degrees, or reprojection first; see catchment areas on a latitude–longitude DEM.

FAQ

Why is my watershed only one pixel?

The outlet is not on the modelled channel, so nothing flows into it. The Coldbrook gauge's coordinates drained 1 cell; snapped 44.8 m onto the creek, the catchment was 491.24 km².

How far should I snap a pour point?

Far enough to reach the channel, not so far that it reaches a different river. At Coldbrook, 50 m was not enough and 100 m was; checking the result against the published area is more reliable than any fixed distance.

What threshold should I use with snap_to_mask?

One that small side channels near the gauge do not exceed. Thresholds of 0.001 and 0.01 km² snapped to a gully; 0.1 km² or more found Esopus Creek.

Why does pysheds say the pour point is out of bounds?

The coordinates are not in the grid's CRS or are in the wrong order. Longitude and latitude on a UTM grid, and swapped x and y, both raised that error.

Can wrong flow direction codes make a watershed small?

Yes, although not always. A WhiteboxTools pointer read with pysheds' default codes rotates every direction; always pass dirmap for a pointer from another tool.

How do I know the delineated watershed is right?

Compare its area with the published drainage area of the gauge. Within a few per cent is a match; orders of magnitude smaller means the outlet missed the channel.