Watershed Modeling Python for hydrology automation

SRTM and LiDAR Data Acquisition

Reliable elevation data is the geometric backbone of every hydrologic simulation, from overland flow routing to floodplain delineation. As part of the Hydrology Data Preparation & DEM Processing pipeline, this stage focuses on replacing manual tile selection and GUI-based downloads with reproducible, API-driven acquisition scripts. Automated retrieval eliminates version drift, ensures spatial consistency, and establishes a clean handoff to DEM Pit Filling Algorithms and coordinate reference system alignment without reprocessing raw data by hand.

The two primary sources covered here serve fundamentally different scales and precision requirements: the NASA Shuttle Radar Topography Mission (SRTM) delivers near-global 30 m raster coverage, while airborne LiDAR point clouds provide sub-meter vertical fidelity for channel morphology, engineered structures, and urban drainage geometry.

Prerequisites & Environment Setup

Python stack: Python 3.9+, isolated in a venv or conda environment.

bash 
pip install earthaccess>=0.8 rasterio>=1.3 geopandas>=0.12 \
            shapely>=2.0 numpy>=1.23 laspy>=2.4 requests>=2.31 tqdm>=4.65

Input data requirements:

Requirement SRTM LiDAR
Format GeoTIFF (.tif) LAS / LAZ point cloud
Native CRS WGS84 EPSG:4326 Varies — typically UTM or state-plane
Vertical datum EGM96 geoid NAVD88 (US) / national geoid
Resolution 1 arc-second (~30 m) 0.25 m – 2 m after rasterization
Nodata encoding −32768 (Int16) N/A (sparse points)

Authentication:

  • NASA Earthdata: free account at urs.earthdata.nasa.gov. Store credentials in ~/.netrc or export EARTHDATA_USERNAME / EARTHDATA_PASSWORD as environment variables — never hardcode tokens in version-controlled scripts.
  • OpenTopography: register at opentopography.org/developers for an API key. Export as OPENTOP_API_KEY.

Storage: budget at least 50 GB SSD for staging raw tiles, intermediate mosaics, and point-cloud LAZ files. SRTM tiles covering a large river basin can reach 10–15 GB uncompressed.

Algorithm Mechanics: How the Acquisition Pipeline Works

Both sources follow the same logical sequence: spatial query → metadata resolution → authenticated retrieval → structural validation. Understanding where each stage can fail prevents silent data corruption downstream.

SRTM and LiDAR Data Acquisition Pipeline Five sequential processing stages — Spatial Query, Metadata Resolve, Authenticated Retrieval, Mosaic and Reproject, Validate — connected by arrows. A dashed feedback arrow from Validate back to Authenticated Retrieval represents the retry path for corrupted or missing tiles. Spatial Query bbox · catalog API Metadata Resolve tile IDs · versions Authenticated Retrieval earthaccess · laspy Mosaic & Reproject rasterio merge · UTM Validate & Stage extent · nodata · CRS retry on corrupt / missing tile

Spatial query converts a watershed boundary polygon into a bounding-box filter compatible with both the NASA CMR API and the OpenTopography catalog. Metadata resolution maps bounding boxes to tile identifiers, dataset versions, and access URLs. Authenticated retrieval streams files directly to local staging directories. Mosaic and reproject stitches multi-tile rasters into a single seamless grid in the project CRS. Validation checks extent, nodata definition, CRS encoding, and — for point clouds — point density before any conditioning occurs.

Dataset Selection: SRTM vs. LiDAR

Choosing between SRTM and LiDAR before writing acquisition code determines every downstream resolution and accuracy decision.

Criterion SRTM GL1 (30 m) Airborne LiDAR (≤1 m)
Geographic coverage Near-global (60°N–56°S) Regional; US coverage via 3DEP
Vertical accuracy (RMSE) ~5–9 m 0.05–0.15 m
Canopy penetration Surface model (DSM) Bare-earth returns available
Typical file size per tile 25 MB 500 MB – 10 GB
Best for Basin-scale routing, continental analysis Urban drainage, channel bathymetry, floodplain mapping
Temporal baseline February 2000 Varies by survey year

Coarser DEMs smooth channel networks and underestimate flow accumulation in incised valleys. Ultra-high-resolution LiDAR introduces significant computational overhead and requires aggressive classification filtering to separate ground returns from vegetation and building artifacts. Align the dataset choice with your hydraulic solver’s cell-size requirements and available compute before initiating downloads.

Anticipate downstream coordinate transformations early. Raw SRTM tiles arrive in WGS84 (EPSG:4326). LiDAR datasets frequently use UTM or local state-plane projections. Establishing a target CRS at this stage prevents geometric distortion during rasterization — a topic covered in depth in coordinate reference system alignment.

Step-by-Step Workflow

Step 1 — Define Spatial Boundary and Query Catalog

Load your watershed boundary as a GeoDataFrame, reproject to EPSG:4326 for catalog compatibility, and extract bounds.

python 
import logging
import geopandas as gpd

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

def get_wgs84_bounds(shapefile: str) -> tuple[float, float, float, float]:
    """Reproject watershed boundary to WGS84 and return (minx, miny, maxx, maxy)."""
    watershed = gpd.read_file(shapefile)
    log.info("Loaded watershed: %d feature(s), CRS=%s", len(watershed), watershed.crs)
    watershed_wgs84 = watershed.to_crs(epsg=4326)
    bbox = tuple(watershed_wgs84.total_bounds)   # (minx, miny, maxx, maxy)
    log.info("Catalog query bbox: %.4f, %.4f, %.4f, %.4f", *bbox)
    return bbox

bbox = get_wgs84_bounds("data/watersheds/basin_04.shp")

For granular control over tile selection, chunked requests, and metadata caching, see how to download SRTM DEMs for Python hydrology workflows.

Step 2 — Authenticated SRTM Tile Retrieval

NASA distributes SRTM data through Earthdata Cloud. The earthaccess library negotiates OAuth tokens, paginates CMR search results, and supports resumable downloads transparently.

python 
import earthaccess
from pathlib import Path

def fetch_srtm_tiles(bbox: tuple, output_dir: str = "data/raw/srtm") -> list[Path]:
    """Download SRTM GL1 (30 m) tiles covering the bounding box."""
    earthaccess.login()   # reads ~/.netrc or EARTHDATA_* env vars
    log.info("Searching SRTM GL1 catalog for bbox %s", bbox)

    results = earthaccess.search_data(
        short_name="SRTMGL1",
        bounding_box=bbox,
        temporal=("2000-02-01", "2000-02-23")  # SRTM acquisition window
    )
    log.info("Found %d SRTM granule(s)", len(results))

    out = Path(output_dir)
    out.mkdir(parents=True, exist_ok=True)

    files = earthaccess.download(results, str(out))
    log.info("Retrieved %d file(s) to %s", len(files), out)
    return [Path(f) for f in files]

srtm_files = fetch_srtm_tiles(bbox)

earthaccess.download is idempotent: files already present and matching the remote checksum are skipped. For large basins spanning dozens of tiles, wrap the download call in a retry decorator with exponential backoff to respect CMR rate limits.

Step 3 — Regional LiDAR Point-Cloud Fetching

LiDAR acquisition uses the OpenTopography REST catalog. The API returns dataset metadata (survey footprint, acquisition year, point density, file URL) as JSON, allowing programmatic selection before any bytes are transferred.

python 
import os
import requests
from tqdm import tqdm

OPENTOP_CATALOG = "https://portal.opentopography.org/API/otCatalog"

def list_lidar_datasets(bbox: tuple) -> list[dict]:
    """Return OpenTopography LiDAR datasets intersecting the bounding box."""
    minx, miny, maxx, maxy = bbox
    params = {
        "productFormat": "PointCloud",
        "minx": minx, "miny": miny, "maxx": maxx, "maxy": maxy,
        "detail": "true",
        "outputFormat": "json",
    }
    api_key = os.getenv("OPENTOP_API_KEY")
    if api_key:
        params["API_Key"] = api_key

    resp = requests.get(OPENTOP_CATALOG, params=params, timeout=30)
    resp.raise_for_status()
    datasets = resp.json().get("Datasets", [])
    log.info("OpenTopography returned %d LiDAR dataset(s)", len(datasets))
    return datasets


def download_lidar(dataset: dict, output_dir: str = "data/raw/lidar") -> Path:
    """Stream-download a single LiDAR LAZ file with progress tracking."""
    out = Path(output_dir)
    out.mkdir(parents=True, exist_ok=True)

    url = dataset.get("downloadURL")
    dataset_id = dataset.get("datasetID", "unknown")
    dest = out / f"{dataset_id}.laz"

    if dest.exists():
        log.info("LAZ already present, skipping: %s", dest)
        return dest

    log.info("Downloading LiDAR dataset %s → %s", dataset_id, dest)
    with requests.get(url, stream=True, timeout=120) as r:
        r.raise_for_status()
        total = int(r.headers.get("content-length", 0))
        with open(dest, "wb") as fh, tqdm(total=total, unit="B", unit_scale=True,
                                           desc=dataset_id) as bar:
            for chunk in r.iter_content(chunk_size=65536):
                fh.write(chunk)
                bar.update(len(chunk))

    log.info("Saved %s (%.1f MB)", dest, dest.stat().st_size / 1e6)
    return dest


datasets = list_lidar_datasets(bbox)
lidar_files = [download_lidar(ds) for ds in datasets[:3]]  # limit for demo

Always verify dataset licensing before commercial deployment — some regional LiDAR surveys carry non-commercial restrictions embedded in the metadata JSON.

Step 4 — Mosaic Assembly and Reprojection

Individual SRTM tiles share the WGS84 geographic coordinate system but must be merged and reprojected to a metric CRS before flow routing. Use rasterio.merge followed by a rasterio.warp.reproject to produce a single Cloud-Optimized GeoTIFF.

python 
import rasterio
from rasterio.merge import merge
from rasterio.warp import calculate_default_transform, reproject, Resampling

def mosaic_and_reproject(tile_paths: list[Path], target_epsg: int,
                          output_path: str = "data/processed/srtm_mosaic.tif") -> Path:
    """Merge SRTM tiles and reproject to target CRS, writing a COG."""
    src_files = [rasterio.open(p) for p in tile_paths]
    log.info("Merging %d SRTM tile(s)", len(src_files))
    mosaic, transform = merge(src_files)

    dst_crs = f"EPSG:{target_epsg}"
    dst_transform, width, height = calculate_default_transform(
        src_files[0].crs, dst_crs,
        mosaic.shape[2], mosaic.shape[1],
        left=src_files[0].bounds.left, bottom=src_files[0].bounds.bottom,
        right=src_files[-1].bounds.right, top=src_files[0].bounds.top
    )

    out = Path(output_path)
    out.parent.mkdir(parents=True, exist_ok=True)
    profile = src_files[0].profile.copy()
    profile.update({
        "crs": dst_crs, "transform": dst_transform,
        "width": width, "height": height,
        "driver": "GTiff",
        "compress": "deflate", "tiled": True,
        "blockxsize": 512, "blockysize": 512,
        "interleave": "band"
    })

    with rasterio.open(out, "w", **profile) as dst:
        reproject(
            source=mosaic[0], destination=rasterio.band(dst, 1),
            src_transform=transform, src_crs=src_files[0].crs,
            dst_transform=dst_transform, dst_crs=dst_crs,
            resampling=Resampling.bilinear
        )

    for src in src_files:
        src.close()

    log.info("Mosaic written: %s (%d×%d px, EPSG:%d)", out, width, height, target_epsg)
    return out

mosaic_path = mosaic_and_reproject(srtm_files, target_epsg=32618)

The bilinear resampling kernel preserves smooth elevation gradients during reprojection. For spatial resolution tradeoffs analysis — such as comparing SRTM at 30 m against a LiDAR-derived 1 m DEM — use nearest-neighbor to avoid interpolating integer class rasters or categorical masks.

Step 5 — Structural Validation and Staging

Raw elevation data frequently contains voids, inconsistent nodata values, or CRS metadata that contradicts the file contents. Validate before passing to any conditioning step.

python 
import laspy

def validate_srtm_tile(filepath: Path) -> bool:
    """Check band count, nodata definition, and CRS encoding."""
    try:
        with rasterio.open(filepath) as src:
            if src.count != 1:
                log.warning("Unexpected band count %d in %s", src.count, filepath)
                return False
            if src.nodata is None:
                log.warning("No nodata value defined in %s", filepath)
                return False
            if src.crs is None:
                log.warning("No CRS defined in %s", filepath)
                return False
            log.debug("SRTM tile valid: %s (%.4f°×%.4f°)", filepath,
                      src.res[1], src.res[0])
            return True
    except Exception as exc:
        log.error("Could not open %s: %s", filepath, exc)
        return False


def validate_lidar(filepath: Path, min_point_count: int = 1000) -> bool:
    """Confirm the LAZ file opens and contains a minimum point density."""
    try:
        with laspy.open(filepath) as las:
            count = las.header.point_count
            if count < min_point_count:
                log.warning("LiDAR %s has only %d points (threshold=%d)",
                            filepath, count, min_point_count)
                return False
            log.debug("LiDAR valid: %s (%d points)", filepath, count)
            return True
    except Exception as exc:
        log.error("Could not open LiDAR %s: %s", filepath, exc)
        return False


srtm_valid = [f for f in Path("data/raw/srtm").glob("*.tif") if validate_srtm_tile(f)]
lidar_valid = [f for f in Path("data/raw/lidar").glob("*.laz") if validate_lidar(f)]
log.info("Validated %d SRTM tile(s) and %d LiDAR file(s)",
         len(srtm_valid), len(lidar_valid))

Validation failures should trigger automated alerts rather than silent skips. Missing tiles or corrupted point clouds propagate errors into flow direction matrices and hydraulic roughness calculations.

Production-Ready Code

The function below integrates all steps into a single, idempotent pipeline entry point suitable for CI/CD or Prefect/Airflow task wrappers.

python 
import logging
from pathlib import Path
import earthaccess
import geopandas as gpd
import rasterio
from rasterio.merge import merge
from rasterio.warp import calculate_default_transform, reproject, Resampling
import laspy
import requests
import os
from tqdm import tqdm

log = logging.getLogger(__name__)


def acquire_elevation_data(
    watershed_shp: str,
    target_epsg: int,
    srtm_dir: str = "data/raw/srtm",
    lidar_dir: str = "data/raw/lidar",
    mosaic_path: str = "data/processed/dem_mosaic.tif",
    max_lidar_datasets: int = 5,
) -> dict[str, Path | list[Path]]:
    """
    End-to-end SRTM and LiDAR acquisition pipeline.

    Parameters
    ----------
    watershed_shp : path to the watershed boundary shapefile
    target_epsg   : EPSG code for the project CRS (e.g. 32618 for UTM 18N)
    srtm_dir      : staging directory for raw SRTM tiles
    lidar_dir     : staging directory for LAZ point clouds
    mosaic_path   : output path for the reprojected SRTM mosaic COG
    max_lidar_datasets : maximum number of LiDAR datasets to download

    Returns
    -------
    dict with keys 'mosaic', 'lidar_files', 'srtm_valid', 'lidar_valid'
    """
    # -- 1. Spatial boundary
    watershed = gpd.read_file(watershed_shp)
    bbox = tuple(watershed.to_crs(epsg=4326).total_bounds)
    log.info("Acquisition bbox (WGS84): %s", bbox)

    # -- 2. SRTM retrieval
    earthaccess.login()
    results = earthaccess.search_data(
        short_name="SRTMGL1",
        bounding_box=bbox,
        temporal=("2000-02-01", "2000-02-23"),
    )
    log.info("SRTM: %d granule(s) found", len(results))
    Path(srtm_dir).mkdir(parents=True, exist_ok=True)
    raw_srtm = [Path(f) for f in earthaccess.download(results, srtm_dir)]

    # -- 3. LiDAR retrieval
    minx, miny, maxx, maxy = bbox
    params = {
        "productFormat": "PointCloud",
        "minx": minx, "miny": miny, "maxx": maxx, "maxy": maxy,
        "detail": "true", "outputFormat": "json",
    }
    api_key = os.getenv("OPENTOP_API_KEY")
    if api_key:
        params["API_Key"] = api_key
    resp = requests.get(
        "https://portal.opentopography.org/API/otCatalog",
        params=params, timeout=30
    )
    resp.raise_for_status()
    datasets = resp.json().get("Datasets", [])[:max_lidar_datasets]
    log.info("LiDAR: %d dataset(s) selected", len(datasets))

    Path(lidar_dir).mkdir(parents=True, exist_ok=True)
    lidar_files: list[Path] = []
    for ds in datasets:
        url = ds.get("downloadURL")
        dsid = ds.get("datasetID", "unknown")
        dest = Path(lidar_dir) / f"{dsid}.laz"
        if dest.exists():
            log.info("Skipping existing LAZ: %s", dest)
            lidar_files.append(dest)
            continue
        with requests.get(url, stream=True, timeout=120) as r:
            r.raise_for_status()
            total = int(r.headers.get("content-length", 0))
            with open(dest, "wb") as fh, tqdm(total=total, unit="B",
                                               unit_scale=True, desc=dsid) as bar:
                for chunk in r.iter_content(chunk_size=65536):
                    fh.write(chunk)
                    bar.update(len(chunk))
        log.info("Downloaded %s (%.1f MB)", dest, dest.stat().st_size / 1e6)
        lidar_files.append(dest)

    # -- 4. Mosaic and reproject SRTM
    srcs = [rasterio.open(p) for p in raw_srtm]
    mosaic_arr, mosaic_tf = merge(srcs)
    dst_crs = f"EPSG:{target_epsg}"
    dst_tf, width, height = calculate_default_transform(
        srcs[0].crs, dst_crs,
        mosaic_arr.shape[2], mosaic_arr.shape[1],
        left=srcs[0].bounds.left, bottom=srcs[0].bounds.bottom,
        right=srcs[-1].bounds.right, top=srcs[0].bounds.top,
    )
    profile = srcs[0].profile.copy()
    profile.update({
        "crs": dst_crs, "transform": dst_tf,
        "width": width, "height": height,
        "driver": "GTiff", "compress": "deflate",
        "tiled": True, "blockxsize": 512, "blockysize": 512,
    })
    mosaic_out = Path(mosaic_path)
    mosaic_out.parent.mkdir(parents=True, exist_ok=True)
    with rasterio.open(mosaic_out, "w", **profile) as dst:
        reproject(
            source=mosaic_arr[0], destination=rasterio.band(dst, 1),
            src_transform=mosaic_tf, src_crs=srcs[0].crs,
            dst_transform=dst_tf, dst_crs=dst_crs,
            resampling=Resampling.bilinear,
        )
    for src in srcs:
        src.close()
    log.info("SRTM mosaic: %s (%d×%d px)", mosaic_out, width, height)

    # -- 5. Validate
    srtm_valid = [p for p in raw_srtm if _validate_raster(p)]
    lidar_valid = [p for p in lidar_files if _validate_lidar(p)]
    log.info("Validation: %d/%d SRTM tiles OK, %d/%d LiDAR files OK",
             len(srtm_valid), len(raw_srtm), len(lidar_valid), len(lidar_files))

    return {
        "mosaic": mosaic_out,
        "lidar_files": lidar_files,
        "srtm_valid": srtm_valid,
        "lidar_valid": lidar_valid,
    }


def _validate_raster(filepath: Path) -> bool:
    try:
        with rasterio.open(filepath) as src:
            return src.count == 1 and src.nodata is not None and src.crs is not None
    except Exception as exc:
        log.error("Raster validation failed for %s: %s", filepath, exc)
        return False


def _validate_lidar(filepath: Path, min_points: int = 1000) -> bool:
    try:
        with laspy.open(filepath) as las:
            return las.header.point_count >= min_points
    except Exception as exc:
        log.error("LiDAR validation failed for %s: %s", filepath, exc)
        return False

Validation Protocol

After acquisition, verify the mosaic and point clouds against independent reference data before any conditioning step.

SRTM mosaic checks:

  1. Open the mosaic in QGIS or with rasterio.open and confirm the spatial extent matches the watershed bounding box — tile seams sometimes leave 1–2 pixel gaps at antimeridian crossings.
  2. Compute a histogram of elevation values (numpy.histogram) and flag any anomalous spikes near the nodata value (−32768) that indicate void clusters.
  3. Compare a 1 km² test patch against known survey benchmarks or Copernicus DEM GLO-30 to catch systematic vertical biases above 10 m RMSE.
  4. Confirm the output CRS with rio info --crs data/processed/dem_mosaic.tif — a missing authority code causes silent failures in downstream rasterio operations.

LiDAR point-cloud checks:

  1. Inspect point density: the LAZ header’s point_count divided by the survey footprint area should match the advertised density (typically 2–20 pts/m²). Values below 0.5 pts/m² suggest the download is incomplete or the wrong classification filter was applied.
  2. Verify the vertical datum entry in the VLR (Variable Length Records). For US datasets, laspy exposes this via las.header.vlrs — look for a WKT or PROJ4 string containing NAVD88 or NGVD29 to confirm vertical reference before merging with SRTM.
  3. Run a ground-return count: (las.classification == 2).sum() / las.header.point_count should exceed 0.1 (10%) for adequately filtered bare-earth data. Ratios below 0.05 indicate heavy canopy or incorrect classification.

Common Failure Modes & Optimization

Authentication token expiry mid-download. earthaccess tokens expire after one hour. For large tile sets, wrap the download loop in a function that calls earthaccess.login() once before each batch of 20 granules rather than once globally.

SRTM void tiles in mountainous terrain. The Shuttle Radar had radar shadow in steep valleys — affected tiles contain contiguous nodata patches that masquerade as valid pixels if nodata is −32768 and the valley floor genuinely reaches that elevation. Detect by counting pixels equal to src.nodata as a fraction of total area; patches exceeding 5% warrant void-filling with an alternative source before DEM Pit Filling Algorithms are applied.

CRS mismatch between LiDAR tiles from different survey years. Successive LiDAR campaigns over the same region sometimes use different horizontal datums (NAD83 vs NAD83(2011)). The difference is small (up to 1 m) but sufficient to create step discontinuities in merged DEMs. Check the WKT for NAD83(2011) identifiers and apply a datum-shift grid before rasterization.

Memory exhaustion during mosaic of high-resolution LiDAR. A 1 m LiDAR DEM covering 500 km² loads as a Float32 array of roughly 500 million pixels — 2 GB minimum. Use rasterio.merge with method="first" and chunked windowed writes, or switch to the spatial resolution tradeoffs strategy of resampling to 5 m before mosaic assembly.

Tile boundary seams after reprojection. Bilinear resampling at tile edges can introduce subtle intensity gradients if source tiles have inconsistent gain calibration. Write each reprojected tile to a virtual raster (.vrt) first and inspect seams at 1:1 zoom before finalising the mosaic.

OpenTopography API rate limits. Unauthenticated requests are capped at 500 per day. Add API_Key to every request and cache catalog JSON responses locally with a 24-hour TTL to avoid re-fetching metadata on repeated pipeline runs.

When to Use This vs. Alternatives

SRTM acquisition via earthaccess is the right choice for regional to continental basin analysis where 30 m resolution is sufficient and global coverage is required. Switch to the USGS 3DEP 1 m LiDAR workflow when the study area is within the contiguous US, Alaska, or Hawaii, and the hydraulic model requires sub-grid channel representation. For coastal and island terrain where neither SRTM nor domestic LiDAR is available, the Copernicus DEM GLO-30 (also 30 m, but with improved void filling in mountainous zones) is a drop-in replacement that uses the same rasterio-based mosaic workflow.

Once acquired, all three elevation sources follow the same conditioning path: apply DEM Pit Filling Algorithms to enforce hydrologic connectivity, then validate the flow network against NHD or comparable authoritative hydrography before moving to flow direction computation.