🔧 Fix Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 19 Notebook By · Updated

Fixing Streams That Stop at Roads, Bridges and Embankments

On this page

Problem statement

A stream network derived from a detailed DEM runs happily down a valley and then stops at a road, or turns and follows the road for a kilometre before finding a way through. The DEM is right: it records the embankment, and the culvert or bridge underneath is invisible to a surface model or to lidar ground returns. Depression filling then turns the land upstream of every embankment into a flat pond, and flow crosses it wherever the fill's arithmetic sends it.

Measured at the 419 places where USGS NHDPlus HR streams cross TIGER/Line 2024 roads in the Esopus Creek basin, New York:

  • On the 10 m 3DEP bare-earth DEM, filling raised the surface by at least 0.5 m within 30 m of 31 crossings (7%) and by at least 2 m at 8. The deepest, 6.78 m, was at State Route 28, which had four of the five deepest.
  • State and secondary roads blocked far more often than local roads: 18 of 70 crossings against 12 of 305.
  • Least-cost breaching within 1 km removed all but two fills of 0.5 m or more at a crossing — the two left were 0.63 m and less — cutting through the embankment at 35 crossings, by up to 7.19 m.
  • WhiteboxTools' burn_streams_at_roads cut crossings filled by 2 m or more from 8 to 5 but left 29 of the 31 above 0.5 m. On the 30 m GLO-30 surface model 60% of crossings sat in fills, the tool made it 63%, and it set one cell to 1,278.2 m where the DEM said 255.6 m.

Quick answer

wbt.breach_depressions_least_cost("dem.tif", "dem_breached.tif", dist=100, fill=True)   # 100 cells of 10 m

Breach rather than fill on a detailed DEM, so embankments are cut at their narrowest point instead of ponding the valley behind them. Then check the crossings: find where mapped streams meet roads and measure how much conditioning changed the DEM there.

Crossings still blocked after filling, road-burning and breaching, at two depth thresholds. 419 crossings, 3DEP 10 mstill blocked ≥ 0.5 mstill blocked ≥ 2 mfill!31!8road burn!29!5breach✓2✓0Filling raises the upstream side; breaching cuts the blockage itself spatialworkflow.io
Only one of the three removes the obstruction rather than drowning it.

Confirm, find every crossing, then breach

1. Confirm the stop is at a road

Overlay derived streams, mapped streams and roads. Streams that end or deflect exactly at a road line, with a flat area immediately upstream in the filled DEM, are blocked by the embankment. Streams that stop in open country are more often a threshold or conditioning problem; see fixing flow accumulation that stops early.

2. Find every crossing

Intersect mapped streams with roads (Example 1). In the Esopus basin that gave 419 crossings after merging points within 20 m: 305 on local roads (TIGER class S1400), 70 on state and secondary roads (S1200), and the rest on private drives, trails and ramps. Most were on small streams: 229 on first-order and 130 on second-order streams.

3. Measure the fill at each crossing

Fill the DEM without flat fixing and subtract the original. The largest rise within a cell or two of each crossing is the depth of the pond the embankment creates. On 3DEP 10 m, 31 crossings had at least 0.5 m and 8 at least 2 m (Example 2).

4. See which roads block

The worst crossings were not random. State Route 28 had four of the five deepest fills — 6.78, 4.51, 4.15 and 2.75 m — and a local road, Evergreen Drive, the other, 2.52 m. 26% of state and secondary road crossings had a fill of half a metre or more, against 4% of local road crossings (Example 3). Large embankments over small streams — highways across side valleys — are where to look first.

5. Breach instead of fill

Least-cost breaching searches for the cheapest cut from each depression to lower ground. Through a narrow embankment that is almost always straight through the road. With a 1 km search it cut at least 0.5 m within 30 m of 35 crossings, the deepest by 7.19 m, and left two crossings sitting in a fill of 0.5 m or more, the deeper at 0.63 m. It took 20.7 s on the 3DEP basin DEM against 1.1 s for filling.

6. Or burn streams at roads explicitly

WhiteboxTools' burn_streams_at_roads lowers the DEM where a stream vector crosses a road vector, within a given width. With a 30 m width, crossings with at least 2 m of fill fell from 8 to 5, but 29 of the 31 crossings with half a metre of fill remained and the deepest stayed at 6.78 m, where breaching had cleared all but two.

7. Inspect the output of any burning tool

On GLO-30, burn_streams_at_roads with a 60 m width changed 802 cells, and raised 205 of them instead of lowering them. The worst came out at 1,278.2 m where the DEM held 255.6 m — the 1,022.62 m "fill" at Glenford-Wittenberg Road in the summary is that one cell, 1 km above the terrain. A tool meant to lower crossings had built a tower in one. Compare any modified DEM with its source before routing it.

8. Use a bare-earth DEM

On the 30 m GLO-30 surface model, 250 of 419 crossings (60%) sat in fills of at least 0.5 m and 174 in fills of 2 m or more, with the deepest at 50.78 m. Breaching more than halved that, to 113, but could not remove it. Road-specific fixes cannot repair a surface that is wrong across whole valley floors; see how DEM resolution and source change a drainage network.

Share of stream crossings with at least half a metre of fill on the 3DEP DEM, by road class. crossings ponded by at least 0.5 mstate and secondary18 of 70local roads12 of 305vehicular trails1 of 16private roads0 of 21030%TIGER/Line 2024 road classes S1200, S1400, S1500 and S1740 spatialworkflow.io
A quarter of state-road crossings ponded the stream behind them; one in twenty-five local-road crossings did.

Fill, burn or breach at road crossings

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; USGS 3DEP 1/3 arc-second elevation, resampled bilinearly to UTM 18N at 10 m over the basin plus 1 km (esopus_3dep13_utm_basin.tif); Copernicus GLO-30 resampled the same way at 30 m (esopus_glo30_utm_basin.tif); NHDPlus HR flowlines from the USGS National Map; Census TIGER 2024 roads for Ulster, Greene and Delaware counties.

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", 1000, "esopus_3dep13_utm_basin.tif")   # 3DEP 1/3 arc-second at 10 m, basin + 1 km
dem_box("glo30", 1000, "esopus_glo30_utm_basin.tif")    # Copernicus GLO-30 at 30 m, basin + 1 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")
for county in ("36111", "36039", "36025"):              # TIGER 2024 roads: Ulster, Greene, Delaware
    if not os.path.exists(f"tl_2024_{county}_roads.zip"):
        Path(f"tl_2024_{county}_roads.zip").write_bytes(get(f"https://www2.census.gov/geo/tiger/TIGER2024/ROADS/tl_2024_{county}_roads.zip").content)

for f in ['nldi_basin_01362500.geojson', 'esopus_3dep13_utm_basin.tif', 'esopus_glo30_utm_basin.tif', 'nhdplushr_flowlines_esopus.gpkg', 'tl_2024_36111_roads.zip', 'tl_2024_36039_roads.zip', 'tl_2024_36025_roads.zip']:
    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_basin.tif: 3,471 x 2,408 cells of 10 m
esopus_glo30_utm_basin.tif: 1,157 x 803 cells of 30 m
nhdplushr_flowlines_esopus.gpkg: 4,739 features
tl_2024_36111_roads.zip: 2.6 MB
tl_2024_36039_roads.zip: 1.2 MB
tl_2024_36025_roads.zip: 2.4 MB

Example 1 — find stream–road crossings

import geopandas as gpd
import pandas as pd

basin = gpd.read_file("nldi_basin_01362500.geojson").to_crs(32618)
outline = basin.geometry.make_valid().iloc[0]
roads = pd.concat([gpd.read_file(f"zip://tl_2024_{county}_roads.zip") for county in ("36111", "36039", "36025")])
roads = roads.to_crs(32618).reset_index(drop=True)
roads = roads[roads.intersects(outline.buffer(2000))]
flowlines = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(32618)
streams = flowlines[flowlines.ftype == 460]                    # stream/river only

pairs = gpd.sjoin(streams[["streamorde", "geometry"]], roads[["MTFCC", "FULLNAME", "geometry"]], predicate="intersects")
found = []
for stream_line, order, road_id in zip(pairs.geometry, pairs.streamorde, pairs.index_right):
    hit = stream_line.intersection(roads.geometry.loc[road_id])
    for point in getattr(hit, "geoms", [hit]):
        if point.geom_type == "Point":
            found.append({"order": order, "road_class": roads.MTFCC.loc[road_id], "road": roads.FULLNAME.loc[road_id], "geometry": point})
crossings = gpd.GeoDataFrame(found, crs=32618)
crossings = crossings[crossings.within(outline)]
crossings = crossings[~crossings.geometry.map(lambda p: (round(p.x / 20), round(p.y / 20))).duplicated()].reset_index(drop=True)
print(f"{len(crossings)} stream-road crossings; road classes {crossings.road_class.value_counts().to_dict()}; "
      f"stream orders {crossings.order.value_counts().sort_index().to_dict()}")
419 stream-road crossings; road classes {'S1400': 305, 'S1200': 70, 'S1740': 21, 'S1500': 16, 'S1710': 6, 'S1780': 1}; stream orders {1: 229, 2: 130, 3: 57, 4: 3}

The three county road files are concatenated, so reset_index is needed before looking roads up by the index that sjoin returns.

Example 2 — fill, burn or breach, measured at the crossings

import os
import time

import numpy as np
import rasterio
import whitebox
from scipy.ndimage import maximum_filter

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("."))
streams[["geometry"]].to_file("streams.shp")
roads[["geometry"]].to_file("roads.shp")


def change_at_crossings(raw_path, conditioned_path, radius_m):
    with rasterio.open(raw_path) as src:
        raw, transform = src.read(1, masked=True), src.transform
    with rasterio.open(conditioned_path) as src:
        conditioned = src.read(1)
    change = np.where(raw.mask, 0.0, conditioned.astype("float64") - raw.data)
    k = max(1, round(radius_m / transform.a))
    dy, dx = np.mgrid[-k:k + 1, -k:k + 1]
    disc = dx ** 2 + dy ** 2 <= k * k
    cols, rows = ~transform * (crossings.geometry.x.values, crossings.geometry.y.values)
    rows, cols = rows.astype(int), cols.astype(int)
    return maximum_filter(change, footprint=disc)[rows, cols], maximum_filter(-change, footprint=disc)[rows, cols]


def summarise(label, raised, seconds):
    deepest = int(np.argmax(raised))
    print(f"  {label:32} {seconds:5.1f} s; crossings filled >= 0.5 m: {int((raised >= 0.5).sum()):3d} ({(raised >= 0.5).mean():.0%}), "
          f">= 2 m: {int((raised >= 2).sum()):3d}; deepest {raised.max():7.2f} m at {crossings.road[deepest]} ({crossings.road_class[deepest]})")


results = {}
for label, dem, radius, width in (("3DEP 10 m", "esopus_3dep13_utm_basin.tif", 30, 30), ("GLO-30", "esopus_glo30_utm_basin.tif", 45, 60)):
    print(label)
    tag = label.split()[0].lower().replace("-", "")
    start = time.perf_counter()
    wbt.fill_depressions(dem, f"{tag}_filled.tif", fix_flats=False)
    filled, _ = change_at_crossings(dem, f"{tag}_filled.tif", radius)
    summarise("fill", filled, time.perf_counter() - start)

    start = time.perf_counter()
    wbt.burn_streams_at_roads(dem, "streams.shp", "roads.shp", f"{tag}_burned.tif", width=width)
    wbt.fill_depressions(f"{tag}_burned.tif", f"{tag}_burned_filled.tif", fix_flats=False)
    after_burn, _ = change_at_crossings(dem, f"{tag}_burned_filled.tif", radius)
    summarise(f"burn_streams_at_roads {width} m + fill", after_burn, time.perf_counter() - start)

    with rasterio.open(dem) as src:
        cell = src.res[0]
    start = time.perf_counter()
    wbt.breach_depressions_least_cost(dem, f"{tag}_breached.tif", dist=int(1000 / cell), fill=True)
    after_breach, cut = change_at_crossings(dem, f"{tag}_breached.tif", radius)
    summarise("breach 1 km + fill", after_breach, time.perf_counter() - start)
    print(f"  {'':32} breaching cut >= 0.5 m within {radius} m of {int((cut >= 0.5).sum())} crossings; deepest cut {cut.max():.2f} m")
    results[label] = dict(filled=filled, after_burn=after_burn, after_breach=after_breach, cut=cut, tag=tag, dem=dem, radius=radius)
3DEP 10 m
  fill                               1.1 s; crossings filled >= 0.5 m:  31 (7%), >= 2 m:   8; deepest    6.78 m at State Rte 28 (S1200)
  burn_streams_at_roads 30 m + fill   1.7 s; crossings filled >= 0.5 m:  29 (7%), >= 2 m:   5; deepest    6.78 m at State Rte 28 (S1200)
  breach 1 km + fill                20.7 s; crossings filled >= 0.5 m:   2 (0%), >= 2 m:   0; deepest    0.63 m at Erikas Walk (S1400)
                                   breaching cut >= 0.5 m within 30 m of 35 crossings; deepest cut 7.19 m
GLO-30
  fill                               0.1 s; crossings filled >= 0.5 m: 250 (60%), >= 2 m: 174; deepest   50.78 m at Coldbrook Rd (S1400)
  burn_streams_at_roads 60 m + fill   0.3 s; crossings filled >= 0.5 m: 263 (63%), >= 2 m: 189; deepest 1022.62 m at Glenford-Wittenberg Rd (S1400)
  breach 1 km + fill                 0.3 s; crossings filled >= 0.5 m: 113 (27%), >= 2 m:  68; deepest   50.82 m at Coldbrook Rd (S1400)
                                   breaching cut >= 0.5 m within 45 m of 98 crossings; deepest cut 15.50 m

The radius covers a cell or two either side of the crossing: 30 m on the 10 m grid, 45 m on the 30 m grid. The breached line's deepest, 0.63 m, is one of the two crossings breaching left above the 0.5 m threshold.

Example 3 — which roads block streams

blocked = crossings.assign(filled_m=results["3DEP 10 m"]["filled"])
print(blocked.groupby("road_class").filled_m.agg(crossings="size", filled_over_half_metre=lambda s: int((s >= 0.5).sum()),
                                                  median_m="median").round(2).to_string())
print(blocked.nlargest(5, "filled_m")[["road", "road_class", "order", "filled_m"]].round(2).to_string(index=False))
            crossings  filled_over_half_metre  median_m
road_class                                             
S1200              70                      18      0.04
S1400             305                      12      0.00
S1500              16                       1      0.00
S1710               6                       0      0.00
S1740              21                       0      0.00
S1780               1                       0      0.17
        road road_class  order  filled_m
State Rte 28      S1200      1      6.78
State Rte 28      S1200      4      4.51
State Rte 28      S1200      3      4.15
State Rte 28      S1200      2      2.75
Evergreen Dr      S1400      1      2.52

TIGER's MTFCC codes: S1200 secondary roads, S1400 local roads, S1500 vehicular trails, S1710 walkways, S1740 private roads, S1780 parking lot roads.

Why a DEM dams streams at roads

Why a DEM dams streams at roads

Lidar and photogrammetric DEMs record the surface the sensor sees. An embankment is solid ground; the culvert pipe through it is not visible. Bare-earth processing removes bridges, which span open space, but not culverts, which are buried. The DEM therefore shows a wall across the valley with no gap, and every cell upstream is a depression.

Why filling makes it worse

Filling raises the depression to the height of the lowest point on its rim. Behind a road, that rim is often the road surface itself, several metres above the stream bed. The valley upstream becomes a flat plateau, flow directions across it are set by flat-resolution rules rather than terrain, and the derived stream wanders across the plateau before crossing the road wherever the rim is lowest.

Why breaching usually fixes it

The cheapest way out of a pond behind a road is straight through the embankment, where the barrier is narrowest. Least-cost breaching finds that path and lowers it, so the derived channel crosses where the culvert is. It can also cut in the wrong place when a longer, shallower route exists, which is why a cost or distance cap is needed.

Why the surface model is different

GLO-30 is a digital surface model, and 98% of the Esopus basin is tree cover, so its elevations include canopy wherever trees stand. The same 419 crossings of the same roads sat in half-metre fills 250 times on GLO-30 and 31 times on the bare-earth DEM, so most GLO-30 fills cannot be blamed on the embankments alone — which is why no road-specific fix cleared them.

Stream-road crossings on the 3DEP 10 m bare-earth DEM and the GLO-30 surface model. the same 419 crossings, two DEMs3DEP 10 m bare earth✓31 crossings in fills (7%)✓deepest 6.78 m, State Rte 28✓breaching left 2, at 0.63 m✓road burning left 29GLO-30 surface model✕250 crossings in fills (60%)✕deepest 50.78 m✕breaching left 113✕road burning: 263, one 1,278 m cellFills of at least half a metre within a cell or two of each crossing spatialworkflow.io
On bare earth, road crossings are a local problem with a local fix; on a surface model they are part of a wider one.

Bridges, culvert inventories and weirs

  • Bridges in lidar DTMs are usually removed; culverts are not.
  • Hydro-enforced DEMs from national agencies may already cut known culverts; check the metadata before breaching again.
  • Culvert inventories give exact locations; burning at surveyed culverts is better than at every mapped crossing.
  • Dams and weirs are real barriers; exclude them from breaching and burning.
  • Roads along valley floors can capture streams for long distances; a larger breach distance may be needed.
  • Derived HAND and flood proxies are too high upstream of unbreached embankments; see calculating HAND in Python.
  • Duplicate crossing points occur where a stream touches a road more than once; merge points within a cell or two.

FAQ

Why do my DEM streams stop at roads?

The DEM records the road embankment but not the culvert beneath it, so the stream is dammed. On the 3DEP 10 m DEM, 31 of 419 stream–road crossings in the Esopus basin sat in fills of at least half a metre.

How do I fix culverts in a DEM for flow routing?

Breach depressions with a least-cost method, or lower the DEM where streams cross roads. Least-cost breaching within 1 km cleared all but two of the 31 crossing fills of half a metre or more.

Does WhiteboxTools BurnStreamsAtRoads work?

Partly. On the 10 m DEM it reduced crossings with 2 m or more of fill from 8 to 5 but left 29 of the 31 half-metre fills; on GLO-30 it created a corrupted cell 1,023 m above the DEM. Check its output.

Which roads block streams most in a DEM?

Large embankments over small streams. A quarter of state and secondary road crossings had a half-metre fill, against one in twenty-five local road crossings.

Should I fill or breach a DEM with roads?

Breach. Filling ponds the valley behind each embankment up to the road surface, while breaching cuts through the embankment where the culvert is.

Can I fix road crossings on a 30 m global DEM?

Only partly. On GLO-30, 60% of crossings sat in fills and breaching left 27%; a bare-earth DEM removes most of the problem at source.