🔧 Fix Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 18 of 20 Updated

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

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.45 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 the DEM turned 204,930 cells of −999999 into a deep pit, and the gauge fell to 448.42 km².
  • A 300 m void across the main river, 5 km upstream, cut the gauge to 54.74 km², whether the void was −999999 or NaN.
  • Cropping the DEM 53.6 km² into the basin cost 88.4 km² at the gauge: valleys cut by the new edge drained out through it.
  • An ESRI-coded pointer read as WhiteboxTools codes gave 492.80 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.

Bar chart of accumulated area at the Coldbrook gauge for the correct DEM, a DEM with its NoData tag removed, a void on the main river, a DEM masked to the basin polygon, two cropped DEMs and a mis-coded pointer.
One number — the area at a gauge — exposed every failure except the display problem.

Step-by-step solution

1. Look at the distribution, not the picture

Accumulation is extremely skewed. On the Esopus DEM the maximum was 5,642,219 cells, the median 10, the 90th percentile 67 and the 99th 2,836; 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,924,517 cells, 492.45 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

The 3DEP file stores missing cells as −999999 with a NoData tag. With the tag removed, those 204,930 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 to 448.42 km². 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.7 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 −999999 changed nothing. 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.6 km² of the basin, but the gauge lost 88.4 km², because valleys cut by the crop drained west out of the new edge — 27.67 km² of flow crossed it at a single cell. Cropping at 552,000 m removed 188.6 km² and cost 198.6 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.45 km² instead of 492.45: 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.80 km² at the gauge; with esri_pntr=True, 492.45 km². The wrong codes rotate every direction by 45°, yet the gauge area moved by only 0.07% — 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.

Bar chart of flow accumulation percentiles on the Esopus DEM — median, 90th, 99th percentile and maximum — on a logarithmic axis.
Five orders of magnitude between the median cell and the outlet: no linear colour ramp can show both.

Code examples

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_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.45 km2; max 5,642,219 cells; median 10, 90th 67, 99th 2836; 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 

with rasterio.open("esopus_3dep13_utm.tif") as src:
    profile, dem = src.profile, src.read(1)
print(f"nodata tag {profile['nodata']}; cells equal to it {(dem == profile['nodata']).sum():,}; "
      f"cells below -1000 m {(dem < -1000).sum():,}")

with rasterio.open("dem_untagged.tif", "w", **dict(profile, nodata=None)) as dst:
    dst.write(dem, 1)
print(f"without the nodata tag: gauge {route('dem_untagged.tif', 'untagged')[2]:.2f} km2 (correct {gauge_km2:.2f})")

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 204,930; cells below -1000 m 204,930
without the nodata tag: gauge 448.42 km2 (correct 492.45)
300 m void where the river drains 373.7 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.6 km2 of the basin, gauge 404.01 km2 (lost 88.4); largest flow out of the new west edge 27.67 km2
cropped at x=552,000: removed 188.6 km2 of the basin, gauge 293.90 km2 (lost 198.6); largest flow out of the new west edge 48.41 km2

Explanation

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.

Bar chart comparing basin area removed by two crops with the area lost at the gauge for each.
Crops lost more area at the gauge than they removed, because cut valleys drained out through the new edge.

Edge cases or notes

  • Passing a DEM where a pointer is expected is not always caught: WhiteboxTools, given the DEM with pntr=True, sat idle for 18 minutes until 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.45 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.6 km² lost 88.4 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.80 km² instead of 492.45 km² at the gauge; always state the coding when passing a pointer between tools.