Height Above Nearest Drainage Explained: A Quick Flood Proxy
On this page
Problem statement
Flood maps come from hydraulic models: flows, cross-sections, roughness, calibration. Most of the world's streams have none. Height Above Nearest Drainage (HAND) offers a shortcut: from a DEM alone, give each cell its height above the stream it drains to, and treat low values as flood-prone. It is fast, repeatable and surprisingly informative. It is also easy to over-read: HAND has no discharge, no return period and no water surface, and its answer depends on two choices that have nothing to do with floods — the stream network and the DEM.
Tested against the FEMA special flood hazard area (zones A and AE: 14.61 km², as computing HAND measures it) in the Esopus Creek basin above Coldbrook, New York, with HAND from WhiteboxTools:
- The best match on the 10 m 3DEP DEM came from streams defined at 10 km² and a cut-off of 6 m: critical success index 0.577, 79.6% of the FEMA area captured, 17.2 km² predicted.
- Streams at 0.1 km² made nearly every small valley floor low: best CSI 0.335, with 32.6 km² predicted to capture 81% of the FEMA area.
- Streams at 50 km² measured tributary floodplains against the main river: the median HAND inside the FEMA area rose from 2.65 m to 10.23 m, and the best CSI fell to 0.398.
- On the 30 m GLO-30 surface model every threshold scored lower — best 0.432, needing an 8 m cut-off and 22.2 km² of predicted area.
Quick answer
HAND is the vertical distance from a cell down to the stream cell its flow reaches:
HAND(cell) = elevation(cell) − elevation(first stream cell downstream of it)
Low HAND means a cell sits little above the channel that drains it. It is a good screening layer for where flooding is plausible, a poor substitute for a flood map, and meaningless without the stream threshold and DEM it was computed from. Calibrate any cut-off against mapped floods before using it.
What HAND normalises, and why the threshold dominates
1. Understand what HAND normalises
Raw elevation cannot separate floodplains from terraces across a basin: a valley floor at 800 m and one at 200 m look nothing alike. HAND subtracts the local drainage level, so both valley floors are near zero and the ridges between them are hundreds of metres. On the Esopus DEM the median HAND was 166.14 m (computing HAND measures it); only 5.0% of cells were within 5 m of their stream.
2. See where the reference level comes from
Each cell is compared with the first stream cell on its D8 flow path. That makes HAND a property of the flow directions and the network together: change either and HAND changes, even with the same DEM.
3. Expect the stream threshold to dominate
The network is defined by an accumulation threshold. On the Esopus 3DEP DEM, against the FEMA area:
streams from median HAND in FEMA area 90th pct best cut-off CSI FEMA captured predicted
0.1 km2 1.20 m 4.06 m 3 m 0.335 81.1% 32.6 km2
1 km2 1.95 m 6.74 m 4 m 0.415 75.0% 22.7 km2
10 km2 2.65 m 9.65 m 6 m 0.577 79.6% 17.2 km2
50 km2 10.23 m 79.80 m 8 m 0.398 46.3% 9.2 km2
Low thresholds flag every hollow; high ones leave mapped floodplains on tributaries measured against a distant river (Example 1).
4. Match the threshold to the rivers the flood map covers
FEMA studies rivers and larger streams; a 10 km² network resembles what it mapped. For a different reference — pluvial flooding, small-stream flash floods — a different threshold will match better. There is no universal value.
5. Prefer a bare-earth DEM
A surface model includes canopy and buildings, which add height to cells on valley floors and distort channels. GLO-30 scored lower than 3DEP at every threshold: CSI 0.217, 0.313, 0.432 and 0.327 at 0.1, 1, 10 and 50 km². Its median HAND inside the FEMA area at 10 km² was 2.00 m against 2.65 m, but its 90th percentile was 14.93 m against 9.65 m: canopy and coarse cells spread the floodplain's heights out (Example 2).
6. Calibrate the cut-off, don't borrow it
With 10 km² streams on 3DEP, CSI rose from 0.248 at 1 m to 0.573 at 5 m and 0.577 at 6 m, then fell to 0.495 at 10 m. Published cut-offs from other regions reflect their rivers, relief and DEMs; calibrate against local flood maps or observed flood extents.
7. Read HAND by flood zone, not only by threshold
Mapped flood zones differ in what they represent: detailed AE studies, approximate A zones and the 0.2%-annual-chance zone. Summarising HAND inside each (Example 3) shows whether a single cut-off can represent all of them. With 10 km² streams, AE zones had a median HAND of 2.93 m and the 0.2% zone 6.57 m, so the rarer flood sits visibly higher. Approximate A zones had a median of 1.21 m but a 90th percentile of 80.38 m: a tenth of their area stood more than 80 m above the nearest 10 km² stream, which no single cut-off can include without flooding hillsides. Outside every zone, only 0.6% of the basin was within 5 m of drainage.
8. Use HAND for what it is good at
Screening where detailed studies are missing, ranking assets by height above their stream, stratifying samples, and adding a terrain covariate to models. Not for insurance zones, design flood levels or anything needing a return period.
HAND against FEMA flood zones
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); Copernicus GLO-30 resampled the same way at 30 m over the basin plus 4 km (esopus_glo30_utm.tif); FEMA's National Flood Hazard Layer zones, from Esri's hosted copy.
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
dem_box("glo30", 4000, "esopus_glo30_utm.tif") # the same, basin + 4 km
if not os.path.exists("fema_flood_hazard_esopus.gpkg"): # FEMA NFHL flood zones, Esri's copy
feats, off = [], 0
while True:
page = get("https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/USA_Flood_Hazard_Reduced_Set_gdb/FeatureServer/0/query",
params={"geometry": f"{W},{S},{E},{N}", "geometryType": "esriGeometryEnvelope", "inSR": 4326,
"spatialRel": "esriSpatialRelIntersects", "outFields": "*", "outSR": 4326, "f": "geojson",
"resultOffset": off, "resultRecordCount": 1000}).json()
feats += page["features"]
# the server caps a page at its own maximum, not at what was asked for: follow the flag
if not page["features"] or not (page.get("exceededTransferLimit") or page.get("properties", {}).get("exceededTransferLimit")):
break
off += len(page["features"])
gpd.GeoDataFrame.from_features(feats, crs=4326).to_file("fema_flood_hazard_esopus.gpkg")
for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm.tif', 'esopus_glo30_utm.tif', 'fema_flood_hazard_esopus.gpkg']:
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_glo30_utm.tif: 1,357 x 1,003 cells of 30 m
fema_flood_hazard_esopus.gpkg: 3,263 features
Example 1 — HAND at four stream thresholds against FEMA zones
import os
import geopandas as gpd
import numpy as np
import rasterio
import whitebox
from rasterio.features import rasterize
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("."))
basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(32618)
fema = gpd.read_file("fema_flood_hazard_esopus.gpkg").to_crs(32618)
fema["geometry"] = fema.geometry.make_valid()
def hand_rasters(dem, tag, thresholds_km2):
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(dem) as src:
transform, shape = src.transform, src.shape
hands = {}
for km2 in thresholds_km2:
wbt.extract_streams(f"{tag}_acc.tif", f"{tag}_streams.tif", threshold=km2 * 1e6 / abs(transform.a * transform.e))
wbt.elevation_above_stream(f"{tag}_filled.tif", f"{tag}_streams.tif", f"{tag}_hand.tif")
with rasterio.open(f"{tag}_hand.tif") as src:
hands[km2] = src.read(1, masked=True)
return hands, transform, shape
def report(hands, transform, shape):
cell_km2 = abs(transform.a * transform.e) / 1e6
inside = rasterize(basin.geometry, out_shape=shape, transform=transform).astype(bool)
sfha = rasterize(fema.loc[fema.SFHA_TF == "T"].geometry, out_shape=shape, transform=transform).astype(bool) & inside
for km2, hand in hands.items():
in_zone = hand.data[sfha & ~hand.mask]
scores = []
for cutoff in (1, 2, 3, 4, 5, 6, 8, 10):
predicted = inside & ~hand.mask & (hand.data <= cutoff)
hits = (predicted & sfha).sum()
scores.append((hits / (predicted.sum() + (sfha & ~predicted).sum()), cutoff, hits / sfha.sum(), predicted.sum() * cell_km2))
csi, cutoff, captured, area = max(scores)
print(f"streams >= {km2:>4} km2: HAND in FEMA area median {np.median(in_zone):5.2f} m, 90th {np.percentile(in_zone, 90):5.2f} m; "
f"best CSI {csi:.3f} at <= {cutoff} m ({captured:.1%} captured, {area:.1f} km2)")
return inside
hand_3dep, transform_3dep, shape_3dep = hand_rasters("esopus_3dep13_utm.tif", "dep", (0.1, 1, 10, 50))
inside_3dep = report(hand_3dep, transform_3dep, shape_3dep)streams >= 0.1 km2: HAND in FEMA area median 1.20 m, 90th 4.06 m; best CSI 0.335 at <= 3 m (81.1% captured, 32.6 km2)
streams >= 1 km2: HAND in FEMA area median 1.95 m, 90th 6.74 m; best CSI 0.415 at <= 4 m (75.0% captured, 22.7 km2)
streams >= 10 km2: HAND in FEMA area median 2.65 m, 90th 9.65 m; best CSI 0.577 at <= 6 m (79.6% captured, 17.2 km2)
streams >= 50 km2: HAND in FEMA area median 10.23 m, 90th 79.80 m; best CSI 0.398 at <= 8 m (46.3% captured, 9.2 km2)
Example 2 — the same test on a 30 m surface model
hand_glo, transform_glo, shape_glo = hand_rasters("esopus_glo30_utm.tif", "glo", (0.1, 1, 10, 50))
report(hand_glo, transform_glo, shape_glo)streams >= 0.1 km2: HAND in FEMA area median 0.45 m, 90th 9.86 m; best CSI 0.217 at <= 4 m (70.0% captured, 43.0 km2)
streams >= 1 km2: HAND in FEMA area median 1.20 m, 90th 11.37 m; best CSI 0.313 at <= 4 m (65.6% captured, 25.7 km2)
streams >= 10 km2: HAND in FEMA area median 2.00 m, 90th 14.93 m; best CSI 0.432 at <= 8 m (75.8% captured, 22.2 km2)
streams >= 50 km2: HAND in FEMA area median 15.14 m, 90th 83.19 m; best CSI 0.327 at <= 10 m (42.2% captured, 10.4 km2)
Example 3 — HAND inside each kind of flood zone
hand = hand_3dep[10]
zones = {
"AE (detailed study)": fema.FLD_ZONE == "AE",
"A (approximate)": fema.FLD_ZONE == "A",
"X, 0.2% annual chance": fema.ZONE_SUBTY == "0.2 Percent Annual Chance Flood Hazard",
}
any_zone = np.zeros(shape_3dep, dtype=bool)
for name, rows in zones.items():
mask = rasterize(fema.loc[rows].geometry, out_shape=shape_3dep, transform=transform_3dep).astype(bool) & inside_3dep
any_zone |= mask
values = hand.data[mask & ~hand.mask]
print(f"{name:24} {mask.sum() / 1e4:6.2f} km2 HAND median {np.median(values):5.2f} m, 90th {np.percentile(values, 90):6.2f} m")
outside = hand.data[inside_3dep & ~any_zone & ~hand.mask]
print(f"{'outside all zones':24} {outside.size / 1e4:6.2f} km2 HAND median {np.median(outside):5.2f} m, share <= 5 m {np.mean(outside <= 5):.1%}")AE (detailed study) 12.23 km2 HAND median 2.93 m, 90th 8.63 m
A (approximate) 2.38 km2 HAND median 1.21 m, 90th 80.38 m
X, 0.2% annual chance 3.35 km2 HAND median 6.57 m, 90th 12.17 m
outside all zones 477.16 km2 HAND median 216.21 m, share <= 5 m 0.6%
A flood proxy, and only a proxy
Why HAND works as a flood proxy at all
Flood water rises from the channel. For a given flood, the cells that get wet are, to a first approximation, those less than some height above the water in the channel next to them. HAND measures exactly that height from the channel bed, so a single cut-off sketches the kind of zone a flood of a certain depth would make — if every river rose by the same amount.
Why it is only a proxy
Rivers do not all rise by the same amount: a large river's flood stage is metres higher than a small stream's. Backwater from confluences, levees, bridges and channel capacity all shape real flooding and none are in HAND. And the D8 path to "nearest drainage" can cross a divide into a different stream where the DEM smooths a ridge, which misattributes the reference level.
Why the stream threshold matters so much
HAND is always relative to the network. At 0.1 km² the network reaches nearly every hollow, so hillside valley floors have near-zero HAND and flood-prone area is overpredicted. At 50 km² tributaries are not streams, so their floodplains are measured against the main river kilometres away, and their HAND is tens of metres — the 90th percentile inside the FEMA area reached 79.80 m.
Why surface models score lower
GLO-30 is a digital surface model: its elevations include forest canopy over most of this basin. Canopy on valley floors adds height, canopy at the channel's edge moves the apparent channel, and 30 m cells blur narrow floodplains into their valley walls. Each adds error to HAND in exactly the cells that matter.
Coasts, deltas, cities and dams
- Coastal and tidal flooding is controlled by sea level, not river level; HAND does not apply.
- Flat deltas and plains have low HAND everywhere; cut-offs lose discrimination.
- Urban drainage follows pipes and streets; DEM streams and HAND miss it.
- Dams and reservoirs create flat water surfaces with zero HAND upstream.
- Variable cut-offs scaled by upstream area represent larger rivers' deeper floods better than one value.
- Validation data should be independent of the calibration data; FEMA zones used to choose a cut-off cannot also confirm it.
- Negative HAND means heights were measured on a different surface from the one routed; see calculating HAND in Python.
Internal links
- How to calculate height above nearest drainage (HAND) in Python — the workflow and the pysheds trap
- How to extract a stream network from a DEM — thresholds and networks
- How DEM resolution and source change a drainage network — why GLO-30 differs
- DEM hydrology explained: from elevation to where water goes — the routing HAND relies on
- LiDAR surfaces explained: DSM, DTM and CHM — surface versus bare-earth models
- Flow direction explained: D8, D-Infinity and MFD compared — the paths to nearest drainage
- How to burn a known river network into a DEM — aligning channels before HAND
- How to rasterize a vector layer in Python — flood zones as masks
- A SAR flood map shows water where there is none — HAND as a mask for radar flood maps
- Flood hazard maps explained — what a modelled flood map adds that a terrain proxy cannot
FAQ
What is Height Above Nearest Drainage?
The elevation of each DEM cell above the stream cell its flow path reaches. Cells with low HAND lie little above their local channel and are candidates for flooding.
Is HAND a flood map?
No. It has no discharge, return period or water surface. Calibrated against FEMA zones in the Esopus basin, the best HAND zone captured 79.6% of the mapped area with a CSI of 0.577.
What HAND threshold indicates flooding?
It depends on the stream network, the DEM and the flood map. With 10 km² streams on a 10 m DEM, 5–6 m matched FEMA best here (CSI 0.573 and 0.577); calibrate locally.
Why does the stream network change HAND?
HAND is measured to the nearest stream. A dense network lowers HAND on small valley floors; a sparse one raises tributary floodplains to tens of metres above the main river.
Can I compute HAND from a 30 m global DEM?
Yes, but expect a weaker match. GLO-30, a surface model, scored lower than the 10 m bare-earth 3DEP DEM at every stream threshold tested.
What is HAND useful for?
Screening flood-prone land where no flood study exists, ranking sites by height above their stream, and as a terrain variable in statistical models.