🔧 Fix Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 18 Notebook By · Updated

Fixing Flow Accumulation That Is All Zeros, NoData or Stops at the Edge

On this page

Problem statement

Flow accumulation should light up a river network: small values on hillslopes, huge values along channels. Instead it looks black, or it has holes, or the main river fades out halfway down the valley, or the value at a gauge is a fraction of its drainage area. The accumulation algorithm is rarely at fault. The display, the DEM's NoData handling, the extent of the DEM or the coding of the flow directions usually is.

Reproduced on the 10 m 3DEP DEM around Esopus Creek, New York, with WhiteboxTools (fill with flat fixing, D8 pointer, accumulation in cells), where the correct area at the Coldbrook gauge is 492.24 km²:

  • The raster looked empty because 99.74% of cells were below 1% of the maximum: median 10 cells, maximum 5.6 million.
  • Removing the NoData tag from a reprojected DEM turned 147,932 cells of −999999 into a deep pit, and the gauge fell from 492.38 to 447.95 km².
  • A 300 m void across the main river, 5 km upstream, cut the gauge to 54.74 km², whether the void was tagged NoData or NaN.
  • Cropping the DEM 53.7 km² into the basin cost 88.9 km² at the gauge: valleys cut by the new edge drained out through it.
  • An ESRI-coded pointer read as WhiteboxTools codes gave 492.67 km² at the gauge — almost right, and therefore easy to miss.

Quick answer

symptom                                        cause                                        fix
raster looks black or all zeros                linear stretch of a skewed range              display log10(acc), or a stretch by percentile
holes or NoData streaks along rivers           voids in the DEM                              fill voids before conditioning
gauge area far too small, flow ends in a pit   NoData not tagged; fill values read as ground set the nodata tag, or mask the fill value
main river fades out towards the edge          DEM cropped through the catchment             extend the DEM beyond the divides
values plausible but slightly wrong            pointer read with another tool's codes        pass esri_pntr or dirmap to match the source

Check the accumulation at a point with a known drainage area before trusting the grid.

Drainage area at the gauge under seven DEM and processing faults. area at the gauge, km²what was doneareaerrorcorrect DEM492.24referenceNoData tag removed447.95−9.0%300 m void on the river54.74−88.9%masked to basin first490.48−0.4%cropped at 545,000 m403.33−18.1%cropped at 552,000 m293.81−40.3%ESRI codes as Whitebox492.67+0.1%WhiteboxTools fill, D8 pointer and accumulation on the 3DEP 10 m DEM; the NoData case against its own tagged run, 492.38 spatialworkflow.io
Only two of the seven are within a percent, and every one of them ran without an error.

Distribution, a known point, tagged NoData, filled voids

1. Look at the distribution, not the picture

Accumulation is extremely skewed. On the Esopus DEM the maximum was 5,624,605 cells, the median 10, the 90th percentile 67 and the 99th 2,855; 10.9% of cells had a value of 1, meaning nothing flowed into them. A linear colour ramp puts 99.74% of cells in its lowest 1%, so the raster looks black with a few faint lines. Display log10(acc) or stretch between percentiles (Example 1).

2. Check the value at a known point

The Coldbrook gauge has a published area of 497.3 km². The largest accumulation within 150 m of it was 4,922,396 cells, 492.24 km². A value orders of magnitude lower means the flow is going somewhere else; see fixing a watershed that comes out as a few pixels if only the snapping is wrong.

3. Make sure the DEM's NoData is tagged

3DEP tiles are on a geographic grid. Reprojected to UTM, the rectangle becomes a slightly rotated one, and the corners of the new grid that fall outside it are stored as −999999 with a NoData tag — 147,932 cells here. With the tag removed, those cells became real elevations 1,000 km below sea level. Filling does not remove a pit that deep; the edge wedges swallowed flow, and the gauge area fell from 492.38 km² to 447.95 km², 9.0% short. Check the tag, and look for common fill values — −9999, −32768, −3.4 × 10³⁸ — in the data (Example 2).

4. Fill voids before routing

A void inside the DEM is a hole in the surface. A 300 m square of NoData placed across Esopus Creek, where it already drained 373.5 km², stopped all of that flow: the gauge dropped to 54.74 km². Flow that reached the void went no further, exactly as if the grid ended there. Encoding the void as NaN instead of tagged NoData made no real difference: 54.71 km². Interpolate voids — rasterio.fill.fillnodata or GDAL's gdal_fillnodata — before conditioning.

5. Use a DEM that covers the whole catchment

Filling treats the grid edge as an outlet. Cropping the DEM at x = 545,000 m removed 53.7 km² of the basin, but the gauge lost 88.9 km², because valleys cut by the crop drained west out of the new edge — 27.66 km² of flow crossed it at a single cell. Cropping at 552,000 m removed 188.8 km² and cost 198.4 km². Clip to a buffered catchment, not a rectangle drawn around the river.

6. Do not mask the DEM to the catchment before routing

Masking to the NLDI basin polygon first gave 490.48 km² instead of 492.24: wherever the polygon's boundary was inside the DEM's true divide, the lost cells could no longer contribute. Route on the full DEM and clip the results afterwards.

7. Match the flow direction codes

ESRI numbers east as 1 and continues clockwise; WhiteboxTools numbers north-east as 1. An ESRI-coded pointer accumulated as WhiteboxTools codes gave 492.67 km² at the gauge; with esri_pntr=True, 492.24 km². The wrong codes rotate every direction by 45°, yet the gauge area moved by only 0.09% — so a mismatch can hide behind a plausible number. See flow direction explained.

8. Check the data type of the output

Accumulation in cells exceeds 65,535 on any sizeable basin. Written as or cast to uint16, the maximum wraps round, and as uint8 almost every channel becomes noise. Keep accumulation as float32 or int32, and use float64 when multiplying by cell area.

Flow accumulation percentiles on the Esopus DEM on a logarithmic scale. accumulation in cells, log scalemedian10 cells90th percentile67 cells99th percentile2,855 cellsmaximum5,624,605 cells07Bar length is log10 of the cell count; 99.74% of cells are below 1% of the maximum spatialworkflow.io
Five orders of magnitude between the median cell and the outlet: no linear colour ramp can show both.

Untagged NoData, voids and a cropped DEM

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); the 3DEP data on its native NAD83 latitude-longitude grid (esopus_3dep13_geo.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
if not os.path.exists("esopus_3dep13_geo.tif"):          # 3DEP on its native NAD83 grid, basin + 1 km
    srcs = [rasterio.open(p) for p in TILES["3dep13"][0]]
    bx = basin.to_crs(32618).total_bounds
    geo, geo_tf = merge(srcs, bounds=transform_bounds("EPSG:32618", srcs[0].crs, bx[0] - 1000, bx[1] - 1000, bx[2] + 1000, bx[3] + 1000))
    with rasterio.open("esopus_3dep13_geo.tif", "w", driver="GTiff", height=geo.shape[1], width=geo.shape[2], count=1,
                       dtype="float32", crs=srcs[0].crs, transform=geo_tf, nodata=-9999, compress="deflate") as dst:
        dst.write(np.where(geo[0] == srcs[0].nodata, -9999, geo[0]).astype("float32"), 1)

for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm.tif', 'esopus_3dep13_geo.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
esopus_3dep13_geo.tif: 4,557 x 2,367 cells of 9.25926e-05 degrees

Example 1 — route a DEM and describe the accumulation

import os

import numpy as np
import rasterio
import whitebox
from pyproj import Transformer

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("."))
gx, gy = Transformer.from_crs("EPSG:4269", "EPSG:32618", always_xy=True).transform(-74.2701944, 42.0144722)


def route(dem, tag):
    wbt.fill_depressions(dem, f"{tag}_filled.tif", fix_flats=True)
    wbt.d8_pointer(f"{tag}_filled.tif", f"{tag}_d8.tif")
    wbt.d8_flow_accumulation(f"{tag}_d8.tif", f"{tag}_acc.tif", out_type="cells", pntr=True)
    with rasterio.open(f"{tag}_acc.tif") as src:
        acc, transform = src.read(1, masked=True), src.transform
    col, row = ~transform * (gx, gy)
    r, c = int(row), int(col)
    return acc, transform, float(acc[r - 15:r + 16, c - 15:c + 16].max()) * abs(transform.a * transform.e) / 1e6


acc, transform, gauge_km2 = route("esopus_3dep13_utm.tif", "base")
cell_km2 = abs(transform.a * transform.e) / 1e6
values = acc.compressed()
print(f"gauge {gauge_km2:.2f} km2; max {values.max():,.0f} cells; median {np.median(values):.0f}, "
      f"90th {np.percentile(values, 90):.0f}, 99th {np.percentile(values, 99):.0f}; "
      f"{(values < 0.01 * values.max()).mean():.2%} below 1% of the maximum; {(values == 1).mean():.1%} equal to 1")
log_acc = np.ma.log10(acc)                                    # what to display instead
print(f"log10 range {log_acc.min():.1f} to {log_acc.max():.1f}")
gauge 492.24 km2; max 5,624,605 cells; median 10, 90th 67, 99th 2855; 99.74% below 1% of the maximum; 10.9% equal to 1
log10 range 0.0 to 6.8

Example 2 — NoData that is not tagged, and a void

from rasterio.warp import Resampling, calculate_default_transform, reproject

# 3DEP is delivered on a geographic grid; reprojected to UTM, the corners outside it become NoData
with rasterio.open("esopus_3dep13_geo.tif") as src:
    utm_transform, width, height = calculate_default_transform(src.crs, "EPSG:32618", src.width, src.height,
                                                               *src.bounds, resolution=10)
    utm_profile = dict(src.profile, crs="EPSG:32618", transform=utm_transform, width=width, height=height,
                       nodata=-999999.0)
    utm = np.full((height, width), -999999.0, "float32")
    reproject(rasterio.band(src, 1), utm, dst_transform=utm_transform, dst_crs="EPSG:32618",
              dst_nodata=-999999.0, resampling=Resampling.bilinear)
with rasterio.open("dem_utm.tif", "w", **utm_profile) as dst:
    dst.write(utm, 1)
print(f"nodata tag {utm_profile['nodata']}; cells equal to it {(utm == -999999).sum():,}; "
      f"cells below -1000 m {(utm < -1000).sum():,}")
tagged_km2 = route("dem_utm.tif", "tagged")[2]
with rasterio.open("dem_untagged.tif", "w", **dict(utm_profile, nodata=None)) as dst:
    dst.write(utm, 1)
print(f"with the tag: gauge {tagged_km2:.2f} km2; without it: {route('dem_untagged.tif', 'untagged')[2]:.2f} km2")

with rasterio.open("esopus_3dep13_utm.tif") as src:
    profile, dem = src.profile, src.read(1)
filled = acc.filled(0)
at_gauge = gauge_km2 / cell_km2
candidates = np.argwhere((filled > 0.6 * at_gauge) & (filled < 0.8 * at_gauge))    # the main river upstream
col, row = ~transform * (gx, gy)
hr, hc = candidates[np.hypot(candidates[:, 0] - row, candidates[:, 1] - col).argmin()]
void = dem.copy()
void[hr - 15:hr + 15, hc - 15:hc + 15] = profile["nodata"]
with rasterio.open("dem_void.tif", "w", **profile) as dst:
    dst.write(void, 1)
print(f"300 m void where the river drains {filled[hr, hc] * cell_km2:.1f} km2: "
      f"gauge {route('dem_void.tif', 'void')[2]:.2f} km2")
nodata tag -999999.0; cells equal to it 147,932; cells below -1000 m 147,932
with the tag: gauge 492.38 km2; without it: 447.95 km2
300 m void where the river drains 373.5 km2: gauge 54.74 km2

Example 3 — a DEM that stops inside the catchment

import geopandas as gpd
from rasterio.features import rasterize
from rasterio.windows import Window
from rasterio.windows import transform as window_transform

basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(profile["crs"])
inside = rasterize(basin.geometry, out_shape=dem.shape, transform=transform).astype(bool)
for cut_x in (545_000, 552_000):
    col0 = int((~transform * (cut_x, transform.f))[0])
    window = Window(col0, 0, dem.shape[1] - col0, dem.shape[0])
    with rasterio.open("dem_cropped.tif", "w", **dict(profile, width=int(window.width),
                                                     transform=window_transform(window, transform))) as dst:
        dst.write(dem[:, col0:], 1)
    cropped_acc, _, cropped_km2 = route("dem_cropped.tif", "cropped")
    print(f"cropped at x={cut_x:,}: removed {inside[:, :col0].sum() * cell_km2:.1f} km2 of the basin, "
          f"gauge {cropped_km2:.2f} km2 (lost {gauge_km2 - cropped_km2:.1f}); "
          f"largest flow out of the new west edge {cropped_acc[:, 0].max() * cell_km2:.2f} km2")
cropped at x=545,000: removed 53.7 km2 of the basin, gauge 403.33 km2 (lost 88.9); largest flow out of the new west edge 27.66 km2
cropped at x=552,000: removed 188.8 km2 of the basin, gauge 293.81 km2 (lost 198.4); largest flow out of the new west edge 48.14 km2

Why accumulation looks empty

Why accumulation looks empty

The number of cells draining to a point grows with distance downstream roughly like area, so a handful of channel cells hold values millions of times larger than the hillslope cells around them. Most viewers stretch colours linearly between minimum and maximum, which leaves everything but the largest river in the darkest colour. Nothing is wrong with the values.

Why untagged NoData is worse than a void

A tagged NoData cell is excluded from routing: flow reaching it leaves the grid, which is wrong but local. An untagged −999999 is treated as ground 1,000 km deep. Depression filling raises pits only to their spill point, and a pit that deep has its spill point at the grid edge, so large areas near the edge are routed into it.

Why a crop loses more than it removes

Depression filling and flat resolution assume the grid edge is an outlet. When a crop cuts across a valley, the valley's lower end is now the edge, and everything upstream of the cut drains out there — including land inside the retained part of the grid that used to flow east into the main river.

Why masking to the basin polygon costs area

A published basin polygon and a DEM's divides never coincide exactly. Cells between the two boundaries are NoData after masking; where the DEM's true divide lies outside the polygon, flow from those cells no longer reaches the channel, and where it lies inside, the masked edge becomes a new outlet.

Area removed by cropping the DEM against area lost at the gauge. removed by the crop against lost at the gaugeremoved by the croplost at the gaugecrop at 545,000 m53.7 km²✕88.9 km²crop at 552,000 m188.8 km²✕198.4 km²Largest flow out of the new west edge: 27.66 and 48.14 km² spatialworkflow.io
You lose more than you cut, because water that used to enter no longer can.

Pointer arguments, tiles and lakes as NoData

  • Passing a DEM where a pointer is expected is not always caught: WhiteboxTools, given the DEM with pntr=True, was still running after 10 minutes, when it was stopped.
  • Accumulation in pysheds includes the cell itself, so the minimum is 1; some tools start at 0.
  • Weighted accumulation (for example runoff) can legitimately contain zeros where weights are zero.
  • Tiled processing without overlap creates artificial edges on every tile; process whole catchments or use overlapping tiles.
  • Lakes as NoData in some DEMs split river networks; fill them with the water surface elevation.
  • Geographic DEMs give correct cell counts but wrong areas; see catchment areas on a latitude–longitude DEM.
  • Compressed output does not change values, but some viewers read the first overview level only.

FAQ

Why does my flow accumulation raster look all black?

The values span five orders of magnitude and a linear stretch hides all but the largest. On the Esopus DEM, 99.74% of cells were below 1% of the maximum; display log10 of the accumulation.

Why is flow accumulation too small at my outlet?

Flow is leaving the grid or disappearing before it gets there: untagged NoData, a void, a DEM cropped through the catchment, or an outlet off the channel. A void across the river cut the Esopus gauge from 492.24 to 54.74 km².

Why does flow accumulation stop at the edge of my DEM?

The edge is an outlet for depression filling. A DEM that ends inside the catchment loses everything upstream of the cut and more: a crop removing 53.7 km² lost 88.9 km² at the gauge.

Should I clip the DEM to the watershed before calculating flow accumulation?

No. Masking to the basin polygon first cost 2 km² at the Esopus gauge; route on a DEM that extends beyond the divides and clip the outputs afterwards.

How do I fix NoData holes in a DEM before flow routing?

Interpolate them with rasterio.fill.fillnodata or gdal_fillnodata, then condition the DEM. Leaving them as NoData makes each one an outlet.

Can wrong flow direction codes still give sensible accumulation?

Yes, which makes them dangerous. An ESRI-coded pointer read as WhiteboxTools codes gave 492.67 km² instead of 492.24 km² at the gauge; always state the coding when passing a pointer between tools.