Spatial Resolution Tradeoffs in Automated Watershed Modeling

conda install -c conda-forge rasterio numpy scipy richdem geopandas scikit-image matplotlib

import logging
import rasterio
from pathlib import Path

logger = logging.getLogger(__name__)

def inspect_dem(path: str) -> dict:
    """Extract resolution, CRS, dtype, and fill status from a source DEM."""
    dem_path = Path(path)
    if not dem_path.exists():
        raise FileNotFoundError(f"DEM not found: {dem_path}")

    with rasterio.open(dem_path) as src:
        if not src.crs:
            raise ValueError(f"DEM lacks CRS metadata: {dem_path}")

        meta = {
            "native_res_m": src.res[0],
            "crs": str(src.crs),
            "is_projected": src.crs.is_projected,
            "dtype": src.dtypes[0],
            "bounds": src.bounds,
            "nodata": src.nodata,
            "shape": src.shape,
            "cell_count": src.width * src.height,
        }

    logger.info(
        "DEM inspection complete: %.1f m resolution, %s, %d x %d cells",
        meta["native_res_m"], meta["crs"], meta["shape"][0], meta["shape"][1],
    )
    if not meta["is_projected"]:
        logger.warning(
            "DEM CRS is geographic (degrees), not projected — "
            "resample in a projected CRS to preserve metric distances."
        )
    return meta

import logging
import numpy as np
import richdem as rd

logger = logging.getLogger(__name__)

def flow_accumulation_profile(dem_path: str) -> dict:
    """
    Compute D8 flow accumulation and return key distribution metrics.
    Used to compare behaviour across candidate resolutions.
    """
    logger.info("Computing flow accumulation profile for: %s", dem_path)

    dem = rd.LoadGDAL(dem_path)
    filled = rd.FillDepressions(dem, in_place=False)
    fa = rd.FlowAccumulation(filled, method="D8")

    fa_arr = np.array(fa, dtype=np.float64)
    fa_arr[fa_arr <= 0] = np.nan  # exclude non-contributing cells

    # Stream density proxy: cells with accumulation above 1000-cell threshold
    stream_cell_pct = float(np.nansum(fa_arr > 1000) / np.isfinite(fa_arr).sum() * 100)

    profile = {
        "p50": float(np.nanpercentile(fa_arr, 50)),
        "p95": float(np.nanpercentile(fa_arr, 95)),
        "p99": float(np.nanpercentile(fa_arr, 99)),
        "stream_density_pct": stream_cell_pct,
    }

    logger.info(
        "Accumulation profile — p50: %.0f, p95: %.0f, p99: %.0f, stream density: %.2f%%",
        profile["p50"], profile["p95"], profile["p99"], profile["stream_density_pct"],
    )
    return profile

import logging
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

logger = logging.getLogger(__name__)

def resample_dem(
    input_path: str,
    output_path: str,
    target_res_m: float,
    method: Resampling = Resampling.average,
) -> dict:
    """
    Resample a DEM to a new grid spacing.
    Defaults to area-weighted average — appropriate for coarsening.
    For refinement, pass Resampling.cubic.
    """
    logger.info("Resampling %s to %.1f m using %s", input_path, target_res_m, method.name)

    with rasterio.open(input_path) as src:
        if not src.crs.is_projected:
            raise ValueError(
                "Resampling requires a projected CRS. "
                "Reproject to UTM before calling resample_dem()."
            )

        transform, width, height = calculate_default_transform(
            src.crs, src.crs, src.width, src.height, *src.bounds,
            resolution=target_res_m,
        )

        kwargs = src.meta.copy()
        kwargs.update({"transform": transform, "width": width, "height": height})

        with rasterio.open(output_path, "w", **kwargs) as dst:
            for band_idx in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band_idx),
                    destination=rasterio.band(dst, band_idx),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=src.crs,
                    resampling=method,
                )

    logger.info("Resampled DEM written to %s (%d x %d cells)", output_path, height, width)
    return {"output_path": output_path, "width": width, "height": height, "res_m": target_res_m}

import logging
import numpy as np
import richdem as rd

logger = logging.getLogger(__name__)

def condition_resampled_dem(dem_path: str) -> np.ndarray:
    """
    Fill depressions and compute D8 flow accumulation on a resampled DEM.
    Returns the flow-accumulation array for downstream threshold tuning.
    """
    logger.info("Conditioning resampled DEM: %s", dem_path)

    dem = rd.LoadGDAL(dem_path)
    pre_fill_stats = {"min": float(np.min(dem)), "max": float(np.max(dem))}
    logger.debug("Pre-fill elevation range: %.2f – %.2f m", pre_fill_stats["min"], pre_fill_stats["max"])

    filled = rd.FillDepressions(dem, in_place=False)
    fill_diff = np.array(filled) - np.array(dem)
    n_filled = int(np.sum(fill_diff > 0.001))  # cells raised > 1 mm
    logger.info("Pit-filling raised %d cells (%.2f%% of grid)", n_filled, n_filled / dem.size * 100)

    fa = rd.FlowAccumulation(filled, method="D8")
    fa_arr = np.array(fa, dtype=np.float64)
    fa_arr[fa_arr <= 0] = np.nan

    logger.info(
        "Post-conditioning flow accumulation — p95: %.0f, p99: %.0f",
        np.nanpercentile(fa_arr, 95), np.nanpercentile(fa_arr, 99),
    )
    return fa_arr

import logging
from pathlib import Path

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

logger = logging.getLogger(__name__)


def build_resolution_calibrated_dem(
    source_path: str,
    output_path: str,
    target_res_m: float,
    coarsen: bool = True,
) -> dict:
    """
    Full resolution-calibration pipeline:
      1. Inspect source DEM metadata.
      2. Resample to target grid spacing.
      3. Fill depressions and compute flow accumulation.
      4. Return summary metrics for validation.

    Parameters
    ----------
    source_path : str
        Path to source DEM (float32/float64, projected CRS, nodata set).
    output_path : str
        Path for the conditioned, resampled output DEM.
    target_res_m : float
        Target grid spacing in metres.
    coarsen : bool
        True (default) for aggregation (average); False for refinement (cubic).

    Returns
    -------
    dict with resolution, cell count, pit stats, and p95 accumulation.
    """
    source = Path(source_path)
    if not source.exists():
        raise FileNotFoundError(f"Source DEM not found: {source}")

    # --- 1. Inspect ---
    with rasterio.open(source) as src:
        native_res = src.res[0]
        if not src.crs.is_projected:
            raise ValueError("Source DEM must use a projected CRS (e.g. UTM).")
        if src.nodata is None:
            logger.warning("Source DEM has no nodata value — fill artifacts may occur.")

        logger.info(
            "Source DEM: %.1f m native resolution, %d × %d cells, dtype=%s",
            native_res, src.height, src.width, src.dtypes[0],
        )

        resample_method = Resampling.average if coarsen else Resampling.cubic
        transform, width, height = calculate_default_transform(
            src.crs, src.crs, src.width, src.height, *src.bounds,
            resolution=target_res_m,
        )
        kwargs = src.meta.copy()
        kwargs.update({"transform": transform, "width": width, "height": height})

        # --- 2. Resample ---
        resampled_path = str(output_path).replace(".tif", "_raw_resample.tif")
        with rasterio.open(resampled_path, "w", **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=src.crs,
                    resampling=resample_method,
                )
        logger.info(
            "Resampled to %.1f m (%d × %d cells) via %s",
            target_res_m, height, width, resample_method.name,
        )

    # --- 3. Recondition ---
    dem = rd.LoadGDAL(resampled_path)
    filled = rd.FillDepressions(dem, in_place=False)

    fill_arr = np.array(filled) - np.array(dem)
    n_filled = int(np.sum(fill_arr > 1e-3))
    logger.info("Depression fill raised %d cells (%.2f%%)", n_filled, n_filled / dem.size * 100)

    rd.SaveGDAL(output_path, filled)
    logger.info("Conditioned DEM written: %s", output_path)

    # --- 4. Flow accumulation for validation handoff ---
    fa = rd.FlowAccumulation(filled, method="D8")
    fa_arr = np.array(fa, dtype=np.float64)
    fa_arr[fa_arr <= 0] = np.nan
    p95 = float(np.nanpercentile(fa_arr, 95))

    logger.info("Post-condition p95 accumulation: %.0f cells", p95)

    return {
        "source_res_m": native_res,
        "target_res_m": target_res_m,
        "output_path": output_path,
        "cell_count": height * width,
        "n_filled_cells": n_filled,
        "fa_p95": p95,
    }

import logging
import numpy as np

logger = logging.getLogger(__name__)

def validate_stream_overlap(modelled_stream: np.ndarray, reference_stream: np.ndarray) -> dict:
    """
    Compare modelled and reference binary stream rasters.
    Both arrays must be boolean (True = stream cell) at the same extent and resolution.
    """
    intersection = np.logical_and(modelled_stream, reference_stream).sum()
    union = np.logical_or(modelled_stream, reference_stream).sum()
    iou = float(intersection / union) if union > 0 else 0.0

    precision = float(intersection / modelled_stream.sum()) if modelled_stream.sum() > 0 else 0.0
    recall = float(intersection / reference_stream.sum()) if reference_stream.sum() > 0 else 0.0
    f1 = (2 * precision * recall / (precision + recall)) if (precision + recall) > 0 else 0.0

    logger.info("Stream overlap — IoU: %.3f, F1: %.3f, Precision: %.3f, Recall: %.3f",
                iou, f1, precision, recall)

    if f1 < 0.70:
        logger.warning("F1 below 0.70 — resolution may be too coarse for the target stream network.")

    return {"iou": iou, "f1": f1, "precision": precision, "recall": recall}