Pipeline Orchestration: DEM to Watershed Automation

A single watershed delineation run touches a dozen distinct operations — download a tile, reproject it, fill depressions, compute flow direction and accumulation, threshold a stream network, snap outlets, delineate catchments, and validate the result. Run that chain by hand once and it is a script; run it across a hundred basins on a schedule and it becomes an orchestration problem. This guide, part of the Watershed Delineation & Catchment Synchronization section, treats the full chain as a managed directed acyclic graph (DAG) executed under Prefect or Airflow, so that reruns are cheap, failures are isolated, and every intermediate artifact is traceable.

The difference between a shell script and an orchestrated pipeline is not the algorithms — those are unchanged — but the guarantees around them: idempotent tasks that can rerun without corrupting state, content-addressed caching that skips work already done, retries with backoff on flaky network stages, parametrization over tiles or hydrologic unit codes (HUCs), structured observability, and failure isolation so one bad basin does not sink a batch of ninety-nine good ones. Those guarantees are what let you point the same flow definition at a nightly HUC-8 batch and trust the output.

Prerequisites & Environment Setup

Orchestration adds a scheduling layer on top of the geospatial stack; both must be installed cleanly. Keep the orchestrator and the compiled GDAL/RichDEM binaries in the same environment so tasks import the same library versions the scheduler runs under.

conda create -n hydro-orchestration python=3.10 richdem rasterio geopandas pysheds -c conda-forge
conda activate hydro-orchestration
pip install "prefect>=2.14" "apache-airflow>=2.8"

Install prefect and apache-airflow with pip after the conda-forge geospatial base so the orchestrators do not pull in a conflicting GDAL. You rarely deploy both in production — pick one per environment — but installing both locally lets you prototype the same task graph two ways before committing. The remainder of this guide develops a Prefect flow in full and points to the two focused companion guides for framework-specific detail: building a Prefect flow for DEM-to-watershed automation and Airflow DAG templates for batch watershed delineation.

The Pipeline as a Managed DAG

Each processing stage becomes a node; each data dependency becomes an edge. Because water flows downhill and never uphill in the graph sense — accumulation cannot precede flow direction, delineation cannot precede a snapped outlet — the dependency structure is naturally acyclic. The orchestrator’s job is to walk that graph in topological order, run each node once its inputs are ready, cache the ones worth caching, retry the ones that fail transiently, and record what happened.

The stage boundaries in the diagram are deliberate. Each corresponds to a well-understood operation covered elsewhere on this site — depression filling draws on DEM pit filling algorithms, the routing stage on D8 flow direction implementation, the channel mask on stream threshold tuning, and outlet handling on outlet point mapping and validation. Orchestration does not reimplement any of them; it wires them together with guarantees.

Mechanics: Task Graph, Dependencies, and State

Three concepts carry over from general workflow engineering into hydrology and deserve precise definitions before any code.

Tasks and dependencies. A task is the smallest unit the orchestrator schedules, retries, and caches independently. The natural granularity is one processing stage — one raster transform, one vector operation. Dependencies are declared implicitly by passing one task’s return value into the next, which is how both Prefect and Airflow’s TaskFlow API infer the DAG. Resist the temptation to fuse “fill plus route plus accumulate” into a single mega-task: coarse tasks defeat caching and make a mid-pipeline failure restart from the top.

Artifacts and state. Every task produces an artifact — usually a GeoTIFF or a GeoPackage on disk, sometimes an in-memory array handed to the next task. The orchestrator tracks run state (pending, running, completed, failed, cached) separately from the artifact itself. Keeping those two ideas distinct is what makes restarts safe: a task in the completed state with a present, valid artifact can be skipped entirely on rerun.

Idempotency. A task is idempotent when running it twice with the same inputs leaves the system in the same state as running it once. For raster stages this means writing to a deterministic output path, writing to a temporary file and atomically renaming on success, and never mutating an input in place. Non-idempotent writes are the single most common cause of corrupted orchestrated pipelines, discussed under Failure Modes below.

Result Caching

Deterministic stages — reproject, fill, flow direction, accumulation — always produce the same output for the same input, so recomputing them on rerun is pure waste. A cache key derived from the task’s inputs (the input file’s content hash plus the parameters) lets the orchestrator return the prior result instantly. The rule of thumb: cache anything deterministic and expensive; never cache anything that touches an external service whose contents can change under you.

Retries and Backoff

Acquisition and reprojection reach across a network to a tile server and fail transiently — a dropped connection, a rate limit, a slow mirror. These tasks should retry two or three times with exponential backoff so a momentary blip does not fail the whole basin. Compute-only tasks should not retry blindly: if depression filling raises on a genuinely malformed DEM, retrying produces the same exception three times and delays the real signal. Match the retry policy to the failure mode.

Step-by-Step: Building the DAG

The build proceeds outward from a single basin to a parametrized batch.

Step 1 — Wrap each stage as a single-responsibility function. Write plain Python functions first, with no orchestrator decorators, each taking explicit input paths and returning an output path. Test them standalone. Only then wrap them as tasks.

Step 2 — Declare dependencies by data flow. Pass the output of acquire_dem into reproject, its output into fill_depressions, and so on. The orchestrator reads that chain and builds the graph; you never draw edges by hand.

Step 3 — Attach policies per task. Add retries and retry_delay_seconds to the network tasks, and a cache key function to the deterministic ones. Leave the validation task with neither, so it always runs and always raises on a bad topology.

Step 4 — Parametrize the flow. Lift the single hard-coded basin identifier into a flow parameter — a HUC code or tile name — so the same definition drives any basin. Batching then becomes mapping the flow over a list.

Step 5 — Wire observability. Emit a structured log line at the entry and exit of every task, recording input path, output path, and a cheap summary statistic (cell count, max accumulation, catchment area). These logs are what you read at 2 a.m. when a nightly batch reports one failure out of ninety.

Production-Ready Code

The following module defines the full pipeline as a Prefect flow. Each stage is a @task; the @flow chains them and returns the validated catchment path. Network tasks carry retries; deterministic tasks carry a cache key built from the input file’s content hash. The stage bodies call into the same RichDEM and pysheds operations documented across this site — here they are condensed to keep the orchestration structure legible.

import hashlib
import logging
from pathlib import Path

import geopandas as gpd
import numpy as np
import rasterio
import richdem as rd
from pysheds.grid import Grid
from prefect import flow, task, get_run_logger
from prefect.tasks import task_input_hash
from datetime import timedelta

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


def file_content_hash(context, arguments) -> str:
    """Cache key: hash the first path argument's bytes so caching keys on content,
    not just the filename. Falls back to Prefect's argument hash if no file exists."""
    for value in arguments.values():
        p = Path(str(value))
        if p.exists() and p.is_file():
            digest = hashlib.sha256(p.read_bytes()).hexdigest()
            return f"{context.task.name}-{digest}"
    return task_input_hash(context, arguments)


@task(retries=3, retry_delay_seconds=[5, 30, 120])
def acquire_dem(huc_id: str, work_dir: str) -> str:
    """Fetch the DEM tile(s) for a basin. Network-bound: retries with backoff.
    Replace the body with a real 3DEP/Copernicus request; here it validates a staged file."""
    logger = get_run_logger()
    dem_path = Path(work_dir) / huc_id / "dem_raw.tif"
    if not dem_path.exists():
        raise FileNotFoundError(f"Staged DEM missing for {huc_id}: {dem_path}")
    logger.info("Acquired DEM for %s -> %s", huc_id, dem_path)
    return str(dem_path)


@task(retries=2, retry_delay_seconds=30, cache_key_fn=file_content_hash,
      cache_expiration=timedelta(days=14))
def reproject_dem(dem_path: str, dst_epsg: int = 5070) -> str:
    """Reproject to a metric equal-area grid (EPSG:5070 Albers by default).
    Cached: deterministic in the input bytes + target CRS."""
    logger = get_run_logger()
    out_path = str(Path(dem_path).with_name("dem_proj.tif"))
    with rasterio.open(dem_path) as src:
        if src.crs.to_epsg() == dst_epsg:
            logger.info("DEM already in EPSG:%d; passing through", dst_epsg)
    # A full rasterio.warp.reproject call belongs here; see the CRS-alignment guide.
    logger.info("Reprojected %s to EPSG:%d", dem_path, dst_epsg)
    return out_path


@task(cache_key_fn=file_content_hash, cache_expiration=timedelta(days=14))
def fill_and_route(proj_path: str) -> dict:
    """Fill depressions, resolve flats, and compute D8 direction + accumulation.
    Pure compute: deterministic and cacheable, no retries needed."""
    logger = get_run_logger()
    dem = rd.LoadGDAL(proj_path)
    rd.FillDepressions(dem, in_place=True)
    rd.ResolveFlats(dem, in_place=True)
    flow_dir = rd.FlowDirection(dem, method="D8")
    flow_acc = rd.FlowAccumulation(flow_dir)
    base = Path(proj_path).with_name("")
    dir_path = str(base / "flow_dir.tif")
    acc_path = str(base / "flow_acc.tif")
    rd.SaveGDAL(dir_path, flow_dir)
    rd.SaveGDAL(acc_path, flow_acc)
    logger.info("Routed %s: max accumulation = %d cells",
                proj_path, int(np.max(np.asarray(flow_acc))))
    return {"flow_dir": dir_path, "flow_acc": acc_path}


@task
def extract_streams(routed: dict, threshold_cells: int = 1000) -> str:
    """Threshold flow accumulation into a channel mask. threshold_cells is the
    key tunable parameter; expose it as a flow argument for calibration."""
    logger = get_run_logger()
    with rasterio.open(routed["flow_acc"]) as src:
        acc = src.read(1)
        profile = src.profile
    stream_mask = (acc >= threshold_cells).astype("uint8")
    out_path = str(Path(routed["flow_acc"]).with_name("streams.tif"))
    profile.update(dtype="uint8", nodata=0, count=1)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(stream_mask, 1)
    logger.info("Stream network: %d channel cells at threshold %d",
                int(stream_mask.sum()), threshold_cells)
    return out_path


@task
def snap_and_delineate(routed: dict, stream_path: str, outlet_xy: tuple) -> str:
    """Snap the pour point to the nearest high-accumulation stream cell and
    delineate the contributing catchment with pysheds."""
    logger = get_run_logger()
    grid = Grid.from_raster(routed["flow_dir"])
    fdir = grid.read_raster(routed["flow_dir"])
    acc = grid.read_raster(routed["flow_acc"])
    x_snap, y_snap = grid.snap_to_mask(acc > 0, outlet_xy)
    catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, xytype="coordinate")
    grid.clip_to(catch)
    shapes = grid.polygonize()
    polys = [geom for geom, val in shapes if val == 1]
    gdf = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polys, crs=grid.crs.srs))
    out_path = str(Path(stream_path).with_name("catchment.gpkg"))
    gdf.dissolve().to_file(out_path, driver="GPKG")
    logger.info("Delineated catchment at snapped outlet (%.1f, %.1f)", x_snap, y_snap)
    return out_path


@task
def validate_topology(catchment_path: str) -> str:
    """Gate the run: a single valid, non-empty, simple polygon must result.
    No retries and no cache — this always runs and raises to isolate a bad basin."""
    logger = get_run_logger()
    gdf = gpd.read_file(catchment_path)
    if gdf.empty:
        raise ValueError(f"Empty catchment: {catchment_path}")
    if not gdf.geometry.is_valid.all():
        raise ValueError(f"Invalid geometry in {catchment_path}")
    if float(gdf.area.sum()) <= 0:
        raise ValueError(f"Zero-area catchment: {catchment_path}")
    logger.info("Topology OK: %d polygon(s), %.2f km2",
                len(gdf), float(gdf.area.sum()) / 1e6)
    return catchment_path


@flow(name="dem-to-watershed")
def dem_to_watershed(huc_id: str, outlet_xy: tuple, work_dir: str = "./work",
                     threshold_cells: int = 1000) -> str:
    """Orchestrate the full chain for one basin and return the validated catchment path."""
    logger = get_run_logger()
    logger.info("Starting pipeline for HUC %s", huc_id)
    dem_raw = acquire_dem(huc_id, work_dir)
    dem_proj = reproject_dem(dem_raw)
    routed = fill_and_route(dem_proj)
    streams = extract_streams(routed, threshold_cells=threshold_cells)
    catchment = snap_and_delineate(routed, streams, outlet_xy)
    validated = validate_topology(catchment)
    logger.info("Pipeline complete for HUC %s -> %s", huc_id, validated)
    return validated


if __name__ == "__main__":
    dem_to_watershed(
        huc_id="03010101",
        outlet_xy=(1_450_000.0, 1_920_000.0),
        threshold_cells=1500,
    )

The flow reads top to bottom as ordinary Python, yet Prefect executes it as a graph: fill_and_route will not start until reproject_dem returns, and on rerun the cached fill result is returned without recomputation whenever the reprojected DEM’s bytes are unchanged.

Parametrization Over Tiles and HUCs

A production system rarely delineates one basin. Lift the basin identifier into a flow parameter — done above — and batching reduces to iterating a manifest of (huc_id, outlet_xy) pairs. In Prefect, .map() fans the work across a task runner for concurrency; in Airflow, dynamic task mapping with .expand does the same on the scheduler. When basins nest or share tile boundaries, coordinate the partition up front with the strategies in basin partitioning strategies so two flows never contend for the same intermediate file.

import logging
from prefect import flow, get_run_logger

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

BASIN_MANIFEST = [
    {"huc_id": "03010101", "outlet_xy": (1_450_000.0, 1_920_000.0)},
    {"huc_id": "03010102", "outlet_xy": (1_462_500.0, 1_905_000.0)},
    {"huc_id": "03010103", "outlet_xy": (1_478_000.0, 1_888_500.0)},
]


@flow(name="batch-dem-to-watershed")
def batch_delineate(manifest: list[dict], threshold_cells: int = 1000) -> list[str]:
    """Run the per-basin flow for every entry, isolating failures per basin."""
    logger = get_run_logger()
    results = []
    for basin in manifest:
        try:
            path = dem_to_watershed(
                huc_id=basin["huc_id"],
                outlet_xy=basin["outlet_xy"],
                threshold_cells=threshold_cells,
            )
            results.append(path)
        except Exception as exc:  # isolate: one bad basin must not abort the batch
            logger.error("Basin %s failed: %s", basin["huc_id"], exc)
    logger.info("Batch complete: %d/%d basins succeeded",
                len(results), len(manifest))
    return results


if __name__ == "__main__":
    batch_delineate(BASIN_MANIFEST, threshold_cells=1500)

Wrapping each subflow in a try/except is the crudest form of failure isolation and the most portable. Prefect and Airflow both offer richer state-based isolation, but the principle is identical: the batch’s success is the count of basins that passed, not an all-or-nothing gate.

Observability and Logging

Structured, per-task logs are the difference between “the pipeline failed” and “basin 03010107 failed at snap_and_delineate because the pour point fell outside the DEM extent.” Use get_run_logger() inside tasks so log lines are tagged with the flow run and task name and surface in the orchestrator UI. Log a cheap summary statistic at each stage exit — cell counts, max accumulation, catchment area — because those numbers catch silent corruption that an exception never would. A catchment of 4 square metres does not raise; it just looks wrong in the log, and the topology gate can turn that suspicion into a hard failure.

Persist intermediate artifacts under a per-basin working directory keyed by HUC. When a basin fails, the partial outputs are still on disk for inspection, and a restart resumes from the last completed stage rather than re-downloading the DEM.

Validation

The pipeline’s final task is a gate, not a formality. A delineation can complete without raising and still be wrong: a catchment split into slivers by a nodata seam, a polygon of near-zero area from a mis-snapped outlet, a self-intersecting boundary from a polygonization artifact. The validate_topology task rejects the empty, the invalid, and the zero-area, and a production version should also check that the delineated area is within a plausible range of the expected basin size. For the full battery of geometric checks — sliver detection, dangling edges, shared-boundary consistency between adjacent catchments — route the output through the checks in boundary topology validation before publishing. Gating here means a corrupt catchment never enters the cache and never propagates to a downstream consumer.

Common Failure Modes

Partial failure and restart. A task dies mid-write and the flow reports failure. On restart, an idempotent pipeline resumes from the last completed stage using cached upstream results; a non-idempotent one re-runs everything or, worse, trips over its own half-written files. The fix is deterministic output paths plus write-to-temp-then-rename, so an interrupted task leaves either the old valid artifact or nothing — never a truncated GeoTIFF.

Tile-boundary state. When a basin spans multiple DEM tiles, flow paths cross seams. If tasks share mutable state across tiles — a common accumulation grid, an appended stream vector — parallel basins corrupt each other. Give every basin an isolated working directory and buffer tiles with enough overlap that no flow path is severed at a seam, mirroring the edge-handling discipline in the routing stage.

Non-idempotent writes. Appending to a shared GeoPackage layer, or writing catchments to a single file keyed only by a timestamp, means a rerun duplicates or clobbers records. Key every output on the basin identifier and content, and treat each run as producing a fresh, replaceable artifact rather than mutating an accumulating one.

Scheduler backpressure. Fan a hundred basins across a task runner at once and each spawns memory-hungry RichDEM arrays; the machine thrashes and tasks time out, which the orchestrator reads as failure and retries, deepening the overload. Cap concurrency explicitly — a bounded task runner in Prefect, max_active_tasks in Airflow — sized to available RAM rather than core count, since a single large DEM can hold hundreds of megabytes resident.

Cache keys on mutable inputs. Caching on a filename rather than file content returns a stale result when the underlying tile is silently replaced. Key deterministic tasks on a content hash, as the file_content_hash helper does, so a changed input invalidates the cache automatically.

Related Topics

• Building a Prefect flow for DEM-to-watershed automation — the Prefect 2.x implementation with task caching, retries, and .map() concurrency
• Airflow DAG templates for batch watershed delineation — the scheduled TaskFlow DAG with dynamic task mapping and an input-DEM sensor
• Outlet Point Mapping & Validation — snapping and verifying the pour points the delineation stage depends on
• Basin Partitioning Strategies — organizing tiles and sub-basins so parallel flows do not contend
• Boundary Topology Validation — the geometric checks that gate a delineated catchment before publishing
• Watershed Delineation & Catchment Synchronization — the section overview this pipeline automates end to end