🔧 Fix Hydrology & watersheds Terrain & elevation Hydrology & Watersheds · step 19 of 20 Updated

Fixing Streams That Stop at Roads, Bridges and Embankments

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 32 crossings (8%) and by at least 2 m at 11. The deepest, 6.27 m, was at State Route 28, and the five deepest were all on Routes 28 and 214.
  • State and secondary roads blocked far more often than local roads: 18 of 70 crossings against 13 of 305.
  • Least-cost breaching within 1 km removed every fill of 0.5 m or more at a crossing, cutting through the embankment at 31 crossings, by up to 7.27 m.
  • WhiteboxTools' burn_streams_at_roads cut crossings filled by 2 m or more from 11 to 6 but left all 32 above 0.5 m. On the 30 m GLO-30 surface model 59% of crossings sat in fills, the tool made it 63%, and it set one cell to 1,278.1 m where the DEM said 225.5 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.

Bar chart of stream-road crossings on the 3DEP 10 m DEM with fills of at least 0.5 m and at least 2 m after plain filling, after WhiteboxTools burn_streams_at_roads and after least-cost breaching.
Breaching removed every pond at a crossing; the road-burning tool mainly made the deep ones shallower.

Step-by-step solution

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, 32 crossings had at least 0.5 m and 11 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.27, 4.50, 4.06 and 2.55 m — and State Route 214 the other, 2.94 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 31 crossings, the deepest by 7.27 m, and left no crossing sitting in a fill of 0.5 m or more. It took 18.4 s on the 3DEP basin DEM against 0.8 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 11 to 6 and the deepest from 6.27 to 5.51 m, but all 32 crossings with half a metre of fill remained, where breaching had cleared them all.

7. Inspect the output of any burning tool

On GLO-30, burn_streams_at_roads with a 60 m width changed 834 cells. One of them, at a private road, came out at 1,278.1 m where the DEM held 225.5 m, and filling the burned DEM then raised the surrounding cells to 242.7 m — a 1,052.63 m "fill" in the summary. A single corrupted cell was enough to create a new pond. Compare any modified DEM with its source before routing it.

8. Use a bare-earth DEM

On the 30 m GLO-30 surface model, 246 of 419 crossings (59%) sat in fills of at least 0.5 m and 177 in fills of 2 m or more, with the deepest at 48.86 m. Breaching halved that to 128 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.

Bar chart of the share of stream crossings with at least half a metre of fill on the 3DEP DEM for state and secondary roads, local roads, private roads and trails.
A quarter of state-road crossings ponded the stream behind them; one in twenty-five local-road crossings did.

Code examples

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_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                               0.8 s; crossings filled >= 0.5 m:  32 (8%), >= 2 m:  11; deepest    6.27 m at State Rte 28 (S1200)
  burn_streams_at_roads 30 m + fill   1.2 s; crossings filled >= 0.5 m:  32 (8%), >= 2 m:   6; deepest    5.51 m at State Rte 28 (S1200)
  breach 1 km + fill                18.4 s; crossings filled >= 0.5 m:   0 (0%), >= 2 m:   0; deepest    0.50 m at Birch Creek Rd (S1400)
                                   breaching cut >= 0.5 m within 30 m of 31 crossings; deepest cut 7.27 m
GLO-30
  fill                               0.1 s; crossings filled >= 0.5 m: 246 (59%), >= 2 m: 177; deepest   48.86 m at Upper Boiceville Rd (S1400)
  burn_streams_at_roads 60 m + fill   0.2 s; crossings filled >= 0.5 m: 265 (63%), >= 2 m: 192; deepest 1052.63 m at (Private Rd) (S1400)
  breach 1 km + fill                 0.3 s; crossings filled >= 0.5 m: 128 (31%), >= 2 m:  68; deepest   48.91 m at Upper Boiceville Rd (S1400)
                                   breaching cut >= 0.5 m within 45 m of 110 crossings; deepest cut 15.72 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.50 m" is just below 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.06
S1400             305                      13      0.00
S1500              16                       1      0.00
S1710               6                       0      0.00
S1740              21                       0      0.00
S1780               1                       0      0.09
         road road_class  order  filled_m
 State Rte 28      S1200      1      6.27
 State Rte 28      S1200      4      4.50
 State Rte 28      S1200      3      4.06
State Rte 214      S1200      3      2.94
 State Rte 28      S1200      2      2.55

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

Explanation

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 246 times on GLO-30 and 32 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.

Two panels comparing stream-road crossings on the 3DEP 10 m bare-earth DEM and the GLO-30 surface model: share of crossings in fills, deepest fill and the effect of breaching and road burning.
On bare earth, road crossings are a local problem with a local fix; on a surface model they are part of a wider one.

Edge cases or notes

  • 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, 32 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 every crossing fill 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 11 to 6 but left all 32 half-metre fills; on GLO-30 it created a corrupted cell 1,052 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, 59% of crossings sat in fills and breaching left 31%; a bare-earth DEM removes most of the problem at source.