🧭 Concept Hydrology & watersheds Terrain & elevation Raster & rasterio Hydrology & Watersheds · step 12 of 20 Updated

How DEM Resolution and Source Change a Drainage Network

Problem statement

The same watershed workflow run on a 1 m lidar DEM, the 10 m USGS 3DEP model and a 30 m global DEM gives three different catchments and three different stream networks. It is tempting to rank them by resolution — finer must be better — but two things change at once when you swap DEMs: the cell size and the source. A 30 m global surface model is not simply a blurred 10 m bare-earth model, and a 1 m lidar surface resolves features that route water somewhere unexpected.

Measured on the Esopus Creek basin above Coldbrook (published 497.3 km²) and its Birch Creek sub-basin at Big Indian (published 32.37 km²), New York, with the same WhiteboxTools fill, pointer and accumulation on every DEM:

  • Averaging the 10 m 3DEP DEM to 30 m barely changed anything: catchment 492.45 → 493.05 km², mean slope 17.01° → 16.72°, and 94.0% → 93.4% of major mapped streams found within 30 m by derived channels.
  • The 30 m Copernicus GLO-30 surface model on the same grid found only 70.5% — the source, not the cell size, cost 23 percentage points of alignment.
  • At Birch Creek the 1 m lidar DEM gave the worst catchment: 25.30 km², 21.8% short, IoU 0.75, against 32.47 km² at 10 m. The largest accumulation within the same 150 m snap radius lay on the edge of the search window: the snap stopped short of a larger channel.
  • A threshold of 1,000 cells meant 0.001 km² at 1 m and 0.9 km² at 30 m: about 1,269 km of channel in the Birch Creek basin against 19 km.

Quick answer 

what you compare              what changes                        what to do
10 m vs 10 m averaged to 30 m resolution only                     little changed; 30 m is fine for basin-scale areas
10 m bare earth vs 30 m DSM   source: canopy, smoothing, artefacts alignment dropped sharply; prefer bare earth
1 m lidar vs 10 m             detail: embankments, ditches, culverts expect local re-routing; condition and check outlets
cell thresholds across DEMs   area per threshold                  always express thresholds as areas

Choose the finest bare-earth DEM your extent and runtime allow, express thresholds in area, and validate outlets and networks on each DEM rather than assuming a finer one is right.

Bar chart of the share of major mapped streams found within 30 metres of derived channels for the 3DEP 10 m DEM, the same DEM averaged to 30 m, and the GLO-30 30 m surface model.
Coarsening the same bare-earth data cost under one point; switching to a 30 m surface model cost twenty-three.

Step-by-step solution

1. Separate resolution from source

Resample your best DEM to the coarse DEM's grid with an average, and run the workflow on both. Differences between the 10 m DEM and its 30 m average are resolution; differences between the 30 m average and the native 30 m product are source (Example 1).

2. Compare catchment areas

At the Coldbrook gauge the three Esopus DEMs gave 492.45, 493.05 and 490.46 km², all within 1.4% of the published 497.3 km², with intersection-over-union against the USGS polygon of 0.9853, 0.9865 and 0.9794. For a large catchment, area is the least sensitive product.

3. Compare channel positions

Derived channels of at least 1 km² found 94.0% of NHDPlus HR's major streams (order 3 and up) within 30 m on the 10 m DEM and 93.4% on its 30 m average, but 70.5% on GLO-30. A surface model includes forest canopy, and this basin is 98% tree cover, so the surface GLO-30 routes over is partly the top of the forest rather than the ground. Burning mapped streams can repair it — see burning a river network into a DEM.

4. Compare network length and slope

Channel length at 1 km² fell only from 336.7 km at 10 m to 323.1 km at 30 m, and 320.8 km on GLO-30. Mean slope fell from 17.01° to 16.72° and 16.70°: coarser cells average across steep valley sides. Slope-dependent products — erosion, wetness indices — are more resolution-sensitive than drainage areas.

5. Expect conditioning to matter more on coarse and noisy DEMs

The GLO-30 catchment was 490.46 km² after filling here, but 450.39 km² after least-cost breaching with a 1 km search in delineating a watershed. On the 10 m DEM the two methods agreed within 0.4 km². A noisier surface gives conditioning more to decide.

6. Do not assume the finest DEM is the most accurate

On the 1 m lidar DEM, the Birch Creek catchment at the gauge was 25.30 km² against 32.47 km² at 10 m and 33.14 km² at 30 m. The largest accumulation within 150 m of the gauge lay on the edge of the search window, 150 m away: the snap had stopped short of a larger channel, so a radius that worked at 10 m and 30 m failed at 1 m. A 1 m DEM also resolves ditches, road embankments and small channels that coarser DEMs smooth over, and any of them can re-route flow near an outlet — check the outlet on each DEM rather than reusing settings (Example 2).

7. Express thresholds in area

A 1,000-cell threshold is 0.001 km² at 1 m, 0.1 km² at 10 m and 0.9 km² at 30 m. In the Birch Creek basin it produced about 1,269 km of channel cells at 1 m, 46.4 km at 10 m and 19.1 km at 30 m. A 1 km² threshold gave 24.5, 18.6 and 18.4 km (Example 3).

8. Budget for runtime and file handling

Routing the 89.7-million-cell 1 m tile took 16.2 s; the 12.5-million-cell 10 m DEM 2.5 s; the 30 m grid 0.2–0.4 s. The 1 m file also needed rewriting first: it used a floating-point predictor, and WhiteboxTools stopped with The GeoTIFF reader does not currently support floating-point predictors (PREDICTOR=3) while its Python wrapper returned normally and wrote nothing.

Bar chart of the Birch Creek catchment area at the gauge from the 1 m lidar DEM, the 10 m 3DEP DEM and the 30 m GLO-30 DEM against the published area.
The same 150 m snap found the right channel at 10 m and 30 m, and stopped short of it at 1 m.

Code examples

Example 1 — resolution against source on the same basin 

import os
import time

import geopandas as gpd
import numpy as np
import rasterio
import whitebox
from exactextract import exact_extract
from pyproj import Transformer
from rasterio.features import rasterize
from rasterio.warp import Resampling, reproject
from scipy.ndimage import distance_transform_edt
from shapely.geometry import Point

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("."))

with rasterio.open("esopus_glo30_utm.tif") as glo:
    grid30 = glo.profile
with rasterio.open("esopus_3dep13_utm.tif") as src:
    averaged = np.full((grid30["height"], grid30["width"]), -9999.0, dtype="float32")
    reproject(rasterio.band(src, 1), averaged, src_nodata=src.nodata, dst_transform=grid30["transform"],
              dst_crs=grid30["crs"], dst_nodata=-9999.0, resampling=Resampling.average)
with rasterio.open("esopus_3dep_30m.tif", "w", **dict(grid30, dtype="float32", nodata=-9999.0)) as dst:
    dst.write(averaged, 1)


def route(dem, tag):
    start = time.perf_counter()
    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)
    seconds = time.perf_counter() - start
    with rasterio.open(f"{tag}_acc.tif") as src:
        return src.read(1), src.transform, src.crs, seconds


def watershed(tag, name, acc, transform, crs, lonlat, lonlat_crs, radius_m=150):
    x, y = Transformer.from_crs(lonlat_crs, crs, always_xy=True).transform(*lonlat)
    col, row = ~transform * (x, y)
    r, c, k = int(row), int(col), int(np.ceil(radius_m / transform.a))
    window = acc[r - k:r + k + 1, c - k:c + k + 1]
    dr, dc = np.unravel_index(window.argmax(), window.shape)
    gpd.GeoDataFrame(geometry=[Point(*(transform * (c - k + dc + 0.5, r - k + dr + 0.5)))], crs=crs).to_file(f"{name}_outlet.shp")
    wbt.watershed(f"{tag}_d8.tif", f"{name}_outlet.shp", f"{name}_ws.tif")
    with rasterio.open(f"{name}_ws.tif") as src:
        ws = src.read(1)
        return (ws != src.nodata) & (ws > 0)


esopus = gpd.read_file("nldi_basin_01362500.geojson").to_crs(32618)
flowlines = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(32618)
major = gpd.clip(flowlines[flowlines.ftype.isin([460, 558]) & (flowlines.streamorde >= 3)], esopus.geometry.make_valid())
samples = [part.interpolate(d) for g in major.geometry for part in getattr(g, "geoms", [g]) for d in np.arange(0, part.length, 10.0)]
px, py = np.array([p.x for p in samples]), np.array([p.y for p in samples])

routed = {}
for label, dem, tag in (("3DEP 10 m", "esopus_3dep13_utm.tif", "dep10"),
                        ("3DEP averaged to 30 m", "esopus_3dep_30m.tif", "dep30"),
                        ("GLO-30 30 m", "esopus_glo30_utm.tif", "glo30")):
    acc, transform, crs, seconds = route(dem, tag)
    routed[tag] = (label, acc, transform, crs)
    cell_km2 = transform.a ** 2 / 1e6
    ws = watershed(tag, tag, acc, transform, crs, (-74.2701944, 42.0144722), "EPSG:4269")
    inside = rasterize(esopus.geometry, out_shape=acc.shape, transform=transform).astype(bool)
    wbt.extract_streams(f"{tag}_acc.tif", f"{tag}_streams.tif", threshold=1.0 / cell_km2, zero_background=True)
    wbt.raster_streams_to_vector(f"{tag}_streams.tif", f"{tag}_d8.tif", f"{tag}_streams.shp")
    lines = gpd.clip(gpd.read_file(f"{tag}_streams.shp").set_crs(crs, allow_override=True), esopus)
    streams = (acc * cell_km2 >= 1.0) & inside
    distance = distance_transform_edt(~streams) * transform.a
    found = distance[((py - transform.f) / transform.e).astype(int), ((px - transform.c) / transform.a).astype(int)] <= 30
    wbt.slope(dem, f"{tag}_slope.tif", units="degrees")
    slope = exact_extract(f"{tag}_slope.tif", esopus, ["mean"], output="pandas")["mean"].iloc[0]
    print(f"{label:22} routed {seconds:4.1f} s; catchment {ws.sum() * cell_km2:6.2f} km2, IoU {(ws & inside).sum() / (ws | inside).sum():.4f}; "
          f"streams >= 1 km2 {lines.length.sum() / 1000:5.1f} km; major mapped streams found {found.mean():.1%}; mean slope {slope:.2f} deg")
3DEP 10 m              routed  2.5 s; catchment 492.45 km2, IoU 0.9853; streams >= 1 km2 336.7 km; major mapped streams found 94.0%; mean slope 17.01 deg
3DEP averaged to 30 m  routed  0.4 s; catchment 493.05 km2, IoU 0.9865; streams >= 1 km2 323.1 km; major mapped streams found 93.4%; mean slope 16.72 deg
GLO-30 30 m            routed  0.2 s; catchment 490.46 km2, IoU 0.9794; streams >= 1 km2 320.8 km; major mapped streams found 70.5%; mean slope 16.70 deg

Example 2 — a 1 m lidar DEM against 10 m and 30 m at a small gauge 

gauges = gpd.read_file("usgs_gauges_esopus.geojson")
birch = gauges.loc[gauges.monitoring_location_number == "013621955"].geometry.iloc[0]
birch_basin = gpd.read_file("nldi_basin_013621955.geojson")
with rasterio.open("birch_3dep1m_utm.tif") as src:
    print("lidar file structure:", src.tags(ns="IMAGE_STRUCTURE"))
    lidar_profile, lidar = src.profile, src.read(1)
with rasterio.open("birch_1m_for_whitebox.tif", "w", **dict(lidar_profile, compress="deflate", predictor=1)) as dst:
    dst.write(lidar, 1)                                            # WhiteboxTools cannot read PREDICTOR=3
acc, transform, crs, seconds = route("birch_1m_for_whitebox.tif", "lidar1")
routed["lidar1"] = ("3DEP 1 m lidar", acc, transform, crs)
print(f"1 m lidar DEM: {acc.size:,} cells, routed in {seconds:.1f} s")
for tag in ("lidar1", "dep10", "glo30"):
    label, acc, transform, crs = routed[tag]
    cell_km2 = transform.a ** 2 / 1e6
    ws = watershed(tag, f"{tag}_birch", acc, transform, crs, (birch.x, birch.y), "EPSG:4326")
    inside = rasterize(birch_basin.to_crs(crs).geometry, out_shape=acc.shape, transform=transform).astype(bool)
    area = ws.sum() * cell_km2
    print(f"{label:16} Birch Creek catchment {area:6.2f} km2 ({area / 32.37 - 1:+.1%} vs published 32.37), "
          f"IoU with NLDI {(ws & inside).sum() / (ws | inside).sum():.4f}, {int(ws.sum()):,} cells")
lidar file structure: {'COMPRESSION': 'DEFLATE', 'INTERLEAVE': 'BAND', 'PREDICTOR': '3'}
1 m lidar DEM: 89,654,775 cells, routed in 16.2 s
3DEP 1 m lidar   Birch Creek catchment  25.30 km2 (-21.8% vs published 32.37), IoU with NLDI 0.7508, 25,300,321 cells
3DEP 10 m        Birch Creek catchment  32.47 km2 (+0.3% vs published 32.37), IoU with NLDI 0.9517, 324,665 cells
GLO-30 30 m      Birch Creek catchment  33.14 km2 (+2.4% vs published 32.37), IoU with NLDI 0.9315, 36,825 cells

The search radius is the same 150 m on every DEM, which is 150 cells at 1 m and 5 cells at 30 m.

Example 3 — the same cell threshold at three resolutions 

for tag in ("lidar1", "dep10", "glo30"):
    label, acc, transform, crs = routed[tag]
    cell_km2 = transform.a ** 2 / 1e6
    inside = rasterize(birch_basin.to_crs(crs).geometry, out_shape=acc.shape, transform=transform).astype(bool)
    for rule, channel in (("1,000 cells", acc >= 1000), ("1 km2", acc * cell_km2 >= 1.0)):
        threshold_km2 = 1000 * cell_km2 if rule == "1,000 cells" else 1.0
        print(f"{label:16} {rule:12} ({threshold_km2:g} km2): about {(channel & inside).sum() * transform.a / 1000:6.1f} km of channel in the Birch Creek basin")
3DEP 1 m lidar   1,000 cells  (0.001 km2): about 1269.2 km of channel in the Birch Creek basin
3DEP 1 m lidar   1 km2        (1 km2): about   24.5 km of channel in the Birch Creek basin
3DEP 10 m        1,000 cells  (0.1 km2): about   46.4 km of channel in the Birch Creek basin
3DEP 10 m        1 km2        (1 km2): about   18.6 km of channel in the Birch Creek basin
GLO-30 30 m      1,000 cells  (0.9 km2): about   19.1 km of channel in the Birch Creek basin
GLO-30 30 m      1 km2        (1 km2): about   18.4 km of channel in the Birch Creek basin

Channel length here is channel cells times cell size, which undercounts diagonal steps; compare the rows with each other, not with mapped lengths.

Explanation

Why resolution alone changed so little

A drainage divide or a large channel is many cells wide at 10 m and still several cells wide at 30 m. Averaging preserves where the valleys and ridges are, so flow routes the same way and catchments barely move. What averaging loses is steepness — valley sides are flattened into their neighbours — so slope fell while area did not.

Why the source changed much more

GLO-30 is derived from radar and measures the top of whatever reflects: tree canopy over most of the Esopus basin. Canopy height varies from cell to cell independently of the ground, which can add false relief across valley floors. Routing on that surface can move a channel by a cell, and at 30 m a cell is the whole width of a small valley — consistent with the drop from 93.4% to 70.5% on the same grid.

Why finer is not automatically better

A 1 m bare-earth model shows the ground as it is, including features that coarser DEMs cannot: road embankments without their culverts, roadside ditches, field drains. Water in the model follows those features; water in reality goes through the culvert. The finer the DEM, the more such local routing decisions it contains, and the more conditioning and checking it needs.

Why cell thresholds do not transfer

Accumulation counts cells, and a cell's area grows with the square of its size. A threshold that draws a plausible network at one resolution draws a far denser or sparser one at another. Area thresholds carry the same meaning at every resolution.

Bar chart of approximate channel length in the Birch Creek basin from a 1,000-cell threshold and a 1 km2 threshold on the 1 m, 10 m and 30 m DEMs.
The same 1,000-cell threshold drew 1,269 km of channel at 1 m and 19 km at 30 m; 1 km² drew similar networks on all three.

Edge cases or notes

  • Mixed-source mosaics change behaviour at tile seams; check networks across boundaries.
  • Resampling method matters: average for elevation; nearest neighbour preserves artefacts; bilinear smooths slightly.
  • Geographic DEMs need area corrections before comparing catchments; see catchment areas on a latitude–longitude DEM.
  • Lidar tiles are often small; the Birch Creek tile covered only 11% of the Esopus basin.
  • GeoTIFF encoding can block tools: rewrite floating-point predictor files before WhiteboxTools reads them.
  • Hydro-enforced DEMs cut culverts and burn streams already; do not burn twice.
  • Runtime grows with the square of resolution; clip to a buffered catchment for 1 m work.

FAQ

Does DEM resolution change watershed area?

Less than you might expect for large basins. The Esopus catchment was 492.45 km² at 10 m and 493.05 km² with the same DEM averaged to 30 m.

Is a 30 m DEM good enough for stream networks?

A 30 m bare-earth DEM placed channels almost as well as 10 m. The 30 m GLO-30 surface model found only 70.5% of major mapped streams against 93.4%, because of its source, not its cell size.

Is a 1 m lidar DEM always better for hydrology?

No. At Birch Creek the same 150 m snap stopped short of the main channel on the 1 m DEM and gave a catchment 21.8% too small, where the 10 m DEM was within 0.3%.

Why does my stream threshold give a different network on another DEM?

A threshold in cells means a different area at each resolution: 1,000 cells is 0.001 km² at 1 m and 0.9 km² at 30 m. Set thresholds in area.

Why does WhiteboxTools fail to read my lidar GeoTIFF?

It does not support the floating-point predictor (PREDICTOR=3) and stopped without raising in Python. Rewrite the file with predictor=1 first.

How does resolution affect slope?

Coarser cells average over steep ground: mean catchment slope fell from 17.01° at 10 m to 16.72° at 30 m.