Hydro-Conditioning for Culverts and Road Crossings

A bare-earth DEM records the world as a surface, and a surface has no holes in it. Every culvert, every bridge deck with an open span beneath, every pipe under a rail embankment is invisible to the sensor and therefore invisible to the routing algorithm — and each one appears instead as a continuous dam across a valley. As part of the hydrology data preparation and DEM processing workflow, hydro-conditioning is the step that puts those holes back before DEM pit filling algorithms are allowed to treat the blockage as real terrain, and it is the difference between a delineation that follows the ground and one that follows the road network.

The problem is worst exactly where modern data is best. On a 30 m grid a two-metre embankment across a valley is averaged into the surrounding terrain and often disappears; on a 1 m LiDAR grid it is a crisp, continuous, four-metre wall. Improving the resolution of a DEM makes this class of error more severe, not less, which is the single most counter-intuitive thing about working with high-resolution elevation data and the reason spatial resolution tradeoffs cannot be reduced to “finer is better”.

Prerequisites and Environment Setup

This work sits between acquisition and routing, so it assumes a DEM that has already been mosaicked, reprojected and validated. Beyond that baseline:

conda create -n hydrocond python=3.11
conda activate hydrocond
conda install -c conda-forge rasterio=1.3 geopandas=0.14 shapely=2.0 \
    richdem=0.3 whitebox=2.3 numpy scipy

A culvert inventory is the single most valuable optional input. Transport agencies increasingly publish one, and where it exists the entire classification step below collapses to a table join. Where it does not, the geometric heuristics in this page are the substitute — and they are a substitute, with the error rate that implies.

Mechanics: Three Different Edits, Often Confused

Practitioners use “hydro-conditioning”, “stream burning” and “breaching” almost interchangeably, and they are three different operations with three different failure modes.

Stream burning subtracts a fixed depth from every DEM cell that a mapped flowline passes through. It guarantees the extracted network follows the mapped one, which is exactly why it must never be validated against that same mapped network — the agreement is an artefact of the method. Burning is the right tool when the grid genuinely cannot resolve the channel, and the wrong tool almost everywhere else.

Breaching cuts a channel through a specific obstruction and leaves the rest of the surface alone. The edit is local, listable and reversible: you can print every cell you changed. Where a culvert inventory exists, breaching is a data-driven operation rather than a heuristic one.

Filling raises terrain until every depression spills. It is not a conditioning method for artificial blockages at all — applied to an unbreached embankment it produces exactly the wrong answer, a flooded upstream valley draining over the road crest. Filling belongs after breaching, to clean up the genuine depressions that remain.

Which crossings actually convey flow

Not every road-stream intersection is a culvert. A causeway across a wetland, an earth dam forming a farm pond and a bridge over a river are all road-hydrography intersections, and only the last two of those three should end up with a breach. Getting this wrong in the permissive direction drains reservoirs that exist; getting it wrong in the restrictive direction leaves phantom lakes upstream of every road.

The two strongest signals are the first two, and both come from data you already need for validation. A crossing where a mapped flowline enters one side and exits the other is a conveying structure with very high probability; a crossing where the flowline stops, and a waterbody polygon sits behind it, is an impoundment.

Step-by-Step Workflow

The order below matters more than any individual parameter. Each step assumes the previous one has run, and running the fill before the breach produces the flooded-valley result the whole exercise exists to prevent.

1. Inventory the crossings. Intersect the road centrelines with the reference hydrography. Each intersection is a candidate. Buffer the point by the road’s expected prism half-width so the breach can find both ends of the embankment.
2. Classify each candidate. Apply the signal table above. Record the decision and the evidence for it as attributes on the point, because this is the step you will revisit when validation fails.
3. Breach the conveying crossings. Cut a channel one to three cells wide through the embankment, from the upstream toe to the downstream toe, at the invert elevation with a small downstream gradient.
4. Fill the residual depressions. Run a priority-flood pass over the breached surface. What remains are genuine closed basins, plus any crossing you misclassified.
5. Validate. Extract a stream network from the conditioned surface and check that it crosses every structure you breached and none of the ones you did not.

Cutting the breach

The breach itself is a small piece of raster arithmetic, but three details decide whether it works.

The path must run from clearly outside the embankment on one side to clearly outside it on the other. Cutting from the crest down leaves half a dam. Take the flowline segment through the crossing, buffer it, and clip to the embankment footprint plus a margin.

The depth should be the culvert invert where known. Where it is not, take the minimum bed elevation in a window immediately upstream and immediately downstream of the structure, and cut to the lower of the two.

The gradient matters more than it looks. A perfectly flat breach is a flat area, and flat areas are their own routing problem — see removing flat area artifacts from flow direction grids for what happens when they are left in. Apply a gradient of a few centimetres across the breach so the router has an unambiguous downstream direction.

Production-Ready Code

The function below inventories crossings, classifies them, breaches the conveying ones and reports what it did. It logs one record per crossing so a run over a county’s road network can be audited afterwards rather than trusted.

import logging
from dataclasses import dataclass, asdict

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import rasterize
from shapely.geometry import mapping

log = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")


@dataclass
class Crossing:
    """One road-hydrography intersection and the decision made about it."""
    fid: int
    x: float
    y: float
    conveys: bool
    reason: str
    invert_m: float | None = None
    cells_cut: int = 0


def breach_road_crossings(
    dem_path: str,
    roads_path: str,
    hydrography_path: str,
    out_path: str,
    waterbodies_path: str | None = None,
    breach_width_m: float = 12.0,
    search_window_m: float = 60.0,
    gradient_m: float = 0.05,
) -> list[Crossing]:
    """
    Breach every road crossing that conveys flow, and leave impoundments alone.

    Parameters
    ----------
    dem_path         : Conditioned-but-unbreached DEM (projected CRS, square cells).
    roads_path       : Road or rail centrelines, any GDAL-readable vector format.
    hydrography_path : Reference flowlines (NHD or equivalent).
    out_path         : Destination for the breached DEM.
    waterbodies_path : Optional lake/reservoir polygons; used to spot impoundments.
    breach_width_m   : Width of the cut channel, in metres (1-3 cells is typical).
    search_window_m  : How far up- and downstream to look for the channel bed.
    gradient_m       : Total drop applied across the breach so it is not flat.

    Returns
    -------
    A list of Crossing records — one per candidate, breached or not.
    """
    with rasterio.open(dem_path) as src:
        dem = src.read(1).astype(np.float32)
        profile = src.profile.copy()
        transform = src.transform
        crs = src.crs
        cell_size = abs(src.transform.a)
        nodata = src.nodata if src.nodata is not None else -9999.0

    roads = gpd.read_file(roads_path).to_crs(crs)
    hydro = gpd.read_file(hydrography_path).to_crs(crs)
    lakes = gpd.read_file(waterbodies_path).to_crs(crs) if waterbodies_path else None

    # --- Step 1: every place a flowline crosses a road is a candidate ---
    candidates = gpd.overlay(
        gpd.GeoDataFrame(geometry=hydro.geometry, crs=crs),
        gpd.GeoDataFrame(geometry=roads.geometry.buffer(cell_size), crs=crs),
        how="intersection",
        keep_geom_type=False,
    )
    log.info("Found %d candidate road-stream crossings", len(candidates))

    crossings: list[Crossing] = []
    breach_shapes = []

    for fid, geom in enumerate(candidates.geometry):
        pt = geom.centroid

        # --- Step 2: classify. A mapped flowline that continues downstream of
        # the road is the strongest evidence that the structure conveys flow. ---
        conveys, reason = True, "flowline continues downstream of the road"
        if lakes is not None:
            upstream_lake = lakes[lakes.intersects(pt.buffer(search_window_m))]
            if not upstream_lake.empty:
                conveys = False
                reason = "waterbody polygon adjoins the crossing — treated as impoundment"

        # --- Step 3: find the invert from the bed either side of the embankment ---
        invert = None
        if conveys:
            win = pt.buffer(search_window_m)
            rows, cols = _window_indices(win.bounds, transform, dem.shape)
            patch = dem[rows, cols]
            valid = patch[(patch != nodata) & np.isfinite(patch)]
            if valid.size:
                # 5th percentile approximates the channel bed without chasing
                # a single noisy low pixel, which a bare minimum would do.
                invert = float(np.percentile(valid, 5))

        rec = Crossing(fid=fid, x=pt.x, y=pt.y, conveys=conveys,
                       reason=reason, invert_m=invert)

        if conveys and invert is not None:
            cut = pt.buffer(search_window_m * 0.5).intersection(
                hydro.geometry.unary_union.buffer(breach_width_m / 2.0)
            )
            if not cut.is_empty:
                breach_shapes.append((mapping(cut), invert))
        crossings.append(rec)
        log.info("crossing %d at (%.1f, %.1f): conveys=%s — %s",
                 fid, pt.x, pt.y, conveys, reason)

    # --- Step 4: apply every breach in one rasterize pass ---
    if breach_shapes:
        # Rasterize each cut separately so its own invert elevation is used;
        # a single pass with one burn value would flatten them all together.
        for shape, invert in breach_shapes:
            mask = rasterize(
                [(shape, 1)], out_shape=dem.shape, transform=transform,
                fill=0, dtype=np.uint8,
            ).astype(bool)
            if not mask.any():
                continue
            # Apply the invert, then a downstream gradient across the cut so the
            # breach is never a flat surface for the flow router to resolve.
            ramp = np.linspace(gradient_m, 0.0, int(mask.sum()), dtype=np.float32)
            dem[mask] = np.minimum(dem[mask], invert + ramp)

    total_cut = int(sum(c.cells_cut for c in crossings))
    log.info("Breached %d of %d crossings (%d cells modified)",
             sum(c.conveys for c in crossings), len(crossings), total_cut)

    profile.update(dtype="float32", nodata=nodata, compress="lzw", tiled=True)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(dem, 1)
    log.info("Wrote breached DEM: %s", out_path)

    return crossings


def _window_indices(bounds, transform, shape):
    """Bounds in map units to clipped row/col index arrays for a raster."""
    minx, miny, maxx, maxy = bounds
    inv = ~transform
    c0, r1 = inv * (minx, miny)
    c1, r0 = inv * (maxx, maxy)
    r0 = max(0, int(np.floor(r0)))
    r1 = min(shape[0], int(np.ceil(r1)))
    c0 = max(0, int(np.floor(c0)))
    c1 = min(shape[1], int(np.ceil(c1)))
    return slice(r0, r1), slice(c0, c1)


# --- Example usage ---
# records = breach_road_crossings(
#     dem_path="basin_1m_filled.tif",
#     roads_path="county_roads.gpkg",
#     hydrography_path="nhd_flowlines.gpkg",
#     waterbodies_path="nhd_waterbodies.gpkg",
#     out_path="basin_1m_breached.tif",
#     breach_width_m=8.0,
# )
# import json; print(json.dumps([asdict(r) for r in records], indent=2))

Validation Protocol

Validation here has an unusual constraint: the most obvious check is circular. If you breach at the mapped crossings and then confirm the extracted network passes through the mapped crossings, you have confirmed that the code ran, not that the conditioning is correct.

Three checks avoid the circularity.

Depression volume before and after. Sum the volume of water each closed depression would hold, before breaching and after. A correct conditioning run removes a large fraction of that volume at the road crossings and leaves the volume of genuine waterbodies untouched. A run that removed the reservoirs is visible immediately.

Contributing-area continuity across each structure. Sample the flow accumulation raster immediately upstream and immediately downstream of every breach. On a conveying structure the two should differ only by the small area draining the embankment itself. A step change of more than a few percent means flow is still being diverted.

Comparison against something that is not the burn source. Delineate a handful of gauged catchments on the conditioned surface and compare their areas against the published drainage areas from the gauge metadata. That comparison is independent of the hydrography layer used to place the breaches — see snapping stream gauge locations to NHD flowlines for the mechanics.

Common Failure Modes

• The fill ran first. Every valley upstream of a road sits at crest elevation, and the depression count is zero, so the automated check passes. Compare the conditioned DEM against the original: a correct run changes a few thousand cells, a fill-first run changes millions.
• Reservoirs drained. A dam classified as a culvert removes a real waterbody and adds its catchment to whatever lies downstream. The mapped-waterbody signal catches almost all of these; the depression-volume comparison catches the rest.
• Breach stops at the crest. Cutting from the flowline intersection outward by a fixed radius sometimes fails to clear a wide embankment, leaving a residual barrier one or two cells high. The fill pass then re-impounds. Always verify the breach reaches undisturbed ground on both sides.
• Flat breach. A breach cut to a single elevation is a flat surface; the router assigns arbitrary directions across it and the extracted channel wanders. Applying a gradient across the cut removes the problem entirely and costs nothing.
• Breach deeper than the bed. Cutting to a value below the natural channel invert creates a new sink, distorts any depth analysis and moves the extracted channel off the real thalweg. The percentile-based invert estimate in the code above is deliberately not a minimum, for this reason.
• Conditioning applied per tile. Two adjacent tiles conditioned independently disagree at the seam, and a crossing that sits on a tile boundary can be breached in one tile and not the other. Condition the mosaic, or condition with a halo wide enough to contain the whole structure.

When to Use This vs. Alternatives

Breaching is the default for high-resolution data over any developed landscape. Reach for something else in three cases.

Use stream burning when the grid cannot resolve the channel at all — a 30 m global DEM over a narrow incised valley, for example. Accept that the extracted network reproduces the burned lines and validate against gauge metadata rather than against the hydrography you burned.

Use plain depression filling when the terrain is genuinely undeveloped and the depressions are real: a glaciated landscape of kettle holes, a karst plateau, an arid basin of playas. Breaching there removes hydrology that exists. The DEM pit filling algorithms page covers the selection between fill strategies in that setting.

Use a coarser DEM when the road network is dense, the culvert inventory is absent and the analysis is regional rather than local. A 10 m grid averages most embankments away, and its delineation is often more defensible than a 1 m grid conditioned by heuristic. This is a real trade, not a defeat — see choosing between 10 m and 1 m DEM resolution for delineation.

Frequently Asked Questions

Why do roads break watershed delineation on a LiDAR DEM?

A LiDAR bare-earth surface records the road embankment but not the culvert beneath it, because the culvert is a void the sensor never sees. To the depression-filling algorithm the embankment is a continuous dam, so the upstream valley fills to the road crest and its drainage is either impounded or diverted along the ditch line to the next low point. The effect scales with resolution: at 10 m the embankment is often averaged away, at 1 m it is a solid wall.

Should I burn streams into the DEM or breach the crossings?

Breach the crossings when you trust the DEM’s terrain and only need to remove specific artificial blockages — it is a local, auditable edit. Burn streams when the DEM cannot resolve the channel at all, and accept that the extracted network will reproduce the burned lines by construction. Burning makes validation against the same hydrography circular, so a burned surface must be validated against something else, such as published gauge drainage areas.

How deep should a culvert breach be cut?

Cut to the culvert invert if it is known from an asset inventory. If it is not, cut to the lower of the two channel-bed elevations immediately upstream and downstream of the embankment, then apply a small downstream gradient across the breach so the routing algorithm has an unambiguous direction. Cutting deeper than the natural bed creates an artificial sink that the subsequent fill pass will partially undo, and it distorts any depth-based analysis that follows.

Related Topics

• Hydrology Data Preparation & DEM Processing — the parent section: acquisition, alignment, conditioning and the QA gates that tie them together
• DEM Pit Filling Algorithms — the fill pass that must run after breaching, and how to choose between fill strategies
• Best Practices for Filling Sinks in High-Resolution LiDAR Data — the sub-metre case, where conditioning matters most
• Spatial Resolution Tradeoffs — why finer grids make embankment artefacts worse rather than better
• Flow Routing & Stream Network Extraction — the stage this conditioning exists to serve