Hydrology Data Preparation & DEM Processing

Digital Elevation Models are the computational substrate of every watershed model. The algorithms that delineate catchments, route stormflow, and calibrate rainfall-runoff parameters all inherit whatever errors, datum shifts, and conditioning choices were baked into the elevation surface at the start. Errors propagate silently: a 1 m vertical offset from a datum mismatch can relocate a drainage divide; an unfilled depression can truncate an accumulation matrix mid-basin; a poorly chosen resampling kernel can smooth away the ridge that separates two sub-catchments. Getting this foundation right is not a preliminary step — it is the work.

This guide covers the full preparation stack as part of the watershed modeling workflow hosted at this site, from SRTM and LiDAR data acquisition through coordinate reference system alignment, DEM pit filling algorithms, and spatial resolution tradeoffs. The conditioned DEM produced here is the direct input to flow routing and stream network extraction, so defects at this stage surface as routing failures downstream.

Foundational Concepts

What a DEM actually encodes

A rasterized DEM stores one elevation value per grid cell. That value is not a precise ground measurement — it is an aggregated, interpolated sample from a sensor that measured something (radar backscatter, LiDAR pulse returns, stereo imagery parallax) and then went through a processing chain before you downloaded it. Understanding what your source sensor captured is prerequisite to deciding which conditioning steps are required.

LiDAR point clouds are classified into ground, vegetation, buildings, and noise returns. A bare-earth DEM uses only ground returns; a Digital Surface Model (DSM) includes canopy tops. For hydrology, bare-earth is almost always the target — flow routing over a DSM will terminate at rooftops and tree canopies. Classification quality degrades in dense low vegetation and on high slopes, introducing isolated high-noise returns that appear as micro-ridges in the DEM.

Radar-derived DEMs (SRTM, TanDEM-X) measure phase difference between two synthetic aperture radar passes. Vegetation canopy adds a positive elevation bias because radar pulses often reflect from canopy rather than ground — measured globally at 2–8 m in forested zones. This bias shifts flow accumulation thresholds and artificially raises ridge elevations relative to valley floors.

Photogrammetric DEMs from stereo imagery require texture and contrast to find correspondences. Water bodies, snow, and featureless sand produce voids or high noise. In semi-arid regions where watershed models are common, dry lake beds and alluvial fans can arrive with significant interpolation artifacts.

Flow physics and terrain topology

Water flows from higher to lower elevation, following the steepest descent gradient. At the grid scale, this translates to comparing each cell’s elevation against its eight neighbors (cardinal and diagonal) and assigning flow to the lowest one — the core principle behind D8 flow direction. The same terrain topology underlies catchment delineation: tracing all cells that drain to a common outlet defines a watershed.

Three topological features disrupt this simple model and must be resolved before routing:

• Depressions (sinks): Local minima with no outflow neighbor below them. Accumulation terminates here rather than propagating downstream.
• Flat areas: Regions where multiple adjacent cells share the same elevation, producing a gradient of zero and no preferred flow direction.
• Ridges and saddles: High-elevation boundaries between drainage basins. Sub-meter errors here can shift watershed boundaries by kilometers when basins are elongated and low-relief.

Resolving these three features is the core task of hydrological conditioning. The approach chosen — filling, breaching, or a hybrid — directly affects the accuracy of every derived product downstream.

Coordinate reference system geometry

Every DEM operation involving distances, slopes, or areas requires a projected coordinate system. Working in geographic coordinates (degrees of latitude/longitude) produces incorrect slope calculations because a degree of longitude is not a fixed distance — it narrows from ~111 km at the equator to zero at the poles. Incorrect slope magnitudes bias flow direction and accumulation.

Vertical datum alignment is equally critical. Two tiles downloaded from different sources may be in NAD83 horizontally but use different vertical datums (NAVD88 vs EGM2008 vs EGM96). The offset between EGM96 and EGM2008 ranges from −0.5 to +1.5 m globally. In low-relief coastal plains, a 1 m vertical shift across a tile boundary can invert the flow direction at the seam and break drainage continuity across the entire mosaic. Full strategies for datum alignment are covered in coordinate reference system alignment.

Algorithm Landscape

Depression handling approaches

The decision between filling and breaching depressions is not purely algorithmic — it depends on terrain type, data source, and what errors you are more willing to tolerate. The diagram below shows the decision path:

The priority-flood algorithm implemented in richdem (FillDepressions) is the most commonly deployed because it handles large grids efficiently via a priority queue. The whitebox library’s BreachDepressionsLeastCost function implements the Lindsay breaching approach. For production pipelines, combining both — fill depressions shallower than 0.5 m (likely sensor noise), then breach the remainder — reduces both artificial storage inflation and artificial channel carving. The full mathematical basis and implementation tradeoffs are covered in DEM pit filling algorithms.

Flow direction encoding

Once the DEM is conditioned, each cell is assigned a flow direction encoding. The choice of encoding affects stream network topology, flow accumulation magnitudes, and the ability to represent divergent flow on hillslopes.

For the conditioning stage specifically, D8 is the standard choice because it produces a single-valued flow direction grid compatible with accumulation algorithms in richdem, whitebox, and TauDEM. Transition to D-Infinity routing patterns after conditioning when modeling divergent hillslope processes or comparing algorithms on steep terrain. For flat agricultural plains, multiple flow direction methods avoid the artificial convergence that D8 introduces across shallow gradients.

Resampling methods for multi-resolution data

For hydrologic use, bilinear is the default for DEM resampling within ±2× the source resolution. When aggregating LiDAR (1 m) to a 10 m product, use average resampling to preserve the statistical representation of each coarser cell rather than sampling a single source pixel. The interaction between resampling choice and watershed delineation accuracy is analyzed in depth in spatial resolution tradeoffs.

Production Python Architecture

The pipeline below handles one DEM tile end-to-end: validation, CRS harmonization, hydrological conditioning, derivative generation, and export. It uses rasterio, richdem, pyproj, and standard library modules. Scale to multi-tile extents by wrapping process_dem_pipeline in a Dask delayed graph or a Prefect flow.

import logging
import hashlib
from pathlib import Path

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
import richdem as rd
from pyproj import CRS, Transformer

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s %(name)s — %(message)s",
)
log = logging.getLogger("dem_pipeline")


def _checksum(path: Path) -> str:
    """Return SHA-256 hex digest of a file for immutable provenance logging."""
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        for chunk in iter(lambda: fh.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def validate_dem(src_path: Path) -> dict:
    """
    Open a DEM tile and verify minimum metadata requirements.
    Returns a metadata dict or raises ValueError on failure.
    """
    log.info("Validating DEM: %s", src_path)
    with rasterio.open(src_path) as src:
        if src.crs is None:
            raise ValueError(f"{src_path}: missing CRS — cannot continue.")
        if src.nodata is None:
            log.warning("%s: nodata value absent; defaulting to -9999.0", src_path)
        nodata_val = src.nodata if src.nodata is not None else -9999.0

        dem = src.read(1).astype(np.float32)
        void_fraction = np.mean(dem == nodata_val)
        if void_fraction > 0.05:
            raise ValueError(
                f"{src_path}: {void_fraction:.1%} void cells — exceeds 5% threshold."
            )

        meta = {
            "path": str(src_path),
            "sha256": _checksum(src_path),
            "crs_wkt": src.crs.to_wkt(),
            "epsg": src.crs.to_epsg(),
            "nodata": nodata_val,
            "shape": src.shape,
            "res_m": src.res,
            "void_fraction": float(void_fraction),
        }
    log.info("Validation passed — EPSG:%s, shape=%s, voids=%.2f%%",
             meta["epsg"], meta["shape"], meta["void_fraction"] * 100)
    return meta


def reproject_dem(
    src_path: Path,
    dst_path: Path,
    target_epsg: int,
    resampling: Resampling = Resampling.bilinear,
) -> None:
    """
    Reproject DEM to target_epsg using bilinear resampling.
    Logs source and target CRS and the PROJ transformation string.
    """
    log.info("Reprojecting %s → EPSG:%d", src_path.name, target_epsg)
    target_crs = CRS.from_epsg(target_epsg)

    with rasterio.open(src_path) as src:
        if src.crs.to_epsg() == target_epsg:
            log.info("CRS already matches target — skipping reproject.")
            import shutil
            shutil.copy2(src_path, dst_path)
            return

        transformer = Transformer.from_crs(src.crs, target_crs, always_xy=True)
        log.info("PROJ pipeline: %s", transformer.to_proj4())

        transform, width, height = calculate_default_transform(
            src.crs, target_crs, src.width, src.height, *src.bounds
        )
        profile = src.profile.copy()
        profile.update(
            crs=target_crs,
            transform=transform,
            width=width,
            height=height,
            dtype="float32",
        )

        with rasterio.open(dst_path, "w", **profile) as dst:
            reproject(
                source=rasterio.band(src, 1),
                destination=rasterio.band(dst, 1),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=target_crs,
                resampling=resampling,
            )
    log.info("Reprojection complete → %s", dst_path)


def condition_dem(
    src_path: Path,
    dst_path: Path,
    nodata: float = -9999.0,
    fill_depth_threshold: float = 0.5,
) -> dict:
    """
    Hydrological conditioning via hybrid fill-then-breach.
    - Depressions shallower than fill_depth_threshold metres are filled.
    - Deeper depressions are left for downstream breaching (whitebox or TauDEM).
    Returns a dict of conditioning diagnostics.
    """
    log.info("Conditioning DEM: %s (fill threshold=%.2f m)", src_path.name, fill_depth_threshold)

    with rasterio.open(src_path) as src:
        dem_array = src.read(1).astype(np.float64)
        profile = src.profile.copy()
        nodata = src.nodata if src.nodata is not None else nodata

    rd_dem = rd.rdarray(dem_array, no_data=nodata)

    # Identify depressions before filling for diagnostics
    depression_depth = rd.FillDepressions(rd_dem, in_place=False) - rd_dem
    shallow_mask = (depression_depth > 0) & (depression_depth <= fill_depth_threshold)
    n_shallow = int(np.sum(shallow_mask))
    n_deep = int(np.sum(depression_depth > fill_depth_threshold))

    log.info("Depressions found — shallow (≤%.2f m): %d cells, deep: %d cells",
             fill_depth_threshold, n_shallow, n_deep)

    # Apply priority-flood fill only to shallow depressions
    filled = rd_dem.copy()
    filled[shallow_mask] = (rd_dem + depression_depth)[shallow_mask]
    filled_rd = rd.rdarray(np.array(filled), no_data=nodata)

    max_delta = float(np.nanmax(depression_depth[shallow_mask])) if n_shallow > 0 else 0.0
    log.info("Max elevation raised: %.4f m over %d cells", max_delta, n_shallow)

    profile.update(dtype="float32", nodata=nodata)
    dst_path.parent.mkdir(parents=True, exist_ok=True)
    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(np.array(filled_rd).astype("float32"), 1)

    log.info("Conditioned DEM written → %s", dst_path)
    return {
        "n_shallow_filled": n_shallow,
        "n_deep_remaining": n_deep,
        "max_fill_delta_m": max_delta,
    }


def derive_flow_products(
    conditioned_path: Path,
    output_dir: Path,
    method: str = "D8",
    stream_threshold: int = 1000,
) -> None:
    """
    Derive flow direction, flow accumulation, and a binary stream mask
    from a conditioned DEM. Writes three GeoTIFF outputs.
    """
    log.info("Deriving flow products from %s (method=%s, threshold=%d)",
             conditioned_path.name, method, stream_threshold)

    with rasterio.open(conditioned_path) as src:
        dem_array = src.read(1).astype(np.float64)
        profile = src.profile.copy()
        nodata = src.nodata if src.nodata is not None else -9999.0

    rd_dem = rd.rdarray(dem_array, no_data=nodata)
    flow_dir = rd.FlowDirection(rd_dem, method=method)
    flow_acc = rd.FlowAccumulation(rd_dem, method=method)

    stream_mask = (np.array(flow_acc) >= stream_threshold).astype(np.int16)
    n_stream_cells = int(np.sum(stream_mask))
    log.info("Stream cells at threshold %d: %d", stream_threshold, n_stream_cells)

    output_dir.mkdir(parents=True, exist_ok=True)
    int16_profile = {**profile, "dtype": "int16", "nodata": -32768}
    float32_profile = {**profile, "dtype": "float32"}

    for arr, fname, prof in [
        (np.array(flow_dir).astype("int16"),  "flow_direction.tif",  int16_profile),
        (np.array(flow_acc).astype("float32"), "flow_accumulation.tif", float32_profile),
        (stream_mask,                           "stream_mask.tif",     int16_profile),
    ]:
        with rasterio.open(output_dir / fname, "w", **prof) as dst:
            dst.write(arr, 1)
        log.info("Written: %s", fname)


def process_dem_pipeline(
    raw_path: Path,
    work_dir: Path,
    target_epsg: int = 32614,
    stream_threshold: int = 1000,
) -> None:
    """
    Orchestrate the full DEM preparation pipeline for one tile.
    Stages: validate → reproject → condition → derive → done.
    """
    log.info("=== Pipeline start: %s ===", raw_path.name)
    work_dir.mkdir(parents=True, exist_ok=True)

    meta = validate_dem(raw_path)

    reproj_path = work_dir / "reprojected.tif"
    reproject_dem(raw_path, reproj_path, target_epsg=target_epsg)

    cond_path = work_dir / "conditioned.tif"
    diag = condition_dem(reproj_path, cond_path)
    log.info("Conditioning diagnostics: %s", diag)

    derive_flow_products(cond_path, work_dir / "derivatives", stream_threshold=stream_threshold)

    log.info("=== Pipeline complete: outputs in %s ===", work_dir)

Key engineering invariants enforced by this implementation:

• Idempotency: Checksums on raw inputs detect tile changes between runs. Intermediate files can be cached and skipped when their source is unchanged.
• Typed nodata handling: Every array operation preserves the nodata value explicitly; no silent NaN substitution.
• Structured diagnostics: condition_dem returns a dict of cell counts and delta magnitudes — feed these into a monitoring dashboard to detect anomalous conditioning across tile batches.
• Separation of concerns: Each function has one responsibility and can be unit-tested independently.

Scaling & Cloud-Native Workflows

Tile-based processing with overlap buffers

Large DEMs must be tiled for parallel processing. The critical constraint is that hydrological conditioning and flow accumulation are not locally bounded — a depression on tile A may drain through tile B. Without overlap buffers, tile boundaries generate false depressions and interrupted flow paths.

The standard approach uses a buffer of at least the diameter of the largest expected depression (a reasonable default is 200 cells). Each tile is processed with its buffer, then the buffer region is cropped from the output before mosaicking. rasterio.windows provides the windowed read/write API; Dask’s dask.array enables delayed, chunked computation across the tiled grid.

import dask
import dask.array as da
import rasterio
from rasterio.windows import Window
import numpy as np
import logging

log = logging.getLogger("dem_tiled")


def tile_paths(src_path, tile_size=2048, buffer=200):
    """Yield (window_with_buffer, crop_window) pairs for a raster."""
    with rasterio.open(src_path) as src:
        rows, cols = src.shape
    for r in range(0, rows, tile_size):
        for c in range(0, cols, tile_size):
            r0 = max(0, r - buffer)
            c0 = max(0, c - buffer)
            r1 = min(rows, r + tile_size + buffer)
            c1 = min(cols, c + tile_size + buffer)
            buffered = Window(c0, r0, c1 - c0, r1 - r0)
            # Inner crop after processing
            inner_r = r - r0
            inner_c = c - c0
            inner_h = min(tile_size, rows - r)
            inner_w = min(tile_size, cols - c)
            inner = Window(inner_c, inner_r, inner_w, inner_h)
            yield buffered, inner


@dask.delayed
def process_tile(src_path, buffered_window, inner_window, target_epsg, work_dir, tile_id):
    """Process one buffered tile and return the path to the cropped output."""
    import richdem as rd
    log.info("Tile %s: starting", tile_id)
    with rasterio.open(src_path) as src:
        arr = src.read(1, window=buffered_window).astype(np.float64)
        nodata = src.nodata or -9999.0
        transform = src.window_transform(buffered_window)
        profile = src.profile.copy()
        profile.update(
            width=buffered_window.width,
            height=buffered_window.height,
            transform=transform,
            dtype="float32",
        )

    rd_dem = rd.rdarray(arr, no_data=nodata)
    filled = rd.FillDepressions(rd_dem, in_place=False)
    out_arr = np.array(filled).astype("float32")

    # Crop buffer zone before writing
    cropped = out_arr[
        inner_window.row_off : inner_window.row_off + inner_window.height,
        inner_window.col_off : inner_window.col_off + inner_window.width,
    ]
    out_path = work_dir / f"tile_{tile_id}_conditioned.tif"
    profile.update(width=inner_window.width, height=inner_window.height)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(cropped, 1)
    log.info("Tile %s: written %s", tile_id, out_path)
    return out_path

COG and Zarr output strategies

Cloud-Optimized GeoTIFF (COG) format is the standard for cloud storage of raster derivatives. Writing COGs requires tiled internal layout and appropriate DEFLATE compression. Pass tiled=True, blockxsize=512, blockysize=512, compress="deflate" in the rasterio write profile. For analysis workflows that read sub-regions repeatedly (e.g., iterative model calibration), Zarr with chunked coordinates offers lower-latency random access than COG because chunks map directly to object storage keys.

For continental-scale pipelines, AWS Batch, Google Cloud Run, or Kubernetes Jobs can execute process_tile for each tile in parallel. The Dask distributed scheduler coordinates dependencies, handles retries on node failure, and provides a real-time task graph dashboard.

Memory estimation

Rule of thumb for a single richdem fill operation: peak RAM ≈ 6–8× the raw array size. A 10,000 × 10,000 float64 DEM occupies 800 MB; FillDepressions peaks at 5–6 GB. Tile to 2,000 × 2,000 cells plus a 200-cell buffer to keep peak RAM under 8 GB per worker on standard cloud instances.

Validation & QA/QC

A conditioned DEM is only trustworthy if it passes independent checks. Run these after every conditioning and derivation run:

Topological checks:

• Flow accumulation maximum equals the number of non-nodata cells (all cells drain to the outlet)
• Zero isolated accumulation values interior to the basin (indicates unfilled depressions)
• Stream network is fully connected from headwaters to outlet (no orphan segments)
• Accumulation-derived watershed area matches reference area within 2%

Elevation fidelity:

• Difference raster (conditioned minus raw) has 95th-percentile delta < 1.0 m
• Maximum elevation change does not exceed the known maximum depression depth in the region
• Mean elevation difference < 0.05 m (large positive bias indicates overfilling)

CRS and geometry:

• Output CRS matches target EPSG (check with rasterio.open().crs.to_epsg())
• Cell size matches expected resolution ± 0.01 m
• Bounding box overlaps expected catchment extent

Comparison against authoritative hydrography:

• Derived stream network overlaps the National Hydrography Dataset (NHD) or HydroSHEDS within 1–2 cell widths
• Watershed boundaries match HUC boundaries within acceptable tolerance for the DEM resolution
• Strahler stream order distribution matches expected regional pattern

import numpy as np
import rasterio
import logging

log = logging.getLogger("dem_qa")


def run_qa_checks(raw_path, conditioned_path, flow_acc_path, expected_area_km2):
    """
    Run post-conditioning QA checks. Logs PASS/FAIL for each criterion.
    Returns True if all checks pass, False otherwise.
    """
    results = {}

    with rasterio.open(raw_path) as raw, rasterio.open(conditioned_path) as cond:
        raw_arr = raw.read(1).astype(np.float32)
        cond_arr = cond.read(1).astype(np.float32)
        nodata = cond.nodata or -9999.0
        valid = cond_arr != nodata
        delta = cond_arr - raw_arr
        delta_valid = delta[valid]

    p95_delta = float(np.percentile(delta_valid, 95))
    mean_delta = float(np.mean(delta_valid))
    results["p95_elevation_delta_lt_1m"] = p95_delta < 1.0
    results["mean_elevation_delta_lt_0.05m"] = abs(mean_delta) < 0.05
    log.info("Elevation delta P95=%.3f m, mean=%.4f m", p95_delta, mean_delta)

    with rasterio.open(flow_acc_path) as acc_src:
        acc_arr = acc_src.read(1)
        cell_area_m2 = abs(acc_src.res[0] * acc_src.res[1])
        max_acc = int(np.nanmax(acc_arr))
        n_valid = int(np.sum(acc_arr >= 0))
        derived_area_km2 = max_acc * cell_area_m2 / 1e6

    results["flow_acc_max_equals_n_cells"] = abs(max_acc - n_valid) / n_valid < 0.02
    area_error_pct = abs(derived_area_km2 - expected_area_km2) / expected_area_km2 * 100
    results["area_within_2pct"] = area_error_pct < 2.0
    log.info("Derived area=%.2f km², reference=%.2f km², error=%.1f%%",
             derived_area_km2, expected_area_km2, area_error_pct)

    all_pass = all(results.values())
    for check, passed in results.items():
        log.info("QA %s: %s", check, "PASS" if passed else "FAIL")
    log.info("Overall QA: %s", "PASS" if all_pass else "FAIL")
    return all_pass

Common Failure Modes

Memory exhaustion on large DEMs. richdem loads the full array into memory before conditioning. A 30 m DEM covering a 50,000 km² basin at float64 occupies ~6 GB; FillDepressions peaks at 40 GB. Mitigation: tile to 2,000 × 2,000 cell blocks with 200-cell buffers, or switch to the whitebox streaming mode which processes tiles internally.

Flat-area stagnation. After filling depressions, the raised cells form large plateaus where all neighbors share the same elevation. D8 assignment fails (no downslope neighbor) and flow accumulation stalls. Mitigation: apply richdem’s ResolveFlats step after FillDepressions. This adds a tiny gradient (epsilon slope) directed toward the plateau’s drainage outlet without altering visible terrain structure. See also removing flat area artifacts from flow direction grids for algorithm-specific strategies.

Projection artifacts at tile seams. When mosaicking tiles reprojected independently, minor floating-point differences in the transform matrix can produce 1-pixel gaps or overlaps at seam lines. Mitigation: reproject all tiles using a shared target transform computed once from the full mosaic bounds (rasterio.merge.merge with a pre-computed output transform).

Vertical datum offset after tile merge. Mixing NAVD88 (CONUS) and EGM2008 (global) tiles without explicit vertical transformation introduces systematic offsets at international or coastal boundaries. Mitigation: use pyproj compound CRS definitions that include the vertical component, and apply the VERTCON or GEOID18 grid shift before merging. The full workflow is detailed in fixing CRS mismatches in watershed shapefiles.

Radar canopy bias propagating into stream initiation. SRTM-derived DEMs in forested areas have 2–8 m positive bias on ridges. This compresses the apparent relief, raises flow accumulation thresholds, and under-generates headwater channels. Mitigation: apply a global or regional canopy-height correction raster (e.g., the ETH Global Canopy Height model) before conditioning, or calibrate accumulation thresholds against observed stream network extent. Threshold calibration strategies are covered in stream threshold tuning.

Nodata boundary truncation. DEMs delivered as rectangular tiles often have nodata cells at irregular coastal or administrative boundaries. If FillDepressions treats nodata as an infinitely deep sink, it will attempt to drain toward it, generating spurious channels along the nodata margin. Mitigation: mask nodata regions explicitly and use richdem’s ignore_nodata flag, or set boundary cells to a high sentinel value before conditioning.

Over-smoothing from resampling. Cubic convolution resampling of a LiDAR DEM from 1 m to 10 m can suppress micro-ridges that separate adjacent headwater basins. In urban or agricultural landscapes with berms and levees, this merges sub-catchments that should remain distinct. Mitigation: use average resampling when aggregating, and validate post-resampling basin boundaries against the original-resolution delineation.

Frequently Asked Questions

Should I fill or breach depressions in my DEM?

Use a hybrid approach: fill small noise-induced depressions (depth < 0.5 m, area < 10 cells) and breach larger, topographically significant depressions. Filling alone overestimates ponded storage and creates flat plateaus that require secondary flat-resolution steps. Breaching alone risks carving artificial channels through legitimate saddles. The whitebox function BreachDepressionsLeastCost with a max_dist parameter limits carve length, preventing implausibly long artificial channels.

What spatial resolution should I use for watershed delineation?

Match resolution to your smallest sub-basin and your modeling objective. Floodplain or urban drainage design requires 1–3 m LiDAR DEMs. Regional HEC-HMS or SWAT calibration typically operates at 10–30 m. As a rule of thumb, the DEM cell size should be no larger than 1/10th the width of your smallest channel of interest. Avoid resampling coarser than this threshold — see spatial resolution tradeoffs for the full analysis including the effect on nested sub-basin delineation.

How do I handle CRS mismatches when merging DEM tiles?

Reproject all tiles to a single metric CRS before mosaicking. Use pyproj compound CRS definitions to handle both horizontal datum and vertical datum shifts in a single transformation. Log the exact EPSG codes and PROJ transformation strings — these form part of your audit trail for regulatory submissions and peer review. Full strategies for datum alignment are covered in coordinate reference system alignment.

Can I run DEM conditioning in parallel across multiple tiles?

Yes, but tile boundaries require overlap buffers. Process each tile with a buffer of at least 200 cells on each side, apply conditioning, then crop the buffer before mosaicking. Without buffers, depressions that straddle tile boundaries remain unfilled and break flow accumulation continuity. Wrap the per-tile function in dask.delayed or a Prefect task to orchestrate parallel execution.

How do I validate that my conditioned DEM is hydrologically correct?

Three independent checks: (1) the maximum flow accumulation value should approximately equal the number of valid cells — if it is significantly lower, unfilled depressions are terminating flow paths; (2) the derived stream network should overlay the National Hydrography Dataset within 1–2 cell widths at your DEM resolution; (3) the difference raster (conditioned minus raw) should have a 95th-percentile delta under 1 m — larger values indicate overfilling or datum errors.

Related Topics

• SRTM and LiDAR Data Acquisition — sourcing, validating, and structuring elevation datasets for Python pipelines
• Coordinate Reference System Alignment — horizontal and vertical datum harmonization strategies
• DEM Pit Filling Algorithms — mathematical foundations and implementation tradeoffs for depression handling
• Spatial Resolution Tradeoffs — how grid scaling decisions affect hydrologic representation and computational cost
• Flow Routing Algorithms & Stream Network Extraction — D8, D-Infinity, and MFD routing on the conditioned DEM produced here
• Watershed Delineation & Catchment Synchronization — translating conditioned DEMs and flow products into delineated basin boundaries