How DEM Resolution and Source Change a Drainage Network
On this page
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.24 → 492.68 km², mean slope 17.01° → 16.72°, and 93.5% → 93.1% of major mapped streams found within 30 m by derived channels.
- The 30 m Copernicus GLO-30 surface model on the same grid lost 44 km² of catchment and found only 76.4% of the streams — the source, not the cell size, cost 9% of the area and 17 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.32 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.
Separate resolution from source, then compare
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 10 m DEM and its 30 m average gave 492.24 and 492.68 km², both within 1% of the published 497.3 km², with intersection-over-union against the USGS polygon of 0.9855 and 0.9869. Coarsening alone left the area alone. GLO-30 gave 448.62 km², 9.8% short, IoU 0.8965: on the same 30 m grid, the different source moved divides far enough to lose 44 km².
3. Compare channel positions
Derived channels of at least 1 km² found 93.5% of NHDPlus HR's major streams (order 3 and up) within 30 m on the 10 m DEM and 93.1% on its 30 m average, but 76.4% 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 335.0 km at 10 m to 322.5 km at 30 m, and was 324.3 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. Do not expect conditioning to rescue a poor source
The GLO-30 catchment was 448.62 km² after filling here, and 451.37 km² after least-cost breaching with a 1 km search in delineating a watershed. Both were about 9% short: when the surface puts a divide in the wrong place, neither method moves it back. The DEM's extent matters more than the method on a surface like this — see the hub guide's note on clipping.
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.32 km² at 10 m and 33.12 km² at 30 m. The largest accumulation within 150 m of the gauge lay on the edge of the search window, 149.6 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 18.8 km at 30 m. A 1 km² threshold gave 24.7, 18.6 and 18.3 km (Example 3).
8. Budget for runtime and file handling
Routing the 84.0-million-cell 1 m tile took 11.3 s; the 12.2-million-cell 10 m DEM 2.1 s; the 30 m grids 0.2–0.3 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.
Lidar against 10 m and 30 m at the same gauge
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 basin of Birch Creek (USGS 013621955), also from the 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); NHDPlus HR flowlines from the USGS National Map; USGS stream gauges from the Water Data API; the newest 1 m 3DEP lidar DEM around Birch Creek.
import os
import time
from pathlib import Path
import geopandas as gpd
import numpy as np
import pandas as pd
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', '013621955'):
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("nhdplushr_flowlines_esopus.gpkg"): # NHDPlus HR flowlines, paged
feats, off = [], 0
while True:
page = get("https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer/3/query", params={
"geometry": f"{W},{S},{E},{N}", "geometryType": "esriGeometryEnvelope", "inSR": 4326,
"spatialRel": "esriSpatialRelIntersects", "outFields": "*", "outSR": 4326, "f": "geojson",
"resultOffset": off, "resultRecordCount": 2000}).json()
feats += page["features"]
if len(page["features"]) < 2000:
break
off += len(page["features"])
flowlines = gpd.GeoDataFrame.from_features(feats, crs=4326)
flowlines.columns = [c.lower() if c != "geometry" else c for c in flowlines.columns]
flowlines.to_file("nhdplushr_flowlines_esopus.gpkg")
if not os.path.exists("usgs_gauges_esopus.geojson"): # USGS stream gauges in the basin box
feats, url = [], (f"https://api.waterdata.usgs.gov/ogcapi/v0/collections/monitoring-locations/items?f=json"
f"&bbox={W},{S},{E},{N}&site_type_code=ST&limit=500")
while url:
page = get(url).json()
feats += page["features"]
url = next((link["href"] for link in page.get("links", []) if link.get("rel") == "next"), None)
gauges = gpd.GeoDataFrame.from_features(feats, crs=4326)
gauges["drainage_area"] = pd.to_numeric(gauges.drainage_area, errors="coerce") # square miles
gauges.to_file("usgs_gauges_esopus.geojson", driver="GeoJSON")
if not os.path.exists("birch_3dep1m_utm.tif"): # newest 1 m 3DEP lidar DEM around Birch Creek
bb = gpd.read_file("nldi_basin_013621955.geojson").total_bounds + np.array([-0.006, -0.006, 0.006, 0.006])
items = get("https://tnmaccess.nationalmap.gov/api/v1/products", params={
"datasets": "Digital Elevation Model (DEM) 1 meter", "bbox": ",".join(map(str, bb)), "outputFormat": "JSON",
"max": 100}).json()["items"]
newest = max({i["title"].split()[-1] for i in items if i["title"].startswith("USGS 1 Meter")})
srcs = [rasterio.open("/vsicurl/" + i["downloadURL"]) for i in items if i["title"].endswith(newest)]
lidar, lidar_tf = merge(srcs, bounds=transform_bounds("EPSG:4326", srcs[0].crs, *bb))
with rasterio.open("birch_3dep1m_utm.tif", "w", driver="GTiff", height=lidar.shape[1], width=lidar.shape[2], count=1,
dtype="float32", crs=srcs[0].crs, transform=lidar_tf, nodata=srcs[0].nodata, compress="deflate",
predictor=3, tiled=True) as dst:
dst.write(lidar[0].astype("float32"), 1)
for f in ['nldi_basin_01362500.geojson', 'nldi_basin_013621955.geojson', 'esopus_3dep13_utm.tif', 'esopus_glo30_utm.tif', 'nhdplushr_flowlines_esopus.gpkg', 'usgs_gauges_esopus.geojson', 'birch_3dep1m_utm.tif']:
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
nldi_basin_013621955.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
nhdplushr_flowlines_esopus.gpkg: 4,739 features
usgs_gauges_esopus.geojson: 541 features
birch_3dep1m_utm.tif: 8,458 x 9,931 cells of 1 m
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_max_procs(1) # one thread: parallel runs give slightly different results
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.1 s; catchment 492.24 km2, IoU 0.9855; streams >= 1 km2 335.0 km; major mapped streams found 93.5%; mean slope 17.01 deg
3DEP averaged to 30 m routed 0.2 s; catchment 492.68 km2, IoU 0.9869; streams >= 1 km2 322.5 km; major mapped streams found 93.1%; mean slope 16.72 deg
GLO-30 30 m routed 0.3 s; catchment 448.62 km2, IoU 0.8965; streams >= 1 km2 324.3 km; major mapped streams found 76.4%; 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: 83,996,398 cells, routed in 11.3 s
3DEP 1 m lidar Birch Creek catchment 25.30 km2 (-21.8% vs published 32.37), IoU with NLDI 0.7507, 25,300,321 cells
3DEP 10 m Birch Creek catchment 32.32 km2 (-0.1% vs published 32.37), IoU with NLDI 0.9531, 323,244 cells
GLO-30 30 m Birch Creek catchment 33.12 km2 (+2.3% vs published 32.37), IoU with NLDI 0.9346, 36,801 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.1 km of channel in the Birch Creek basin
3DEP 1 m lidar 1 km2 (1 km2): about 24.7 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 18.8 km of channel in the Birch Creek basin
GLO-30 30 m 1 km2 (1 km2): about 18.3 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.
Why the source matters more than the resolution
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.1% to 76.4% 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.
Mosaics, resampling and hydro-enforced DEMs
- 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 about a tenth 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.
Internal links
- How to delineate a watershed from a pour point in Python — the workflow run on each DEM
- How to burn a known river network into a DEM — repairing channel alignment on coarse DEMs
- Fixing streams that stop at roads, bridges and embankments — what fine DEMs resolve
- How to extract a stream network from a DEM — thresholds as areas
- How to summarise catchment attributes: area, slope and land cover — slope at 10 m and 30 m
- Digital elevation models explained: DEM, DSM and DTM — surface versus bare-earth models
- How to reproject a raster with rasterio — averaging to a coarser grid
- Height above nearest drainage explained: a quick flood proxy — HAND on 10 m and 30 m DEMs
FAQ
Does DEM resolution change watershed area?
Less than you might expect for large basins. The Esopus catchment was 492.24 km² at 10 m and 492.68 km² with the same DEM averaged to 30 m. A different 30 m source, GLO-30, lost 9.8%.
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 76.4% of major mapped streams against 93.1%, 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.2%.
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.