⚙️ How-To Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 14 of 20 Updated

How to Calculate Height Above Nearest Drainage (HAND) in Python

Problem statement

Height Above Nearest Drainage gives every cell of a DEM its elevation above the stream cell it drains to. Cells a metre or two above their stream are candidates for flooding; cells 50 m above it are not. It needs only a DEM, it runs in seconds, and it is widely used for screening flood exposure where no flood model exists. The calculation is short, but three choices change the answer: the surface it is measured on, the stream network it is measured to, and the cut-off that turns heights into a flood-prone zone.

Calculated for the Esopus Creek basin above Coldbrook, New York, on the 10 m 3DEP DEM, and checked against the FEMA special flood hazard area (zones A and AE, 14.61 km² inside the basin):

  • WhiteboxTools computed HAND for 11.4 million cells in about a second per stream network, after a 1.8 s fill, pointer and accumulation.
  • pysheds' compute_hand agreed with WhiteboxTools to within 0.5 m on 98.6% of cells, with a median difference of 0.0 m — but measured on the raw DEM rather than the conditioned surface, it produced 1,734 negative heights down to −6.76 m.
  • With streams defined at 10 km² and a cut-off of 5 m, the HAND zone matched the FEMA area best: critical success index 0.578, 73.8% of the FEMA area captured, 14.8 km² predicted against 14.6 km².
  • The same cut-off on streams at 0.1 km² predicted a zone more than twice the size, because every small valley floor became "near drainage".

Quick answer 

import whitebox

wbt = whitebox.WhiteboxTools()
wbt.fill_depressions("/data/dem.tif", "/data/filled.tif", fix_flats=True)
wbt.d8_pointer("/data/filled.tif", "/data/d8.tif")
wbt.d8_flow_accumulation("/data/d8.tif", "/data/acc.tif", out_type="cells", pntr=True)
wbt.extract_streams("/data/acc.tif", "/data/streams.tif", threshold=100_000)      # 10 km2 of 10 m cells
wbt.elevation_above_stream("/data/filled.tif", "/data/streams.tif", "/data/hand.tif")

Condition the DEM, choose a stream threshold that matches the rivers you care about, and compute HAND on the same conditioned surface used for routing.

Flow of a HAND calculation: fill the DEM, compute D8 flow directions and accumulation, extract streams at a threshold, then measure each cell's height above the stream cell it drains to.
Five steps, and the stream threshold in step four matters more than any of the others.

Step-by-step solution

1. Use a projected, bare-earth DEM

HAND measures vertical differences along flow paths, so the DEM's surface matters: tree canopy and buildings in a surface model add metres of false height. The 3DEP 10 m DEM here is bare earth in UTM zone 18N. See how DEM resolution and source change drainage.

2. Condition the DEM

Pits and flats break flow paths, leaving cells that drain nowhere and have no HAND. WhiteboxTools' fill_depressions with fix_flats=True took 1.8 s together with the pointer and accumulation; see breaching or filling.

3. Route flow and accumulate

HAND follows D8 flow directions from each cell to the first stream cell downstream. Compute the pointer and the accumulation in cells on the conditioned DEM.

4. Choose the stream threshold deliberately

The stream network defines "drainage". At 10 m cells, a threshold of 100,000 cells is 10 km². Lower thresholds reach up small side valleys and make their floors low-HAND; higher ones leave floodplains of mid-sized tributaries measured against the main river. Against the FEMA area, 10 km² worked best here (Example 3).

5. Compute HAND on the conditioned surface

elevation_above_stream(filled, streams, hand) subtracts the stream cell's elevation from each cell's elevation along the flow path. Use the same filled DEM you routed on. Measuring the raw DEM against flow paths from a conditioned one gives negative HAND wherever filling raised a cell: 1,734 cells, as low as −6.76 m, in pysheds.

6. Check the output

HAND should be zero on the streams, non-negative everywhere, and NoData only where the DEM is NoData or flow leaves the grid before reaching a stream. The share of cells within a few metres of drainage is a useful sanity check: most of a mountain basin should not be.

7. Calibrate a cut-off against a reference

A HAND value is a height, not a flood extent. To turn it into a zone, compare HAND ≤ t with mapped flood areas for several t and choose by a score that penalises both misses and false alarms, such as the critical success index (CSI = hits / (hits + misses + false alarms)).

8. Save the zone with its parameters

Record the DEM, conditioning, stream threshold and cut-off with the output. "HAND ≤ 5 m from 10 km² streams on 3DEP 10 m" is reproducible; "HAND flood zone" is not.

Bar chart of the best critical success index against the FEMA flood hazard area for HAND zones built from streams at 0.1, 1, 10 and 50 square kilometres on the 3DEP 10 m DEM.
The stream threshold decided the match more than the HAND cut-off did.

Code examples

Example 1 — HAND with WhiteboxTools 

import os

import numpy as np
import rasterio
import whitebox

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", out_type="cells", pntr=True)
wbt.extract_streams("acc.tif", "streams_10km2.tif", threshold=100_000)          # 10 km2 of 10 m cells
wbt.elevation_above_stream("filled.tif", "streams_10km2.tif", "hand_10km2.tif")

with rasterio.open("hand_10km2.tif") as src:
    hand = src.read(1, masked=True)
print(f"{hand.count():,} cells; HAND {hand.min():.2f} to {hand.max():.2f} m, median {np.ma.median(hand):.2f} m; "
      f"within 5 m of drainage {(hand <= 5).sum() / hand.count():.1%}")
11,427,564 cells; HAND 0.00 to 882.97 m, median 164.50 m; within 5 m of drainage 5.0%

Example 2 — pysheds, and the raw-DEM trap 

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 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))
fdir = grid.flowdir(surface)
streams = grid.accumulation(fdir) > 100_000

for label, elevation in (("raw DEM", dem), ("conditioned surface", surface)):
    values = np.asarray(grid.compute_hand(fdir, elevation, streams))
    valid = np.isfinite(values) & (np.asarray(dem) != dem.nodata)
    print(f"pysheds HAND on the {label}: {int((values[valid] < 0).sum()):,} negative cells, minimum {values[valid].min():.2f} m")

pysheds_hand = np.asarray(grid.compute_hand(fdir, surface, streams))
both = np.isfinite(pysheds_hand) & ~hand.mask
difference = np.abs(pysheds_hand[both] - hand.data[both])
print(f"against WhiteboxTools: {both.sum():,} cells, median |difference| {np.median(difference):.2f} m, "
      f"90th percentile {np.percentile(difference, 90):.2f} m, within 0.5 m {np.mean(difference <= 0.5):.1%}")
pysheds HAND on the raw DEM: 1,734 negative cells, minimum -6.76 m
pysheds HAND on the conditioned surface: 0 negative cells, minimum 0.00 m
against WhiteboxTools: 10,900,023 cells, median |difference| 0.00 m, 90th percentile 0.00 m, within 0.5 m 98.6%

Example 3 — calibrate a flood-prone zone against FEMA 

import geopandas as gpd
from rasterio.features import rasterize

with rasterio.open("hand_10km2.tif") as src:
    profile, transform, shape = src.profile, src.transform, src.shape
cell_km2 = abs(transform.a * transform.e) / 1e6
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(profile["crs"])
fema = gpd.read_file("fema_flood_hazard_esopus.gpkg").to_crs(profile["crs"])
inside = rasterize(basin.geometry, out_shape=shape, transform=transform).astype(bool)
sfha = rasterize(fema.loc[fema.SFHA_TF == "T"].geometry.make_valid(), out_shape=shape, transform=transform).astype(bool) & inside

scores = []
for cutoff in (1, 2, 3, 4, 5, 6, 8, 10):
    predicted = inside & ~hand.mask & (hand.data <= cutoff)
    hits, false_alarms, misses = (predicted & sfha).sum(), (predicted & ~sfha).sum(), (sfha & ~predicted).sum()
    scores.append((cutoff, hits / (hits + false_alarms + misses), hits / sfha.sum(), predicted.sum() * cell_km2))
    print(f"HAND <= {cutoff:>2} m: CSI {scores[-1][1]:.3f}, FEMA area captured {scores[-1][2]:.1%}, zone {scores[-1][3]:.1f} km2")

best = max(scores, key=lambda s: s[1])[0]
zone = (inside & ~hand.mask & (hand.data <= best)).astype("uint8")
profile.update(dtype="uint8", nodata=0, compress="deflate")
with rasterio.open(f"hand_le_{best}m_streams_10km2.tif", "w", **profile) as dst:
    dst.write(zone, 1)
    dst.update_tags(dem="3DEP 1/3 arc-second, UTM 18N", conditioning="fill_depressions fix_flats",
                    streams="D8 accumulation >= 10 km2", cutoff_m=str(best))
print(f"FEMA area {sfha.sum() * cell_km2:.2f} km2; saved HAND <= {best} m zone, {zone.sum() * cell_km2:.1f} km2")
HAND <=  1 m: CSI 0.250, FEMA area captured 25.7%, zone 4.2 km2
HAND <=  2 m: CSI 0.393, FEMA area captured 41.7%, zone 7.0 km2
HAND <=  3 m: CSI 0.494, FEMA area captured 55.0%, zone 9.7 km2
HAND <=  4 m: CSI 0.554, FEMA area captured 65.9%, zone 12.4 km2
HAND <=  5 m: CSI 0.578, FEMA area captured 73.8%, zone 14.8 km2
HAND <=  6 m: CSI 0.577, FEMA area captured 79.9%, zone 17.3 km2
HAND <=  8 m: CSI 0.542, FEMA area captured 87.1%, zone 21.6 km2
HAND <= 10 m: CSI 0.496, FEMA area captured 90.8%, zone 25.4 km2
FEMA area 14.61 km2; saved HAND <= 5 m zone, 14.8 km2

Explanation

What HAND measures

Every cell drains, along its D8 flow path, to a first stream cell. HAND is the elevation difference between the two. It normalises terrain by the local drainage: a valley floor 800 m above sea level and one at 200 m both have HAND near zero, while the ridge between them has HAND of hundreds of metres. That is why a single cut-off can sketch floodplains across a whole basin.

Why the conditioned surface is needed

Flow paths come from the conditioned DEM, which raised pits to their spill level. On the raw DEM, a cell at the bottom of a filled pit can be lower than the stream cell its path reaches, giving a negative height. Using the conditioned surface keeps the arithmetic consistent with the paths; the pits themselves have HAND equal to their fill level above the stream.

Why the stream threshold dominates

HAND is relative to the nearest drainage, so the network defines what counts as the reference level. A dense network makes many small valley floors low; a sparse one measures tributary floodplains against a distant main river and makes them high. The FEMA area in the Esopus basin follows streams of roughly 10 km² and larger, which is why that threshold scored best.

Why WhiteboxTools and pysheds agree

Both follow D8 paths to the first stream cell and subtract elevations. Differences come from how each resolves flats and which cells each treats as streams, which at 10 m affected a small share of cells: 98.6% agreed to within 0.5 m.

Two panels contrasting HAND computed on the raw DEM, with negative values in filled pits, and on the conditioned surface, with none.
Negative HAND is a mismatch between the surface measured and the surface routed, not a real depression below the river.

Edge cases or notes

  • Streams leaving the grid leave upstream cells without a stream to reach; extend the DEM beyond the basin.
  • Lakes and reservoirs are flat and near zero HAND; mask them if the zone should exclude open water.
  • Culverts and embankments block flow paths and raise HAND upstream of roads; see streams blocked by roads.
  • Surface models such as GLO-30 put canopy into HAND and scored lower against FEMA at every threshold.
  • Large rivers can flood well above a fixed cut-off; vary the cut-off with stream size for regional work.
  • HAND is not depth: a zone from HAND ≤ 5 m has no return period and no hydraulics.
  • Output data type: store HAND as float32 with a NoData value; many tools write float64 by default.

FAQ

How do I calculate HAND in Python?

Condition the DEM, compute D8 flow directions and accumulation, extract streams at a threshold, then run WhiteboxTools' elevation_above_stream or pysheds' compute_hand. On the Esopus 10 m DEM each HAND run took about a second.

Why does my HAND raster have negative values?

The heights were measured on the raw DEM while flow paths came from a conditioned surface. pysheds produced 1,734 negative cells that way and none on the conditioned surface.

What stream threshold should I use for HAND?

One that matches the rivers whose floodplains you want. Against FEMA's flood hazard area, 10 km² gave the best match on the Esopus basin; 0.1 km² predicted twice the area.

What HAND value indicates flood risk?

There is no universal value. Calibrated against FEMA zones A and AE, HAND ≤ 5 m from 10 km² streams scored best here, capturing 73.8% of the mapped area.

Is pysheds or WhiteboxTools better for HAND?

They agree closely: 98.6% of cells within 0.5 m. WhiteboxTools conditions faster; pysheds keeps everything in NumPy arrays.

Can HAND replace a flood model?

No. It ranks terrain by height above drainage and has no flows, return periods or hydraulics; use it for screening and prioritising.