Watershed Modeling Python for hydrology automation

Converting NOAA Atlas 14 IDF Data into Model Forcing

The depth-duration-frequency table is the entry point to every design storm, and the step from the published values to a model input carries three quiet traps: the series type, the units, and the confidence limits that get discarded on the way. This guide covers the retrieval and conversion, as part of the precipitation forcing and design storms topic within rainfall-runoff modeling and hydrologic simulation.

Prerequisites

  • A catchment polygon in a projected CRS, so the centroid is area-weighted correctly.
  • requests and pandas.
  • A decision on the series type — partial duration or annual maximum — taken from the applicable design guidance rather than from the default.

Core Technique: Parse Once, Structure Properly

The point precipitation service returns a compact CSV whose rows are durations and whose columns are return periods, with the estimate and both confidence bounds in separate blocks. The useful representation for a model pipeline is the transpose: a tidy frame of duration, return period, estimate, lower and upper.

From Three Blocks to One Tidy Table The service response contains three matrices of duration by return period: the central estimate, the lower confidence bound and the upper bound. Reshaping them into one long table with duration, return period, estimate, lower and upper columns makes every downstream use a filter rather than a lookup. estimate block 19 × 9 matrix lower bound block same shape upper bound block same shape melt duration return period estimate lower upper 60 min 100 49.8 41.2 60.7 6 h 100 104.6 88.1 125.3 24 h 100 152.0 129.4 180.6 Keeping the bounds in the same row is what makes a three-run sensitivity analysis a one-line filter rather than a second retrieval.

Annotated Code Example

python 
import io
import logging

import geopandas as gpd
import pandas as pd
import requests

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

PFDS_URL = "https://hdsc.nws.noaa.gov/cgi-bin/hdsc/new/fe_text_mean.csv"
IN_TO_MM = 25.4

# Duration labels as the service returns them, mapped to minutes.
DURATION_MINUTES = {
    "5-min": 5, "10-min": 10, "15-min": 15, "30-min": 30, "60-min": 60,
    "2-hr": 120, "3-hr": 180, "6-hr": 360, "12-hr": 720, "24-hr": 1440,
    "2-day": 2880, "3-day": 4320, "4-day": 5760, "7-day": 10080,
}


def fetch_ddf(
    lat: float,
    lon: float,
    series: str = "pds",
    units: str = "english",
    timeout: int = 30,
) -> pd.DataFrame:
    """
    Retrieve depth-duration-frequency depths for one point, tidied to mm.

    Parameters
    ----------
    lat, lon : Point location in decimal degrees, WGS84.
    series   : "pds" for partial duration, "ams" for annual maximum.
    units    : Service units; the response is converted to mm regardless.

    Returns
    -------
    Tidy frame: duration_min, return_period_yr, estimate_mm, lower_mm, upper_mm.
    """
    params = {"lat": lat, "lon": lon, "type": "pf", "data": "depth",
              "units": units, "series": series}
    log.info("Requesting %s depths at (%.5f, %.5f)", series.upper(), lat, lon)
    resp = requests.get(PFDS_URL, params=params, timeout=timeout)
    resp.raise_for_status()

    text = resp.text
    if "PRECIPITATION FREQUENCY ESTIMATES" not in text.upper():
        raise ValueError("unexpected response — the point may be outside the "
                         "published coverage")

    # --- The response is three labelled blocks with the same shape. Split on
    # the block headers rather than on fixed line numbers, which change
    # between regional volumes. ---
    blocks = {}
    current, buf = None, []
    for line in text.splitlines():
        upper = line.upper()
        if "PRECIPITATION FREQUENCY ESTIMATES" in upper:
            current, buf = "estimate", []
        elif "UPPER BOUND" in upper:
            blocks[current] = buf
            current, buf = "upper", []
        elif "LOWER BOUND" in upper:
            blocks[current] = buf
            current, buf = "lower", []
        elif current and line.strip() and line[0].isdigit() or (
                current and line.strip().split(":")[0] in DURATION_MINUTES):
            buf.append(line)
    if current:
        blocks[current] = buf

    frames = []
    for name, lines in blocks.items():
        if not lines:
            continue
        df = pd.read_csv(io.StringIO("\n".join(lines)), header=None)
        df = df.rename(columns={0: "duration"})
        df["duration"] = df["duration"].astype(str).str.strip().str.rstrip(":")
        long = df.melt(id_vars="duration", var_name="rp_index", value_name=name)
        frames.append(long.set_index(["duration", "rp_index"]))

    merged = pd.concat(frames, axis=1).reset_index()

    # Return periods are in a fixed order in every regional volume.
    rp_order = [1, 2, 5, 10, 25, 50, 100, 200, 500, 1000]
    merged["return_period_yr"] = merged["rp_index"].map(
        {i + 1: rp for i, rp in enumerate(rp_order)}
    )
    merged["duration_min"] = merged["duration"].map(DURATION_MINUTES)
    merged = merged.dropna(subset=["duration_min", "return_period_yr"])

    # --- Convert to millimetres. The service returns inches under
    # units="english"; a pipeline that mixes the two produces depths off by
    # a factor of 25.4, which is large enough to notice and small enough
    # to be mistaken for a very wet region. ---
    scale = IN_TO_MM if units == "english" else 1.0
    for col in ("estimate", "lower", "upper"):
        if col in merged:
            merged[f"{col}_mm"] = merged[col].astype(float) * scale

    out = merged[["duration_min", "return_period_yr",
                  "estimate_mm", "lower_mm", "upper_mm"]].copy()
    out = out.sort_values(["return_period_yr", "duration_min"]).reset_index(drop=True)

    row = out[(out.return_period_yr == 100) & (out.duration_min == 1440)]
    if len(row):
        r = row.iloc[0]
        spread = 100.0 * (r.upper_mm - r.lower_mm) / r.estimate_mm
        log.info("100-yr 24-hr: %.1f mm [%.1f, %.1f] — a %.0f %% spread",
                 r.estimate_mm, r.lower_mm, r.upper_mm, spread)
    log.info("Retrieved %d duration/return-period pairs", len(out))
    return out


def basin_centroid(catchment_path: str) -> tuple[float, float]:
    """Area-weighted centroid of a catchment, returned as WGS84 lat/lon."""
    gdf = gpd.read_file(catchment_path)
    if gdf.crs is None or gdf.crs.is_geographic:
        raise ValueError("compute the centroid in a projected CRS — a centroid "
                         "in degrees is not the area-weighted centre")
    pt = gdf.geometry.union_all().centroid
    lonlat = gpd.GeoSeries([pt], crs=gdf.crs).to_crs(4326).iloc[0]
    log.info("Catchment centroid: %.5f, %.5f", lonlat.y, lonlat.x)
    return float(lonlat.y), float(lonlat.x)


# --- Example usage ---
# lat, lon = basin_centroid("catchment.gpkg")
# ddf = fetch_ddf(lat, lon, series="pds")
# design = ddf[ddf.return_period_yr == 100].set_index("duration_min")["estimate_mm"]
# print(design.round(1))

Parameter Reference

Parameter Options Effect
series pds / ams Partial duration is ~10 % higher at a 2-year return period; they converge above 10 years
units english / metric The service default is inches; failing to convert gives depths 25.4× too small
Centroid CRS Projected A centroid computed in degrees is not the area-weighted centre
Bounds Retrieved and kept A 100-year depth routinely carries a ±20 % published interval

Worked Example: Carrying the Uncertainty Through

Running the same model three times — at the lower bound, the estimate and the upper bound — costs almost nothing and changes the conversation about the result.

Forcing 24-hr depth Modelled peak Relative to central
Lower bound 129.4 mm 41.2 m³/s −18 %
Central estimate 152.0 mm 50.4 m³/s —
Upper bound 180.6 mm 62.9 m³/s +25 %
Forcing Uncertainty Against Parameter Uncertainty The published depth interval produces a design peak range from 41 to 63 cubic metres per second, a spread of 43 percent. The spread caused by plausible curve number and lag choices is 12 percent. The forcing interval dominates. 35 45 55 65 75 modelled design peak (m³/s) published depth interval — 41 to 63 forcing uncertainty 47 to 53 parameter uncertainty Calibrating parameters to three significant figures while ignoring a ±20 % forcing interval is effort spent on the smaller of the two uncertainties. same catchment, same model, 100-year 24-hour storm

Depth-duration curves are close to straight in log-log space and distinctly curved in linear space, so the interpolation choice matters most exactly where the published table is sparsest.

Interpolate the Depth-Duration Curve in Log Space Published depths at 60, 120 and 180 minutes. A linear interpolation between them cuts across the curve and underestimates the 90-minute depth by about 6 percent. Log-log interpolation follows the curve and lands within one percent. log-log: 46.2 mm linear: 43.4 mm — 6 % low 60 min 90 120 180 min duration — published points marked depth

Gotchas and Edge Cases

  • Units left in inches. A 152 mm storm becomes a 6 mm storm; the model runs and produces almost no runoff, which is often misdiagnosed as a curve number problem.
  • Annual maximum where partial duration is specified. Roughly 10 % low at a 2-year return period, converging to nothing above 10 years. State the series in the report.
  • Centroid in degrees. A geographic centroid is not the area-weighted centre and can sit several kilometres from it on an elongated basin.
  • One lookup for an orographic basin. Where depths vary strongly across the basin, sample several points or use the gridded product.
  • Bounds discarded on retrieval. Recovering them later means a second request; keeping them costs two columns.
  • Interpolating between published durations linearly. Depth-duration curves are close to linear in log-log space and distinctly not in linear space. Interpolate in logs.
  • Assuming coverage. The published volumes do not cover every jurisdiction, and a request outside coverage returns a page rather than data. Check the response before parsing.