🔧 Fix Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 16 Notebook By · Updated

Fixing a Watershed That Comes Out as a Few Pixels

On this page

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 cell received flow from nothing but itself, and the catchment traced from the gauge's coordinates was 3 cells, 0.0003 km².
  • Snapping to cells draining more than 0.001 or 0.01 km² moved the outlet 30–37 m onto gullies, 0.01 km² at most; snapping to cells above 0.1 km² moved it 47.2 m onto Esopus Creek and gave 491.02 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 62.2 m and gave 492.21 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.

Causes of a watershed that comes out as a few pixels, and the fix for each. why the watershed is tinycausefixoutlet on the banksnap to a channelthreshold finds a gullyraise it: >= 0.1 km2snap distance too short100 m reached the creekdegrees on a UTM gridtransform the outletx and y swappedalways_xy=Trueforeign D8 codespass dirmapColdbrook gauge, conditioned 3DEP 10 m DEM spatialworkflow.io
Tiny watersheds almost always come from the outlet, not the DEM.

Read the accumulation at the outlet, then snap properly

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.05 km² — was 158 m away down the creek, and the nearest cell of Esopus Creek itself was 47.2 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 cells (0.001 km²) found a gully 30.2 m away draining a few cells, 100 cells (0.01 km²) another 37.2 m away with a 0.01 km² catchment; 1,000 cells (0.1 km²) or more found the creek, 47.2 m away, and 491.02 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 4.9 m to its cell centre; at 30 m and 50 m it found small flow lines 13.5 m and 33.0 m away. All four gave 0.00 km². At 100 m it reached the creek 62.2 m away: 492.21 km². At 300 m, 154.7 m away: 492.24 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 (535740.0, 4645810.0, 576450.0, 4675890.0).. 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.87 km², and with WhiteboxTools' codes, dirmap=(128, 1, 2, 4, 8, 16, 32, 64), 492.21 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.02 km² against 497.3 km², −1.3%, 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.

Catchment area at the Coldbrook gauge for the unsnapped outlet, pysheds snapping thresholds and WhiteboxTools snap distances. catchment area by snapping choiceunsnapped0.0003 km2pysheds > 0.01 km20.01 km2, gullypysheds > 0.1 km2491.02 km2whitebox 50 m0.00 km2whitebox 100 m492.21 km2whitebox 300 m492.24 km20500 km2Published drainage area 497.3 km2 spatialworkflow.io
Areas jump from almost nothing to the full catchment once snapping reaches the creek; there is no middle ground.

Checking outlets and snap distances

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; the same 3DEP data at 10 m over the basin plus 4 km (esopus_3dep13_utm.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", 4000, "esopus_3dep13_utm.tif")         # the same, basin + 4 km

for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm.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.tif: 4,071 x 3,008 cells of 10 m

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.05 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.0003 km2
snap to channels > 0.001 km2: moved  30.2 m, catchment    0.00 km2
snap to channels >  0.01 km2: moved  37.2 m, catchment    0.01 km2
snap to channels >   0.1 km2: moved  47.2 m, catchment  491.02 km2
snap to channels >    10 km2: moved  47.2 m, catchment  491.02 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_max_procs(1)                      # one thread: parallel runs give slightly different results
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   4.9 m, watershed    0.00 km2
snap_dist  30 m: moved  13.5 m, watershed    0.00 km2
snap_dist  50 m: moved  33.0 m, watershed    0.00 km2
snap_dist 100 m: moved  62.2 m, watershed  492.21 km2
snap_dist 300 m: moved 154.7 m, watershed  492.24 km2
pysheds default codes  492.87 km2
WhiteboxTools codes    492.21 km2

Why the gauge cell drains nothing

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

D8 direction codes in the ESRI and pysheds, WhiteboxTools and TauDEM conventions. one direction, three numbersdirectionESRI, pyshedsWhiteboxToolsTauDEMeast121south-east248south487south-west8166west16325north-west32644north641283north-east12812pysheds dirmap order is N, NE, E, SE, S, SW, W, NW spatialworkflow.io
Same eight directions, three numbering schemes; a pointer is meaningless without knowing which.

Confluences, dams and old datums

  • 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.
  • One conditioning pass may not be enough. pysheds' fill_pits, fill_depressions and resolve_flats, run once on the Esopus DEM, gave a snapped watershed of 459.0 km², 7.7% short and plausible-looking; running fill_depressions and resolve_flats a second time, as Example 1 does, gave 491.0 km².
  • 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 own cell drained only itself; snapped 47.2 m onto the creek, the catchment was 491.02 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.