⚙️ How-To Terrain & elevation Raster & rasterio Terrain & Elevation · step 11 of 12 Updated

How to Fill Sinks and Derive Flow Direction from a DEM

Problem statement

Every hydrological product — flow direction, accumulation, stream networks, watersheds — assumes water can always run downhill to the edge of the grid. The standard advice is to "fill the sinks" first.

Measure your DEM before following it, because the advice assumes a problem you may not have:

import numpy as np
from scipy.ndimage import label, minimum_filter

neighbour_min = minimum_filter(dem, size=3, mode="nearest")
interior = np.zeros(dem.shape, bool)
interior[1:-1, 1:-1] = True

true_pits = (dem < neighbour_min) & interior          # strictly lower than all 8
no_outlet = (dem <= neighbour_min) & interior         # no strictly-lower neighbour

print(f"true pits:            {true_pits.sum():,} in {label(true_pits)[1]:,} regions")
print(f"cells with no outlet: {no_outlet.sum():,} in {label(no_outlet)[1]:,} regions")
true pits:            0 in 0 regions
cells with no outlet: 1,393 in 130 regions

Zero true pits. Copernicus DEM is void-filled and conditioned, so no cell is strictly lower than all its neighbours. What it does have is 1,393 cells with nowhere to drain — cells on flats, where at least one neighbour is at exactly the same height.

That distinction matters, because filling fixes pits and makes flats worse.

Quick answer 

from skimage.morphology import reconstruction


def fill_sinks(dem):
    """Priority-flood: raise every depression to the level of its lowest outlet."""
    seed = np.full_like(dem, dem.max())
    seed[0, :], seed[-1, :] = dem[0, :], dem[-1, :]      # boundary stays put
    seed[:, 0], seed[:, -1] = dem[:, 0], dem[:, -1]
    return reconstruction(seed, dem, method="erosion")


filled = fill_sinks(dem)
raised = filled - dem
print(f"{(raised > 0).sum():,} cells raised, median {np.median(raised[raised > 0]):.2f} m, "
      f"max {raised.max():.2f} m")

after = ((filled <= minimum_filter(filled, size=3, mode="nearest")) & interior)
print(f"cells with no outlet after filling: {after.sum():,} in {label(after)[1]:,} regions")
2,057 cells raised, median 1.07 m, max 10.51 m
cells with no outlet after filling: 2,067 in 114 regions

The count went up, from 1,393 to 2,067. Filling did its job — the 130 depressions became 114 — and every depression it raised became a flat surface at its outlet level, where D8 has no steepest descent.

Stage What it fixes What it leaves
fill sinks cells strictly lower than all neighbours flats at every outlet level
resolve flats arbitrary direction on flat ground nothing
flow direction which neighbour each cell drains to —
accumulation how many cells drain through each —
A true pit strictly lower than all neighbours, and a flat where several cells share the same height, with filling converting the first into the second.
Two different problems. Filling converts the left one into the right one, and only the right one remains.

Step-by-step solution

1. Distinguish pits from flats before doing anything 

def drainage_report(surface, label_text):
    neighbour_min = minimum_filter(surface, size=3, mode="nearest")
    interior = np.zeros(surface.shape, bool)
    interior[1:-1, 1:-1] = True

    strict = (surface < neighbour_min) & interior
    stuck = (surface <= neighbour_min) & interior
    print(f"{label_text}")
    print(f"   true pits            {strict.sum():>6,} in {label(strict)[1]:>5,} regions")
    print(f"   no downhill neighbour {stuck.sum():>6,} in {label(stuck)[1]:>5,} regions")


drainage_report(dem, "RAW DEM")
RAW DEM
   true pits                 0 in     0 regions
   no downhill neighbour  1,393 in   130 regions

Two very different diagnoses:

  • Many true pits — a noisy DEM. Filling is exactly right and will remove almost all of them.
  • Zero true pits, many flats — a conditioned product like Copernicus DEM. Filling changes little and flat-resolution is your real problem.

Running the strict test alone would have reported "no pits, nothing to do", and then produced a fragmented stream network anyway.

2. Fill with priority-flood 

seed = np.full_like(dem, dem.max())
seed[0, :], seed[-1, :] = dem[0, :], dem[-1, :]
seed[:, 0], seed[:, -1] = dem[:, 0], dem[:, -1]
filled = reconstruction(seed, dem, method="erosion")

Morphological reconstruction by erosion is priority-flood: it lowers a surface starting at the maximum everywhere except the boundary until it rests on the DEM. The result is the lowest surface above the DEM with no interior minima — each depression raised to exactly its lowest outlet, and no higher.

raised = filled - dem
print(f"{(raised > 0).sum():,} raised, median {np.median(raised[raised > 0]):.2f} m, "
      f"max {raised.max():.2f} m")
2,057 raised, median 1.07 m, max 10.51 m

A median of 1.07 m is consistent with noise in a DEM of about 4 m vertical accuracy. The 10.51 m maximum deserves a look at the hillshade — that is either a real depression or a void-fill artefact.

Filling is destructive. Keep the original and never use the filled surface for slope, contours or elevation reporting.

3. Expect the flats to grow 

drainage_report(filled, "FILLED")
FILLED
   true pits                 0 in     0 regions
   no downhill neighbour  2,067 in   114 regions

Regions down from 130 to 114 — depressions merged and resolved. Cells up from 1,393 to 2,067, because every filled depression is now a plateau at its outlet level.

This is not a failure. It is what filling does, and it moves the problem from "water stops here" to "water does not know which way to go here".

4. Compute D8 flow direction 

OFFSETS = [(-1, 0), (-1, 1), (0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1)]
CODES = [64, 128, 1, 2, 4, 8, 16, 32]      # ESRI convention, clockwise from north


def d8_direction(surface, cell_x, cell_y):
    best_drop = np.full(surface.shape, -np.inf)
    direction = np.zeros(surface.shape, dtype="uint8")

    for (dr, dc), code in zip(OFFSETS, CODES):
        shifted = np.roll(np.roll(surface, -dr, axis=0), -dc, axis=1)
        distance = np.hypot(dr * cell_y, dc * cell_x)       # diagonals are further
        drop = (surface - shifted) / distance
        better = drop > best_drop
        best_drop = np.where(better, drop, best_drop)
        direction = np.where(better, code, direction).astype("uint8")

    direction[best_drop <= 0] = 0                            # flats and edges
    direction[0, :] = direction[-1, :] = 0                   # np.roll wraps — mask it
    direction[:, 0] = direction[:, -1] = 0
    return direction, best_drop

Dividing by the distance is not optional. Comparing raw height drops instead of gradients makes diagonal neighbours win whenever they are within √2 of the orthogonal drop, and the result is a distinctive herringbone flow pattern across every hillside.

5. Accumulate from the top down 

def d8_accumulation(surface, direction):
    rows, cols = surface.shape
    order = np.argsort(surface.ravel())[::-1]        # highest cells first
    accumulation = np.ones(surface.size)
    lookup = dict(zip(CODES, OFFSETS))
    flat_dir = direction.ravel()

    for idx in order:
        code = flat_dir[idx]
        if code == 0:
            continue
        dr, dc = lookup[code]
        r, c = divmod(int(idx), cols)
        nr, nc = r + dr, c + dc
        if 0 <= nr < rows and 0 <= nc < cols:
            accumulation[nr * cols + nc] += accumulation[idx]

    return accumulation.reshape(surface.shape)

Processing highest-first guarantees every upstream contribution has arrived before a cell is passed on. It works only because D8 gives each cell exactly one downstream neighbour, so the flow graph is a forest.

The eight D8 direction codes arranged around a centre cell, with diagonal distances scaled by root two.
Eight candidates, and the drop must be divided by the distance — otherwise diagonals always win.

Code examples

Example 1 — the full chain with both diagnoses reported 

import math

import numpy as np
import rasterio
from scipy.ndimage import label, minimum_filter
from skimage.morphology import reconstruction

OFFSETS = [(-1, 0), (-1, 1), (0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1)]
CODES = [64, 128, 1, 2, 4, 8, 16, 32]


def drainage_counts(surface):
    neighbour_min = minimum_filter(surface, size=3, mode="nearest")
    interior = np.zeros(surface.shape, bool)
    interior[1:-1, 1:-1] = True
    strict = (surface < neighbour_min) & interior
    stuck = (surface <= neighbour_min) & interior
    return {"pits": int(strict.sum()), "pit_regions": int(label(strict)[1]),
            "stuck": int(stuck.sum()), "stuck_regions": int(label(stuck)[1])}


def hydrology(path, *, stream_cells=200):
    with rasterio.open(path) as src:
        dem = src.read(1).astype("float64")
        if src.nodata is not None:
            dem = np.where(dem == src.nodata, np.nan, dem)
        if src.crs.is_geographic:
            lat = (src.bounds.bottom + src.bounds.top) / 2
            cell_x = abs(src.transform.a) * 111_320 * math.cos(math.radians(lat))
            cell_y = abs(src.transform.e) * 110_574
        else:
            cell_x, cell_y = abs(src.transform.a), abs(src.transform.e)
        profile = src.profile.copy()

    if np.isnan(dem).any():
        raise ValueError("fill voids before hydrology — NaN breaks flow routing")

    before = drainage_counts(dem)
    print(f"  before: {before['pits']:,} true pits, "
          f"{before['stuck']:,} cells with no downhill neighbour "
          f"({before['stuck_regions']:,} regions)")

    seed = np.full_like(dem, dem.max())
    seed[0, :], seed[-1, :] = dem[0, :], dem[-1, :]
    seed[:, 0], seed[:, -1] = dem[:, 0], dem[:, -1]
    filled = reconstruction(seed, dem, method="erosion")

    after = drainage_counts(filled)
    if after["pits"]:
        raise RuntimeError(f"{after['pits']} true pits survived filling")
    raised = filled - dem
    print(f"  filled: {(raised > 0).sum():,} cells raised, "
          f"median {np.median(raised[raised > 0]):.2f} m, max {raised.max():.2f} m")
    print(f"  after:  {after['pits']:,} true pits, {after['stuck']:,} on flats "
          f"({after['stuck_regions']:,} regions) — filling converts pits into flats")

    best_drop = np.full(filled.shape, -np.inf)
    direction = np.zeros(filled.shape, dtype="uint8")
    for (dr, dc), code in zip(OFFSETS, CODES):
        shifted = np.roll(np.roll(filled, -dr, axis=0), -dc, axis=1)
        drop = (filled - shifted) / np.hypot(dr * cell_y, dc * cell_x)
        better = drop > best_drop
        best_drop = np.where(better, drop, best_drop)
        direction = np.where(better, code, direction).astype("uint8")
    direction[best_drop <= 0] = 0
    direction[0, :] = direction[-1, :] = direction[:, 0] = direction[:, -1] = 0
    print(f"  {(direction == 0).sum() - 2 * (sum(filled.shape) - 2) - 4:,} interior cells "
          f"have no assigned direction (flats)")

    rows, cols = filled.shape
    order = np.argsort(filled.ravel())[::-1]
    accumulation = np.ones(filled.size)
    lookup = dict(zip(CODES, OFFSETS))
    flat_dir = direction.ravel()
    for idx in order:
        code = flat_dir[idx]
        if code == 0:
            continue
        dr, dc = lookup[code]
        r, c = divmod(int(idx), cols)
        nr, nc = r + dr, c + dc
        if 0 <= nr < rows and 0 <= nc < cols:
            accumulation[nr * cols + nc] += accumulation[idx]
    accumulation = accumulation.reshape(filled.shape)

    cell_km2 = cell_x * cell_y / 1e6
    print(f"  accumulation: max {accumulation.max():,.0f} cells "
          f"({accumulation.max() * cell_km2:.1f} km²)")
    return filled, direction, accumulation, (cell_x, cell_y), profile


filled, direction, accumulation, cell, profile = hydrology("snowdonia_glo30.tif")
  before: 0 true pits, 1,393 cells with no downhill neighbour (130 regions)
  filled: 2,057 cells raised, median 1.07 m, max 10.51 m
  after:  0 true pits, 2,067 on flats (114 regions) — filling converts pits into flats
  accumulation: max 3,866 cells (3.3 km²)

That maximum accumulation is the largest single basin the flow graph found: 3,866 cells, or 3.3 km². On an 8 × 7 km extent that is small, and the reason is exactly the flats — 114 regions where flow has no assigned direction, each truncating the paths that reach it.

This is the number to watch. A conditioned DEM with unresolved flats produces many small basins; a properly flat-resolved one produces a few large ones.

Example 2 — choosing the stream threshold 

import pandas as pd
from scipy.ndimage import label as connected


def stream_threshold_sweep(accumulation, cell, thresholds=(50, 100, 200, 500, 1000, 5000)):
    cell_km2 = cell[0] * cell[1] / 1e6
    rows = []
    for t in thresholds:
        streams = accumulation >= t
        rows.append({
            "cells": t,
            "area_km2": round(t * cell_km2, 3),
            "share": f"{streams.mean():.2%}",
            "segments": int(connected(streams)[1]),
        })
    print(pd.DataFrame(rows).to_string(index=False))


stream_threshold_sweep(accumulation, cell)
 cells  area_km2 share  segments
    50     0.043 6.53%      1312
   100     0.086 3.47%       780
   200     0.171 1.82%       412
   500     0.428 0.77%       205
  1000     0.857 0.37%       100
  5000     4.283 0.00%         0

Read the segments column as a warning rather than a result. A real drainage network over this area is a handful of connected systems, not 412 fragments — and at 5,000 cells nothing survives at all, because no basin in the flow graph is that large.

That pattern, many segments at every threshold and nothing above a moderate one, is the signature of unresolved flats. The area_km2 column is still the right way to express the threshold: "a channel begins where 0.17 km² drains into it" is defensible; "200 cells" is meaningless at a different resolution.

Example 3 — delineating a watershed, and checking it 

def watershed(direction, outlet):
    """Walk the flow graph upstream from one cell."""
    rows, cols = direction.shape
    basin = np.zeros(direction.shape, bool)
    basin[outlet] = True
    stack = [outlet]

    while stack:
        r, c = stack.pop()
        for (dr, dc), code in zip(OFFSETS, CODES):
            nr, nc = r - dr, c - dc                   # cells that drain INTO (r, c)
            if not (0 <= nr < rows and 0 <= nc < cols) or basin[nr, nc]:
                continue
            if direction[nr, nc] == code:
                basin[nr, nc] = True
                stack.append((nr, nc))
    return basin


outlet = np.unravel_index(accumulation.argmax(), accumulation.shape)
basin = watershed(direction, outlet)

cell_km2 = cell[0] * cell[1] / 1e6
print(f"outlet {outlet}, accumulation {accumulation[outlet]:,.0f}")
print(f"basin {basin.sum():,} cells = {basin.sum() * cell_km2:.2f} km²")
print(f"consistent: {abs(basin.sum() - accumulation[outlet]) < 1}")
outlet (216, 124), accumulation 3,866
basin 3,866 cells = 3.31 km²
consistent: True

The last line is a real check, not decoration. Flow accumulation at a cell is by definition the number of cells draining through it, so an independent upstream walk must return exactly that count.

If the two disagree, the direction grid contains a cycle or an inconsistency — which is what an incomplete fill, or a direction assigned into a flat, produces. It is worth asserting.

Explanation

Why filling and flat-resolution are separate problems

Filling guarantees no cell is strictly lower than all its neighbours. It says nothing about cells that are equal to a neighbour.

On a filled depression every cell is at the outlet level — a plateau. D8 asks "which neighbour gives the steepest descent" and the answer is none, because the descent is zero in every direction. The cell gets direction 0 and flow stops.

The measured effect on this DEM: filling raised 2,057 cells, and the count of cells with nowhere to drain rose from 1,393 to 2,067. The problem did not go away; it changed shape.

Proper flat-resolution imposes a small artificial gradient across each flat — away from higher ground and toward the outlet — so that D8 has something to follow. That is the step hand-rolled implementations usually omit, and it is why pysheds, richdem and whitebox produce connected networks where a fifty-line NumPy version produces fragments.

Why zero true pits is a clue about your DEM

A raw photogrammetric or radar DEM is noisy and has thousands of genuine single-cell pits. Copernicus DEM has none, because it is a conditioned product — voids filled, noise smoothed, and in places already hydrologically enforced.

So the strict test returning zero is informative rather than reassuring. It tells you the DEM has been processed, that filling will change little, and that flats are where your effort should go.

Run both tests. The strict one characterises the DEM; the <= one tells you what will actually break.

A depression before filling with one pit cell, and after filling as a plateau of six cells with no downhill neighbour.
One pit becomes six stuck cells. Filling trades a small hard problem for a larger easy-looking one.

Why D8 is crude and still standard

D8 sends all of a cell's flow to one neighbour. Real water spreads, especially on hillslopes and fans, and D8 cannot represent that — a slope draining south-south-east becomes a zigzag of south and south-east steps.

D-infinity apportions flow between two neighbours by angle; MFD spreads it among all downhill neighbours. Both give smoother, more realistic accumulation.

D8 survives because it gives each cell exactly one downstream neighbour, which makes the flow graph a forest. That is what allows the single-pass accumulation in Example 1 and the simple upstream walk in Example 3. With MFD neither works, and watershed delineation needs a fundamentally different algorithm.

Why the diagonal correction matters

Comparing raw height drops rather than gradients biases every decision toward diagonals: a diagonal neighbour is √2 further away, so it can drop 41% more without being steeper.

On a uniform planar slope this produces alternating diagonal and orthogonal steps in a regular pattern — a herringbone across the whole hillside, visible immediately in the direction raster and inherited by every stream. Dividing by hypot(dr * cell_y, dc * cell_x) fixes it, and also handles the non-square cells that a geographic DEM produces.

Edge cases or notes

  • Run both pit tests. Strictly-lower finds true pits; less-than-or-equal finds everything that cannot drain, which is what actually breaks flow.
  • Filling increases the flat area. Expect the "no downhill neighbour" count to go up, not down.
  • Never use the filled DEM for elevation, slope or contours. It is a hydrological working surface.
  • Real closed basins are destroyed by filling. Many single-cell depressions are noise; a few large ones may be real.
  • NaN breaks flow routing entirely. Fill voids first.
  • np.roll wraps at the edges. Mask the boundary rows and columns.
  • Divide diagonal drops by the true distance, or diagonals win and you get a herringbone.
  • pysheds, richdem and whitebox implement flat-resolution and D-infinity properly. For real hydrology, use them — this article is for understanding what they do.

FAQ

What is the difference between a pit and a flat?

A pit is strictly lower than all eight neighbours. A flat has at least one neighbour at exactly the same height. Filling removes pits and creates flats.

My DEM has zero pits. Do I still need to fill it?

Probably not, but check the other test. Copernicus DEM has zero true pits and 1,393 cells with no downhill neighbour — all of them flats, which filling will not help with.

Why did filling make the problem worse?

It did not make it worse, it changed its shape. Every filled depression becomes a plateau at the outlet level, so the count of cells with nowhere to drain rises while the number of separate problem regions falls.

Why is my stream network fragmented?

Unresolved flats. Each one truncates the flow paths reaching it. Use a library that implements flat-resolution — pysheds, richdem or whitebox.

What stream threshold should I use?

Express it as a contributing area in km² rather than a cell count, and calibrate against a mapped network. A cell count is meaningless at a different resolution.

Why divide diagonal drops by √2?

Because diagonal neighbours are further away. Comparing raw drops makes diagonals win 41% of the time they should not, producing a herringbone flow pattern.

Can I use the filled DEM for anything else?

No. Keep the original for slope, contours and any elevation reporting. The filled surface exists only to make flow routing possible.