Watershed Modeling Python for hydrology automation

Tuning Flow Accumulation Thresholds for Ephemeral Streams

Ephemeral and intermittent channels activate only during storm events and are absent from standard software defaults tuned for perennial rivers. As part of stream threshold tuning within the broader flow routing and stream network extraction workflow, selecting the right accumulation threshold for these transient channels requires replacing a static cell-count guess with an iterative, validation-driven calibration loop that measures spatial agreement against a reference network.

Prerequisites

Beyond the base environment (a conditioned DEM, a working WhiteboxTools or richdem install, and rasterio/geopandas), this technique requires:

  • A reference stream network — NHD flowlines, mapped ephemeral channels from field surveys, or high-resolution aerial photo-interpreted lines — rasterized or vector, covering the calibration basin.
  • scipy and numpy for the pixel-comparison arithmetic.
  • A DEM that has already passed DEM pit filling so depression artifacts do not create false headwater initiation points.
  • Python logging configured at INFO level; the sweep loop writes one log record per threshold step so long runs can be monitored without blocking.

Why Static Defaults Fail for Intermittent Channels

Ephemeral streams lack continuous baseflow. Their channel initiation points are governed by short-duration, high-intensity storms rather than sustained groundwater discharge. Standard thresholds optimised for perennial rivers typically:

  • Over-delineate in humid regions — gentle topographic hollows and agricultural swales accumulate drainage area quickly, producing phantom networks where no geomorphically active channel exists.
  • Under-delineate in arid and semi-arid basins — missing critical headwater conveyance paths that only activate during monsoonal pulses or convective bursts.

The threshold selection must also account for the routing algorithm in use. D-Infinity routing patterns disperse flow across two downslope neighbours rather than committing it to one cell, so peak accumulation values are systematically lower than D8 on the same DEM — the calibrated threshold will be 20–40 % smaller when switching algorithms.

Core Technique: The Threshold Sweep Loop

The fundamental idea is straightforward: instead of picking one threshold and hoping, you extract a binary stream mask for every candidate value in a range and score each mask against a reference network using overlap metrics. The value with the best F1-score becomes the operating threshold.

Two spatial metrics drive the decision:

  • Jaccard index (intersection over union): TP / (TP + FP + FN). Harsh toward both false positives and false negatives; useful for comparing across basins of different sizes.
  • F1-score (harmonic mean of precision and recall): 2·P·R / (P + R). More forgiving of class imbalance (most pixels are non-stream) and therefore the primary selection criterion.

The diagram below shows how F1 and Jaccard typically behave across a threshold sweep on a semi-arid 10 m DEM. Both metrics peak near the true optimum, but F1 has a broader, more stable plateau that makes the choice more robust to reference network imperfections.

Threshold Sweep Metric Curves Line chart showing F1-score and Jaccard index plotted against flow accumulation threshold from 100 to 1200 cells. Both curves rise steeply then decline; F1 peaks near 500 cells at 0.77 and Jaccard peaks at the same point at 0.62. 100 300 500 700 1000 1200 Threshold (cells) 0.0 0.2 0.4 0.6 0.8 Score optimal: 500 cells F1-score Jaccard

The sweep is fast on a typical 10 m DEM: comparing binary rasters with NumPy at 20 threshold steps takes under two seconds. The bottleneck is the upstream work — computing the accumulation raster once with WhiteboxTools, which runs in seconds to minutes depending on basin size.

Threshold Calibration Workflow

The calibration pipeline has five deterministic steps.

Threshold Calibration Pipeline Five sequential steps connected by arrows: Preprocess DEM, Compute Flow Accumulation, Define Candidate Range, Sweep and Validate, Select Optimal Threshold. Preprocess DEM Flow Accumulation Define Candidate Range Sweep & Validate Select Optimal Threshold ① Condition ② D8 / D∞ ③ Resolution-scaled ④ Jaccard / F1 ⑤ + slope mask
  1. Preprocess the DEM. Apply DEM pit filling algorithms (breach or fill sinks) and compute D8 flow direction. Conditioning errors propagate into every downstream threshold test, so this step must be clean before sweeping.
  2. Generate the accumulation raster. Compute upslope contributing area in cell counts (out_type="cells") or square metres. Cell counts simplify cross-resolution comparison; area units make thresholds portable across DEM sources.
  3. Define the candidate range. Anchor bounds to DEM resolution and regional precipitation intensity (see the parameter table below). A 10 m DEM typically needs a range of 100–2 000 cells for semi-arid conditions; a 1 m LiDAR DEM needs 20–300 cells for the same basin.
  4. Sweep and validate. For each candidate value, extract a binary stream mask and compare it against the reference network using Jaccard index and F1-score. Log every result.
  5. Select the optimal value. Pick the threshold that maximises F1. Apply a minimum slope mask (≥ 2–3 %) before finalising to suppress false positives in flat depositional zones.

Annotated Code Example

The function below integrates all five steps. It uses whitebox for hydrologic processing, rasterio for I/O, and numpy/scipy for spatial validation, and writes a log record for each iteration so long sweeps remain observable.

python 
import logging
import os
import numpy as np
import rasterio
from rasterio.features import rasterize
from whitebox import WhiteboxTools
from shapely.geometry import mapping
import geopandas as gpd

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


def tune_ephemeral_threshold(
    dem_path: str,
    ref_shp_path: str,
    threshold_range: tuple = (100, 2000),
    step: int = 100,
    buffer_cells: int = 2,
    out_dir: str = "./tuning_output",
) -> dict:
    """
    Sweep a range of flow accumulation thresholds and return the value
    that maximises F1-score against a reference ephemeral stream network.

    Parameters
    ----------
    dem_path        : Path to the input DEM (GeoTIFF, projected CRS).
    ref_shp_path    : Path to reference stream vector (any GDAL-readable format).
    threshold_range : (min, max) candidate cell-count thresholds.
    step            : Increment between candidates.
    buffer_cells    : Spatial tolerance — buffer applied to reference lines
                      before rasterizing (cells, not metres).
    out_dir         : Directory for intermediate WhiteboxTools outputs.

    Returns
    -------
    dict with keys: threshold, jaccard, f1, precision, recall
    """
    os.makedirs(out_dir, exist_ok=True)

    # --- Step 1: condition DEM and compute D8 flow accumulation ---
    wbt = WhiteboxTools()
    wbt.set_working_dir(out_dir)
    wbt.verbose = False  # suppress per-cell console output during sweep

    filled_path = os.path.join(out_dir, "dem_filled.tif")
    acc_path = os.path.join(out_dir, "flow_acc.tif")

    log.info("Filling depressions: %s", dem_path)
    wbt.fill_depressions(dem_path, filled_path)

    log.info("Computing D8 flow accumulation")
    wbt.d8_flow_accumulation(filled_path, acc_path, out_type="cells")

    # --- Step 2: rasterize reference network with spatial tolerance buffer ---
    with rasterio.open(dem_path) as src:
        meta = src.meta.copy()
        transform = src.transform
        cell_size = src.res[0]                    # metres per cell (projected CRS assumed)
        shape = (src.height, src.width)

        ref_gdf = gpd.read_file(ref_shp_path).to_crs(src.crs)
        buffer_dist = buffer_cells * cell_size    # convert cells → metres
        ref_buffered = ref_gdf.geometry.buffer(buffer_dist)

        # rasterize buffered reference lines to a binary truth mask
        ref_shapes = [
            (mapping(geom), 1)
            for geom in ref_buffered
            if geom is not None and not geom.is_empty
        ]
        ref_raster = rasterize(
            ref_shapes,
            out_shape=shape,
            transform=transform,
            fill=0,
            dtype=np.uint8,
        )
    log.info(
        "Reference network rasterized — positive cells: %d (%.2f %%)",
        ref_raster.sum(),
        100 * ref_raster.mean(),
    )

    # --- Step 3: load flow accumulation array ---
    with rasterio.open(acc_path) as acc_src:
        flow_acc = acc_src.read(1).astype(np.float32)

    # --- Step 4: threshold sweep ---
    results = []
    candidates = range(threshold_range[0], threshold_range[1] + step, step)

    for thresh in candidates:
        pred = (flow_acc >= thresh).astype(np.uint8).ravel()
        truth = ref_raster.ravel()

        tp = int(np.sum((pred == 1) & (truth == 1)))
        fp = int(np.sum((pred == 1) & (truth == 0)))
        fn = int(np.sum((pred == 0) & (truth == 1)))
        union = tp + fp + fn

        jaccard   = tp / union if union > 0 else 0.0
        precision = tp / (tp + fp) if (tp + fp) > 0 else 0.0
        recall    = tp / (tp + fn) if (tp + fn) > 0 else 0.0
        f1 = (
            2 * precision * recall / (precision + recall)
            if (precision + recall) > 0
            else 0.0
        )

        log.info(
            "thresh=%d  Jaccard=%.3f  F1=%.3f  prec=%.3f  rec=%.3f",
            thresh, jaccard, f1, precision, recall,
        )
        results.append(
            {"threshold": thresh, "jaccard": jaccard, "f1": f1,
             "precision": precision, "recall": recall}
        )

    # --- Step 5: select optimal threshold (max F1) ---
    best = max(results, key=lambda r: r["f1"])
    log.info(
        "Optimal threshold: %d cells  (F1=%.3f, Jaccard=%.3f)",
        best["threshold"], best["f1"], best["jaccard"],
    )
    return best


# --- Example usage ---
# result = tune_ephemeral_threshold(
#     dem_path="basin_10m.tif",
#     ref_shp_path="nhd_ephemeral_flowlines.shp",
#     threshold_range=(200, 1500),
#     step=50,
#     buffer_cells=2,
# )
# print(f"Use threshold {result['threshold']} cells (F1={result['f1']:.3f})")

Parameter Reference Table

Parameter Typical range Effect on extracted network
threshold_range lower bound 20–200 cells (1 m LiDAR) / 100–500 cells (10 m DEM) Too low → over-delineation; phantom channels in flat uplands
threshold_range upper bound 300–500 cells (1 m) / 1 500–3 000 cells (10 m) Too high → missed headwater channels; under-delineation in arid basins
step 10–50 (fine pass) / 50–200 (coarse pass) Smaller step = finer optimum but longer runtime
buffer_cells 1–3 Larger buffer forgives georeferencing offsets; >3 cells inflates recall artificially
Minimum slope mask ≥ 2 % (depositional zones) / ≥ 5 % (arid piedmonts) Suppresses false positives in playas, floodplains, irrigated fields
Routing algorithm D8 (default) / D∞ (divergent hillslopes) D∞ reduces peak accumulation by ~20–40 %; lower optimal threshold expected

Worked Example: Output Interpretation

A well-calibrated sweep on a 10 m DEM over a 50 km² semi-arid basin typically yields a curve where F1 rises steeply from the lower bound, peaks near 400–700 cells, then declines slowly as recall drops faster than precision improves. Concretely:

Threshold Jaccard F1 Precision Recall
100 0.41 0.58 0.44 0.86
300 0.55 0.71 0.63 0.82
500 0.62 0.77 0.74 0.80
700 0.58 0.73 0.81 0.67
1 000 0.44 0.61 0.88 0.47

In this example the optimal threshold is 500 cells (F1 = 0.77). The precision of 0.74 indicates that 26 % of extracted pixels are false positives — acceptable for planning-level analysis but worth reducing with a slope mask if the output feeds a hydraulic model. A recall of 0.80 means 20 % of reference channel pixels are missed, mainly at headwater tips where the DEM resolution cannot resolve the incision.

Anomalies to watch for:

  • F1 plateau across a wide range — the reference network may be coarser than the DEM; collect finer reference data or apply a skeletonisation step.
  • Precision > 0.90 but recall < 0.40 — the reference network is incomplete (common with NHD in headwater areas); supplement with field-mapped lines before trusting the metric.
  • Jaccard and F1 diverge sharply — usually caused by a large spatial imbalance (many non-stream cells); the Jaccard denominator grows with false positives, making it harsher than F1 in these cases.

Gotchas and Edge Cases

  • DEM resolution scaling. A 10 m DEM requires thresholds 4–9× higher than a 1 m LiDAR DEM for equivalent channel initiation, not 10×, because flow accumulation scales with the square of cell size when expressed in cells but linearly when expressed in area. If you switch DEM sources, rerun the sweep — do not rescale the old threshold arithmetically.
  • CRS mismatch between DEM and reference network. geopandas silently reprojects in the code above, but if coordinate reference system alignment is not verified first, the buffer distance (buffer_cells * cell_size) will be computed in the wrong units and produce an incorrectly sized tolerance zone.
  • Flat-area artefacts inflating accumulation. Depression filling creates flat plateaus where D8 assigns arbitrary flow directions, producing accumulation spikes that generate spurious channel initiations. Apply a flat-area routing correction (WhiteboxTools fix_flats or richdem’s FlatSurface resolver) before sweeping. See removing flat area artifacts from flow direction grids for the mechanics.
  • Two-pass refinement. Run a coarse sweep (step = 200) first, identify the approximate optimum, then run a fine sweep (step = 25) over a ±25 % window around it. Two-pass calibration typically reduces total iterations by 60 % while achieving equivalent or better resolution on the final threshold.
  • Algorithm-dependent re-calibration. If you switch from D8 to D-infinity routing to handle divergent flow on steep hillslopes, the entire sweep must be rerun on the new accumulation grid. D∞ peak accumulation values differ non-linearly from D8 and the transfer function is terrain-specific.

FAQ

What flow accumulation threshold should I use for ephemeral streams?

There is no universal value. Thresholds depend on DEM resolution, precipitation regime, and terrain type. For a 10 m DEM in an arid basin you might start at 200–500 cells; for 1 m LiDAR the equivalent is often 50–150 cells. Always calibrate against field-verified or NHD-derived reference lines using Jaccard index or F1-score rather than relying on software defaults.

Why do static cell-count thresholds over-delineate in humid regions?

In humid regions, gentle slopes accumulate drainage area quickly in topographic hollows and agricultural swales. A low threshold marks every such hollow as a stream channel even though no geomorphically active channel exists. Raising the threshold, combined with a minimum slope mask, filters these false positives.

How does D-Infinity routing change threshold selection for ephemeral channels?

D-Infinity routing patterns disperse flow across two downslope neighbours rather than committing it to a single cell. This lowers peak accumulation values at any given point compared to D8, so the optimal threshold is typically 20–40 % lower when switching from D8 to D-Infinity routing on the same DEM. Always rerun the full sweep on the new accumulation grid rather than scaling the old threshold.