WPC QPF DSS-to-HEC-RAS Rain-on-Grid Qualification

This notebook downloads a current Weather Prediction Center (WPC) six-hour quantitative precipitation forecast, translates the published GRIB2 grids to HEC-DSS with qpkit, and verifies the complete HEC-RAS forcing chain on a real RasExamples model: source GRIB, DSS read-back, temporary plan HDF, completed plan HDF, and hydraulic response.

WPC's published GRIB2 compression is not accepted by native HEC-RAS import through 7.0.1. The DSS translation demonstrated here is therefore mandatory, not an optional alternate route. Notebook 916 provides the analogous archived HRRR workflow.

Qualification scope and reproducibility

QPFRequest 0.1.0 exposes the current UTC day's WPC directory but no historical date argument. This is consequently an operational live-data example: it probes eligible cycles newest-to-oldest, rejects incomplete cycles, downloads the full f006 through f168 sequence, and then selects the wettest consecutive 24-hour window over the model. The committed outputs record the exact cycle used for this qualification run.

qpkit==0.1.0 uses a pydsstools wheel that requires NumPy 1.x. Keep this workflow in a dedicated environment rather than downgrading a general ras-commander environment:

Bash Session

uv venv working/.venv-wpc-qpf
uv pip install --python working/.venv-wpc-qpf/Scripts/python.exe "numpy<2" "scipy<1.17" -e ".[full,dss]" nbclient ipykernel
uv pip install --python working/.venv-wpc-qpf/Scripts/python.exe "git+https://github.com/gyanz/[email protected]"

The model compute uses HEC-RAS 7.0 and the Bald Eagle Creek RasExamples fixture. Use RAS_COMMANDER_EXAMPLE_RUN_ROOT to provide a short local working path when HEC-Vortex/DSS or Wine path handling requires one. Generated model and forecast files remain outside the committed example.

Manual review remains required: inspect the WPC product/cycle, DSS catalog and maps, HEC-RAS runtime and BCO messages, precipitation coverage/timing/units, terrain coverage, mass balance, convergence, and hydraulic plausibility.

Python

# Development mode: import repository source when run from a clone.
USE_LOCAL_SOURCE = True

import sys
from pathlib import Path

repo_root = Path.cwd().parent if Path.cwd().name.lower() == "examples" else Path.cwd()
if USE_LOCAL_SOURCE and str(repo_root) not in sys.path:
    sys.path.insert(0, str(repo_root))

import ras_commander

print(f"ras-commander {ras_commander.__version__} ({'repository source' if USE_LOCAL_SOURCE else 'installed package'})")

Text Only

ras-commander 0.102.0 (repository source)

Python

from contextlib import redirect_stderr, redirect_stdout
from datetime import datetime, timedelta, timezone
from importlib.metadata import version as package_version
from io import StringIO
import logging
import os
from pathlib import Path
import re
import warnings

import cfgrib
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import rasterio
from rasterio.windows import Window, from_bounds
from rasterio.warp import transform_bounds
from qpkit import BoundingBox, QPFGridOptions, QPFRequest, QPKit

from ras_commander import (
    HdfBndry,
    HdfInfiltration,
    HdfMesh,
    HdfProject,
    HdfResultsMesh,
    HdfResultsPlan,
    HdfStruc,
    RasCmdr,
    RasDss,
    RasExamples,
    RasPlan,
    RasPreprocess,
    RasUnsteady,
    get_ras_exe,
    init_ras_project,
)

ipython = get_ipython()
if ipython is not None:
    ipython.run_line_magic("matplotlib", "inline")

for logger_name in (
    "httpx", "pydsstools", "pydsstools.heclib.logging", "qpkit",
    "rasterio", "ras_commander",
):
    logging.getLogger(logger_name).setLevel(logging.WARNING)

QPKIT_VERSION = package_version("qpkit")
if QPKIT_VERSION != "0.1.0":
    raise RuntimeError(f"This qualified workflow requires qpkit 0.1.0; found {QPKIT_VERSION}")

run_root = Path(os.environ.get(
    "RAS_COMMANDER_EXAMPLE_RUN_ROOT",
    repo_root / "working" / "example_926_wpc_qpf_qualification",
)).resolve()
run_root.mkdir(parents=True, exist_ok=True)

def rel(path: Path) -> str:
    path = Path(path).resolve()
    try:
        return path.relative_to(repo_root.resolve()).as_posix()
    except ValueError:
        return path.name

RAS_VERSION = "7.0"
WINDOW_RECORDS = 4
QPF_INTERVAL_HOURS = 6
EXPECTED_QPF_FORECAST_HOURS = tuple(range(6, 169, QPF_INTERVAL_HOURS))
WPC_CYCLE_HOURS = (0, 6, 12, 18)
PRECIPITATION_BUFFER_PERCENT = 50.0
COMPUTATION_INTERVAL = "5SEC"
MAX_SOURCE_DSS_MEAN_DIFFERENCE_MM = 0.30
MAX_WSEL_ERROR_FT = 1.00
MAX_CELLS_WSEL_ERROR_GT_0_1_FT = 5
MESH_NAME = "BaldEagleCr"
TEMPLATE_PLAN = "06"
TEMPLATE_UNSTEADY = "03"
TEMPLATE_GEOMETRY = "09"
ras_exe = get_ras_exe(RAS_VERSION)

display(pd.DataFrame([{
    "product": "WPC 6-hour QPF",
    "forecast horizon": "f006-f168",
    "selected model window": f"{WINDOW_RECORDS * QPF_INTERVAL_HOURS} hours",
    "HEC-RAS": RAS_VERSION,
    "qpkit": QPKIT_VERSION,
}]))

Text Only

2026-09-26 04:43:19 - pydsstools.heclib.logging - INFO - DssLogger(GLOBAL, level=TERSE): level set — method=GLOBAL, level=TERSE.

	product	forecast horizon	selected model window	HEC-RAS	qpkit
0	WPC 6-hour QPF	f006-f168	24 hours	7.0	0.1.0

Extract and inspect the real model fixture

The precipitation crop is derived from the geometry HDF rather than a hard-coded region. The retained inflow and the no-rain comparison isolate the incremental hydraulic effect of the WPC forcing.

Python

project_path = RasExamples.extract_project(
    "BaldEagleCrkMulti2D",
    output_path=run_root / "model",
    suffix="wpc_qpf_qualification",
)
ras = init_ras_project(
    project_path,
    ras_version=RAS_VERSION,
    load_results_summary=False,
    hide_intro=True,
)

template_plan_row = ras.plan_df.loc[
    ras.plan_df["plan_number"].astype(str).str.zfill(2) == TEMPLATE_PLAN
].iloc[0]
assert str(template_plan_row["unsteady_number"]).zfill(2) == TEMPLATE_UNSTEADY
assert str(template_plan_row["geometry_number"]).zfill(2) == TEMPLATE_GEOMETRY

geometry_hdf = ras.project_folder / f"{ras.project_name}.g{TEMPLATE_GEOMETRY}.hdf"
template_unsteady_path = ras.project_folder / f"{ras.project_name}.u{TEMPLATE_UNSTEADY}"
project_bounds = tuple(float(v) for v in HdfProject.get_project_bounds_latlon(
    geometry_hdf, buffer_percent=0.0
))
BOUNDS = tuple(float(v) for v in HdfProject.get_project_bounds_latlon(
    geometry_hdf, buffer_percent=PRECIPITATION_BUFFER_PERCENT
))
WPC_BBOX = BoundingBox(
    left_lon=BOUNDS[0], bottom_lat=BOUNDS[1],
    right_lon=BOUNDS[2], top_lat=BOUNDS[3],
)
mesh_areas = HdfMesh.get_mesh_areas(geometry_hdf)
mesh_attributes = HdfMesh.get_mesh_area_attributes(geometry_hdf)
mesh_cell_count = int(np.asarray(mesh_attributes.loc["Cell Count", "Value"]).reshape(-1)[0])
mesh_cell_points = HdfMesh.get_mesh_cell_points(geometry_hdf)
bc_lines = HdfBndry.get_bc_lines(geometry_hdf)
structures = HdfStruc.get_structures(geometry_hdf)
infiltration_stats = HdfInfiltration.get_preprocessed_infiltration_stats(
    geometry_hdf, mesh_name=MESH_NAME, variable="Curve Number"
)
flow_rows = RasUnsteady.get_inline_hydrograph_boundaries(
    template_unsteady_path, ras_object=ras
)
flow_rows = flow_rows.loc[flow_rows["bc_type"] == "Flow Hydrograph"]
assert not flow_rows.empty
flow_values = np.asarray(flow_rows.iloc[0]["values"], dtype=float)

display(pd.DataFrame([{
    "model": ras.project_name,
    "mesh": MESH_NAME,
    "computational cells": mesh_cell_count,
    "project WGS84 bounds": project_bounds,
    "buffered WPC bounds": BOUNDS,
    "terrain": Path(ras.rasmap_df.iloc[0]["terrain_hdf_path"][0]).name,
    "infiltration": Path(ras.rasmap_df.iloc[0]["infiltration_hdf_path"][0]).name,
}]).T.rename(columns={0: "value"}))
display(infiltration_stats.round(3))

fig, (ax_model, ax_flow) = plt.subplots(1, 2, figsize=(15, 5.8), constrained_layout=True)
mesh_areas.boundary.plot(ax=ax_model, color="#263238", linewidth=1.2, label="2D flow area")
if not bc_lines.empty:
    bc_lines.plot(ax=ax_model, color="#d95f02", linewidth=2.0, label="Boundary-condition lines")
if not structures.empty:
    structures.plot(ax=ax_model, color="#6a3d9a", linewidth=2.0, label="SA/2D structure")
ax_model.set(title="Bald Eagle Creek Model", xlabel="Easting (ft)", ylabel="Northing (ft)")
ax_model.set_aspect("equal")
ax_model.legend(loc="best")

ax_flow.plot(np.arange(len(flow_values)), flow_values, color="#2b6f9f", linewidth=2)
ax_flow.set(title="Retained Synthetic Inflow", xlabel="Hours from simulation start", ylabel="Flow (cfs)")
ax_flow.grid(alpha=0.25)
plt.show()

	value
model	BaldEagleDamBrk
mesh	BaldEagleCr
computational cells	18066
project WGS84 bounds	(-77.75914883720999, 40.96067310349849, -77.32...
buffered WPC bounds	(-77.86715797330825, 40.9043529357912, -77.219...
terrain	Terrain50.hdf
infiltration	Infiltration.hdf

	mesh_name	n_cells	min_val	max_val	mean_val	median_val	std_val	variable
0	BaldEagleCr	18066	36.0	100.0	73.645	73.0	14.803	Curve Number

Select a complete live WPC cycle and translate to DSS

The newest published cycle can be incomplete while WPC is uploading it. A candidate is accepted only when its forecast-hour set is exactly f006 through f168. The source GRIB files are downloaded at their published CONUS extent.

qpkit 0.1.0's QPFGridOptions.extents mask does not preserve the requested WGS84 crop reliably for this product. The workflow therefore performs an explicit native-projection crop with rasterio, verifies every cropped GRIB, and then calls the public kit.qpf_dss.write() route without the faulty extent mask. This avoids both a huge CONUS DSS and a silently undersized raster.

Python

qpf_grib_dir = run_root / "grib"
cropped_grib_dir = run_root / "cropped_grib"
precipitation_dir = ras.project_folder / "Precipitation" / "wpc_qpf"
qpf_grib_dir.mkdir(parents=True, exist_ok=True)
cropped_grib_dir.mkdir(parents=True, exist_ok=True)
precipitation_dir.mkdir(parents=True, exist_ok=True)
qpf_dss_file = precipitation_dir / "wpc_qpf_bald_eagle.dss"
if qpf_dss_file.exists():
    qpf_dss_file.unlink()

def qpf_forecast_hour(item) -> int:
    filename = item.filename if hasattr(item, "filename") else Path(item).name
    match = re.search(r"f(\d{3})\.grb$", Path(filename).name)
    if match is None:
        raise RuntimeError(f"Could not parse WPC QPF forecast hour: {filename}")
    return int(match.group(1))

current_utc = datetime.now(timezone.utc)
candidate_cycles = tuple(c for c in reversed(WPC_CYCLE_HOURS) if c <= current_utc.hour)
rejected_cycles = []
with QPKit(logger=logging.getLogger("qpkit")) as kit:
    qpf_request = None
    for cycle_hour in candidate_cycles:
        candidate = QPFRequest(interval=QPF_INTERVAL_HOURS, cycle_hour=cycle_hour)
        candidate_plan = kit.build_plan_qpf(candidate, qpf_grib_dir)
        candidate_hours = tuple(sorted(qpf_forecast_hour(item) for item in candidate_plan.items))
        if candidate_hours == EXPECTED_QPF_FORECAST_HOURS:
            qpf_request = candidate
            selected_cycle_hour = cycle_hour
            break
        rejected_cycles.append({
            "cycle": f"{cycle_hour:02d}Z",
            "missing": sorted(set(EXPECTED_QPF_FORECAST_HOURS) - set(candidate_hours)),
            "unexpected": sorted(set(candidate_hours) - set(EXPECTED_QPF_FORECAST_HOURS)),
        })

    if qpf_request is None:
        raise RuntimeError(
            f"No complete current-day WPC QPF cycle on {current_utc:%Y-%m-%d} UTC; "
            f"candidates={candidate_cycles}, rejected={rejected_cycles}"
        )

    download_result = kit.download_qpf(qpf_request, qpf_grib_dir, force=False)

downloaded_files = tuple(sorted(download_result.all_files, key=qpf_forecast_hour))
downloaded_hours = tuple(qpf_forecast_hour(path) for path in downloaded_files)
if downloaded_hours != EXPECTED_QPF_FORECAST_HOURS:
    raise RuntimeError(f"Incomplete WPC download: {downloaded_hours}")
if download_result.failed:
    raise RuntimeError(f"WPC acquisition failures: {download_result.failed}")

cropped_grib_files = []
crop_rows = []
for source_path in downloaded_files:
    target_path = cropped_grib_dir / source_path.name
    with rasterio.open(source_path) as src:
        projected_bounds = transform_bounds("EPSG:4326", src.crs, *BOUNDS)
        window = from_bounds(*projected_bounds, transform=src.transform)
        window = Window(
            max(0, int(np.floor(window.col_off)) - 3),
            max(0, int(np.floor(window.row_off)) - 3),
            min(src.width, int(np.ceil(window.width)) + 6),
            min(src.height, int(np.ceil(window.height)) + 6),
        )
        values = src.read(1, window=window).astype("float64")
        transform = src.window_transform(window)
        crs = src.crs
    if target_path.exists():
        target_path.unlink()
    with rasterio.open(
        target_path, "w", driver="GRIB", dtype="float64", count=1,
        width=values.shape[1], height=values.shape[0], crs=crs, transform=transform,
    ) as dst:
        dst.write(values, 1)
    with rasterio.open(target_path) as check:
        check_values = check.read(1)
        assert check.shape == values.shape and np.allclose(check_values, values, atol=0.01)
    cropped_grib_files.append(target_path)
    crop_rows.append({
        "forecast_hour": qpf_forecast_hour(source_path),
        "rows": values.shape[0], "columns": values.shape[1],
        "mean_mm": float(np.nanmean(values)), "maximum_mm": float(np.nanmax(values)),
    })

qpf_grid_options = QPFGridOptions(
    part_a="WPC", part_b="BALD_EAGLE_CREEK", part_c="PRECIPITATION",
    part_f="QPF-6H", cell_size=2500.0, overwrite=True,
)
with QPKit(logger=logging.getLogger("qpkit")) as kit:
    with warnings.catch_warnings():
        warnings.filterwarnings(
            "ignore", message="You will likely lose important projection information.*",
            category=UserWarning, module=r"pyproj\.crs\.crs",
        )
        dss_result = kit.qpf_dss.write(
            tuple(cropped_grib_files), qpf_dss_file,
            grid_options=qpf_grid_options, dss_version=7,
        )

if dss_result.failed:
    raise RuntimeError(f"WPC DSS failures: {dss_result.failed}")
if len(dss_result.written) != len(downloaded_files):
    raise RuntimeError(
        f"Downloaded {len(downloaded_files)} files but wrote {len(dss_result.written)} DSS grids"
    )
assert qpf_dss_file.is_file() and qpf_dss_file.stat().st_size > 0

print(f"WPC cycle: {current_utc:%Y-%m-%d} {selected_cycle_hour:02d}Z")
print(f"Rejected incomplete cycles: {rejected_cycles or 'none'}")
print(f"Downloaded GRIB files: {len(downloaded_files)}")
print(f"Native-projection cropped GRIB files: {len(cropped_grib_files)}")
print(f"Written DSS grids: {len(dss_result.written)}")
print(f"DSS file: {rel(qpf_dss_file)}")
display(pd.DataFrame(crop_rows).round(3))

Text Only

WPC cycle: 2026-09-26 00Z
Rejected incomplete cycles: [{'cycle': '06Z', 'missing': [], 'unexpected': [174]}]
Downloaded GRIB files: 28
Native-projection cropped GRIB files: 28
Written DSS grids: 28
DSS file: working/926_u/model/BaldEagleCrkMulti2D_wpc_qpf_qualification/Precipitation/wpc_qpf/wpc_qpf_bald_eagle.dss

	forecast_hour	rows	columns	mean_mm	maximum_mm
0	6	25	31	0.000	0.000
1	12	25	31	0.000	0.000
2	18	25	31	0.000	0.000
3	24	25	31	0.097	1.000
4	30	25	31	0.311	1.250
5	36	25	31	0.072	0.500
6	42	25	31	2.086	3.750
7	48	25	31	1.958	4.875
8	54	25	31	0.181	0.250
9	60	25	31	0.162	0.750
10	66	25	31	0.144	0.250
11	72	25	31	0.057	0.250
12	78	25	31	0.000	0.000
13	84	25	31	0.000	0.000
14	90	25	31	0.000	0.000
15	96	25	31	0.000	0.000
16	102	25	31	0.000	0.000
17	108	25	31	0.000	0.000
18	114	25	31	0.000	0.000
19	120	25	31	0.000	0.000
20	126	25	31	0.000	0.000
21	132	25	31	0.000	0.000
22	138	25	31	0.000	0.000
23	144	25	31	0.000	0.000
24	150	25	31	0.011	0.250
25	156	25	31	0.000	0.000
26	162	25	31	0.037	0.750
27	168	25	31	0.000	0.000

Read back all DSS grids and select a 24-hour model window

Every grid is read through RasDss, checked for units/type/timing, and included in the catalog audit. The model window is chosen by the largest four-record spatial-mean total. This keeps the hydraulic qualification bounded while using an objective, data-driven wet interval from the live forecast.

Python

qpf_catalog = RasDss.get_catalog(qpf_dss_file)
assert len(qpf_catalog) == len(EXPECTED_QPF_FORECAST_HOURS)

grid_rows = []
dss_arrays = []
first_grid = None
for pathname in qpf_catalog["pathname"].astype(str):
    grid = RasDss.read_grid(qpf_dss_file, pathname)
    values = np.asarray(grid["data"], dtype=float)
    finite = values[np.isfinite(values)]
    assert finite.size and float(finite.min()) >= -1e-9
    assert grid["units"].upper() == "MM"
    assert grid["data_type"].upper() == "PER-CUM"
    if first_grid is None:
        first_grid = grid
    dss_arrays.append(values)
    grid_rows.append({
        "pathname": pathname,
        "_values": values,
        "start_time": pd.Timestamp(grid["start_time"]),
        "end_time": pd.Timestamp(grid["end_time"]),
        "mean_mm": float(np.nanmean(values)),
        "maximum_mm": float(np.nanmax(values)),
        "rows": int(grid["shape"][0]),
        "columns": int(grid["shape"][1]),
        "cell_size_m": float(grid["cell_size"]),
    })

dss_audit = pd.DataFrame(grid_rows).sort_values("start_time").reset_index(drop=True)
dss_arrays = dss_audit.pop("_values").tolist()
expected_duration = pd.Timedelta(hours=QPF_INTERVAL_HOURS)
assert ((dss_audit["end_time"] - dss_audit["start_time"]) == expected_duration).all()
assert (dss_audit["start_time"].iloc[1:].reset_index(drop=True)
        == dss_audit["end_time"].iloc[:-1].reset_index(drop=True)).all()

window_mean_totals = np.array([
    float(np.nanmean(np.nansum(np.stack(dss_arrays[i:i + WINDOW_RECORDS]), axis=0)))
    for i in range(len(dss_arrays) - WINDOW_RECORDS + 1)
])
selected_start_index = int(np.argmax(window_mean_totals))
selected_stop_index = selected_start_index + WINDOW_RECORDS
selected_catalog = dss_audit.iloc[selected_start_index:selected_stop_index].copy().reset_index(drop=True)
selected_arrays = dss_arrays[selected_start_index:selected_stop_index]
selected_dss_total = np.nansum(np.stack(selected_arrays), axis=0)
simulation_start = selected_catalog.iloc[0]["start_time"].to_pydatetime()
simulation_end = selected_catalog.iloc[-1]["end_time"].to_pydatetime()
assert simulation_end - simulation_start == timedelta(hours=24)
assert float(np.nanmean(selected_dss_total)) > 0.0

source_rows = []
raw_maps = []
for local_index, source_path in enumerate(downloaded_files[selected_start_index:selected_stop_index]):
    datasets = cfgrib.open_datasets(source_path, backend_kwargs={"indexpath": ""})
    dataset = next(ds for ds in datasets if ds.data_vars)
    preferred = ("tp", "unknown", "APCP_surface")
    variable = next((name for name in preferred if name in dataset.data_vars), next(iter(dataset.data_vars)))
    grid = dataset[variable].squeeze()
    values = np.asarray(grid.values, dtype=float)
    lat = np.asarray(grid.coords["latitude"].values, dtype=float)
    lon = np.asarray(grid.coords["longitude"].values, dtype=float)
    if lat.ndim == 1 and lon.ndim == 1:
        lon, lat = np.meshgrid(lon, lat)
    lon = np.where(lon > 180.0, lon - 360.0, lon)
    mask = (
        (lon >= WPC_BBOX.left_lon) & (lon <= WPC_BBOX.right_lon)
        & (lat >= WPC_BBOX.bottom_lat) & (lat <= WPC_BBOX.top_lat)
    )
    if not mask.any():
        raise RuntimeError(f"Source WPC grid does not intersect model bounds: {source_path.name}")
    raw_maps.append((values, lon, lat, mask, grid.attrs.get("units", "source units")))
    source_rows.append({
        "forecast_hour": qpf_forecast_hour(source_path),
        "source_mean_mm": float(np.nanmean(values[mask])),
        "DSS_mean_mm": float(selected_catalog.iloc[local_index]["mean_mm"]),
    })
    for ds in datasets:
        ds.close()

source_dss_audit = pd.DataFrame(source_rows)
source_dss_audit["difference_mm"] = source_dss_audit["DSS_mean_mm"] - source_dss_audit["source_mean_mm"]
assert np.isfinite(source_dss_audit[["source_mean_mm", "DSS_mean_mm"]].to_numpy()).all()
assert float(source_dss_audit["difference_mm"].abs().max()) <= MAX_SOURCE_DSS_MEAN_DIFFERENCE_MM

dss_display = dss_audit.drop(columns="pathname").copy()
for column in ("mean_mm", "maximum_mm", "cell_size_m"):
    dss_display[column] = dss_display[column].round(3)
display(dss_display)
display(source_dss_audit.round(4))

wettest_local = int(np.argmax(selected_catalog["mean_mm"].to_numpy()))
raw_values, raw_lon, raw_lat, raw_mask, raw_units = raw_maps[wettest_local]
rows, cols = np.where(raw_mask)
row_slice = slice(rows.min(), rows.max() + 1)
col_slice = slice(cols.min(), cols.max() + 1)

fig, (ax_series, ax_raw, ax_dss) = plt.subplots(1, 3, figsize=(18, 5.2), constrained_layout=True)
ax_series.bar(dss_audit["end_time"], dss_audit["mean_mm"] / 25.4,
              width=pd.Timedelta(hours=5), color="#b8c7d9", label="Full forecast")
ax_series.bar(selected_catalog["end_time"], selected_catalog["mean_mm"] / 25.4,
              width=pd.Timedelta(hours=5), color="#2b6f9f", label="Selected 24 hours")
ax_series.set(title="WPC QPF DSS Interval Means", xlabel="Interval end (UTC)", ylabel="Six-hour depth (in)")
ax_series.tick_params(axis="x", rotation=30)
ax_series.grid(axis="y", alpha=0.25)
ax_series.legend()

raw_image = ax_raw.pcolormesh(
    raw_lon[row_slice, col_slice], raw_lat[row_slice, col_slice],
    raw_values[row_slice, col_slice], shading="auto", cmap="Blues",
)
ax_raw.set(title=f"Published WPC GRIB | f{qpf_forecast_hour(downloaded_files[selected_start_index + wettest_local]):03d}",
           xlabel="Longitude", ylabel="Latitude")
fig.colorbar(raw_image, ax=ax_raw, label=f"Six-hour QPF ({raw_units})")

dss_image = ax_dss.imshow(selected_dss_total / 25.4, origin="upper", cmap="Blues")
ax_dss.set(title="Selected DSS 24-Hour Total", xlabel="SHG column", ylabel="SHG row")
fig.colorbar(dss_image, ax=ax_dss, label="Cumulative depth (in)")
plt.show()

print(f"Selected model window: {simulation_start:%Y-%m-%d %H:%M} to {simulation_end:%Y-%m-%d %H:%M} UTC")
print(f"Selected DSS mean total: {np.nanmean(selected_dss_total) / 25.4:.4f} inches")
print(f"Selected DSS maximum total: {np.nanmax(selected_dss_total) / 25.4:.4f} inches")

	start_time	end_time	mean_mm	maximum_mm	rows	columns	cell_size_m
0	2026-09-26 00:00:00	2026-09-26 06:00:00	0.000	0.000	27	32	2500.0
1	2026-09-26 06:00:00	2026-09-26 12:00:00	0.000	0.000	27	32	2500.0
2	2026-09-26 12:00:00	2026-09-26 18:00:00	0.000	0.000	27	32	2500.0
3	2026-09-26 18:00:00	2026-09-27 00:00:00	0.097	1.000	27	32	2500.0
4	2026-09-27 00:00:00	2026-09-27 06:00:00	0.313	1.250	27	32	2500.0
5	2026-09-27 06:00:00	2026-09-27 12:00:00	0.074	0.500	27	32	2500.0
6	2026-09-27 12:00:00	2026-09-27 18:00:00	2.086	3.750	27	32	2500.0
7	2026-09-27 18:00:00	2026-09-28 00:00:00	1.962	4.875	27	32	2500.0
8	2026-09-28 00:00:00	2026-09-28 06:00:00	0.181	0.250	27	32	2500.0
9	2026-09-28 06:00:00	2026-09-28 12:00:00	0.160	0.750	27	32	2500.0
10	2026-09-28 12:00:00	2026-09-28 18:00:00	0.143	0.250	27	32	2500.0
11	2026-09-28 18:00:00	2026-09-29 00:00:00	0.058	0.250	27	32	2500.0
12	2026-09-29 00:00:00	2026-09-29 06:00:00	0.000	0.000	27	32	2500.0
13	2026-09-29 06:00:00	2026-09-29 12:00:00	0.000	0.000	27	32	2500.0
14	2026-09-29 12:00:00	2026-09-29 18:00:00	0.000	0.000	27	32	2500.0
15	2026-09-29 18:00:00	2026-09-30 00:00:00	0.000	0.000	27	32	2500.0
16	2026-09-30 00:00:00	2026-09-30 06:00:00	0.000	0.000	27	32	2500.0
17	2026-09-30 06:00:00	2026-09-30 12:00:00	0.000	0.000	27	32	2500.0
18	2026-09-30 12:00:00	2026-09-30 18:00:00	0.000	0.000	27	32	2500.0
19	2026-09-30 18:00:00	2026-10-01 00:00:00	0.000	0.000	27	32	2500.0
20	2026-10-01 00:00:00	2026-10-01 06:00:00	0.000	0.000	27	32	2500.0
21	2026-10-01 06:00:00	2026-10-01 12:00:00	0.000	0.000	27	32	2500.0
22	2026-10-01 12:00:00	2026-10-01 18:00:00	0.000	0.000	27	32	2500.0
23	2026-10-01 18:00:00	2026-10-02 00:00:00	0.000	0.000	27	32	2500.0
24	2026-10-02 00:00:00	2026-10-02 06:00:00	0.011	0.250	27	32	2500.0
25	2026-10-02 06:00:00	2026-10-02 12:00:00	0.000	0.000	27	32	2500.0
26	2026-10-02 12:00:00	2026-10-02 18:00:00	0.037	0.750	27	32	2500.0
27	2026-10-02 18:00:00	2026-10-03 00:00:00	0.000	0.000	27	32	2500.0

	forecast_hour	source_mean_mm	DSS_mean_mm	difference_mm
0	30	0.2434	0.3128	0.0694
1	36	0.0456	0.0735	0.0280
2	42	2.2219	2.0859	-0.1361
3	48	1.7391	1.9617	0.2227

Text Only

Selected model window: 2026-09-27 00:00 to 2026-09-28 00:00 UTC
Selected DSS mean total: 0.1491 inches
Selected DSS maximum total: 0.3986 inches

Author the baseline and WPC rain-on-grid plans

Both plans retain the same geometry, boundaries, infiltration, and inflow. The forecast plan references the first selected DSS pathname; HEC-RAS discovers the rest of the same A/B/C/F record series by time. ratio=1.0 is explicit because HEC-RAS 6.0-6.1 can otherwise retain an ineffective ratio. The DSS setter completes the required evapotranspiration, air-density, and air-pressure modes without overriding explicit nonblank choices. Runtime-message and HDF review below remains required.

Python

baseline_unsteady = RasPlan.clone_unsteady(
    TEMPLATE_UNSTEADY, new_title="WPC No-Rain Baseline", ras_object=ras
)
forecast_unsteady = RasPlan.clone_unsteady(
    TEMPLATE_UNSTEADY, new_title="WPC QPF Forecast", ras_object=ras
)
plan_intervals = {"computation": COMPUTATION_INTERVAL, "output": "10MIN", "mapping": "10MIN"}
baseline_plan = RasPlan.clone_plan(
    TEMPLATE_PLAN, new_plan_shortid="WPC No Rain", new_title="WPC No-Rain Baseline",
    geometry=TEMPLATE_GEOMETRY, unsteady_flow=baseline_unsteady, num_cores=4,
    intervals=plan_intervals, ras_object=ras,
)
forecast_plan = RasPlan.clone_plan(
    TEMPLATE_PLAN, new_plan_shortid="WPC Forecast", new_title="WPC QPF Forecast",
    geometry=TEMPLATE_GEOMETRY, unsteady_flow=forecast_unsteady, num_cores=4,
    intervals=plan_intervals, ras_object=ras,
)
for plan_number in (baseline_plan, forecast_plan):
    RasPlan.update_simulation_date(plan_number, simulation_start, simulation_end, ras_object=ras)

RasUnsteady.set_constant_precipitation(
    baseline_unsteady, value=0.0, units="in/hr", ras_object=ras
)
first_selected_pathname = str(selected_catalog.iloc[0]["pathname"])
RasUnsteady.configure_gridded_dss_precipitation(
    forecast_unsteady,
    dss_filename=str(qpf_dss_file),
    dss_pathname=first_selected_pathname,
    ras_object=ras,
    ratio=1.0,
    interpolation="Bilinear",
)
baseline_config = RasUnsteady.get_met_precipitation_config(baseline_unsteady, ras_object=ras)
forecast_config = RasUnsteady.get_met_precipitation_config(forecast_unsteady, ras_object=ras)
assert baseline_config["mode"] == "Constant" and baseline_config["constant_value"] == 0.0
assert forecast_config["mode"] == "Gridded" and forecast_config["source"] == "DSS"

display(pd.DataFrame([
    {"condition": "No-rain baseline", "plan": f"p{baseline_plan}", "unsteady": f"u{baseline_unsteady}",
     "precipitation": "Constant 0.0 in/hr"},
    {"condition": "WPC QPF forecast", "plan": f"p{forecast_plan}", "unsteady": f"u{forecast_unsteady}",
     "precipitation": f"{WINDOW_RECORDS} selected six-hour DSS grids"},
]))

	condition	plan	unsteady	precipitation
0	No-rain baseline	p07	u04	Constant 0.0 in/hr
1	WPC QPF forecast	p08	u05	4 selected six-hour DSS grids

Gold-standard precompute validation

RasPreprocess.preprocess_plan() creates the temporary plan HDF consumed by the solver. The values are checked against the selected DSS grids before either hydraulic run begins.

Python

prepared = RasPreprocess.preprocess_plan(forecast_plan, ras_object=ras, max_wait=600)
if not bool(prepared):
    raise RuntimeError(prepared.error)
with h5py.File(prepared.tmp_hdf_path, "r") as hdf:
    solver_rainfall = hdf["Event Conditions/Meteorology/Precipitation/Values"][:]
    solver_timestamps = hdf["Event Conditions/Meteorology/Precipitation/Timestamp"][:]
assert solver_rainfall.ndim == 2 and solver_rainfall.shape[0] == len(solver_timestamps)
assert np.isfinite(solver_rainfall).all() and float(np.nanmin(solver_rainfall)) >= 0.0
solver_time = pd.DatetimeIndex(pd.to_datetime(
    [value.decode("ascii") for value in solver_timestamps],
    format="%d%b%Y %H:%M:%S.%f",
))
# HEC-RAS retains preceding DSS rows as interpolation context. Only rows whose
# interval endpoints are after the plan start and through the plan end are
# accumulated for the solver-window depth comparison.
solver_window_mask = (solver_time > pd.Timestamp(simulation_start)) & (solver_time <= pd.Timestamp(simulation_end))
solver_window_rainfall = solver_rainfall[solver_window_mask]
solver_window_time = solver_time[solver_window_mask]
assert solver_window_rainfall.shape[0] == WINDOW_RECORDS
solver_cumulative = np.cumsum(solver_window_rainfall, axis=0)
solver_total = solver_cumulative[-1]
dss_to_precompute_error = solver_total - np.flipud(selected_dss_total).reshape(-1) / 25.4
assert float(np.nanmax(np.abs(dss_to_precompute_error))) <= 5e-6

fig, ax = plt.subplots(figsize=(11, 4.2))
ax.plot(solver_window_time, np.nanmean(solver_cumulative, axis=1), marker="o", label="Raster mean", color="#2b6f9f")
ax.plot(solver_window_time, np.nanmax(solver_cumulative, axis=1), marker="o", label="Raster maximum", color="#b5453c")
ax.set(title="HEC-RAS Temporary-HDF Cumulative WPC Precipitation", xlabel="Model time (UTC)", ylabel="Cumulative depth (in)")
ax.grid(alpha=0.25)
ax.legend()
fig.autofmt_xdate()
plt.show()

display(pd.DataFrame([{
    "temporary-HDF rows including context": solver_rainfall.shape[0],
    "solver-window rows": solver_window_rainfall.shape[0],
    "DSS raster positions": solver_rainfall.shape[1],
    "final mean depth (in)": float(np.nanmean(solver_total)),
    "final maximum depth (in)": float(np.nanmax(solver_total)),
    "maximum DSS-to-precompute error (in)": float(np.nanmax(np.abs(dss_to_precompute_error))),
}]).round(7))

	temporary-HDF rows including context	solver-window rows	DSS raster positions	final mean depth (in)	final maximum depth (in)	maximum DSS-to-precompute error (in)
0	6	4	864	0.149107	0.398622	2.000000e-07

Compute the no-rain and WPC plans

The BCO files and plan-HDF compute messages are both audited. Automated checks are gates and evidence, not a substitute for reviewing the full messages in HEC-RAS and checking the model manually.

Python

hdf_output_variables = [
    "Cell Hydraulic Depth", "Water Surface", "Cell Precipitation Rate",
    "Cell Cumulative Precipitation Depth",
]
print(f"Computing p{baseline_plan} no-rain baseline in HEC-RAS {RAS_VERSION}...")
baseline_result = RasCmdr.compute_plan(
    baseline_plan, ras_object=ras, force_rerun=True, force_geompre=True, verify=True,
    use_optimal_hdf_settings=True, hdf_output_variables=hdf_output_variables,
)
if not bool(baseline_result):
    raise RuntimeError("No-rain baseline compute failed")
print(f"Baseline p{baseline_plan}: complete")

print(f"Computing p{forecast_plan} WPC QPF forecast in HEC-RAS {RAS_VERSION}...")
forecast_result = RasCmdr.compute_plan(
    forecast_plan, ras_object=ras, force_rerun=True, force_geompre=True, verify=True,
    use_optimal_hdf_settings=True, hdf_output_variables=hdf_output_variables,
)
if not bool(forecast_result):
    raise RuntimeError("WPC forecast compute failed")
print(f"Forecast p{forecast_plan}: complete")

baseline_hdf = ras.project_folder / f"{ras.project_name}.p{baseline_plan}.hdf"
forecast_hdf = ras.project_folder / f"{ras.project_name}.p{forecast_plan}.hdf"
diagnostic_pattern = re.compile(
    r"\b(?:warning|fatal|failed|missing)\b|"
    r"\berror(?::|\s+(?:reading|opening|processing|computing|occurred|encountered))",
    flags=re.IGNORECASE,
)

def runtime_audit(condition: str, plan_number: str, hdf_path: Path) -> dict:
    messages = HdfResultsPlan.get_compute_messages_hdf_only(hdf_path)
    ws_error = HdfResultsMesh.get_mesh_max_ws_err(hdf_path)
    error_column = "cell_maximum_water_surface_error"
    volume = HdfResultsPlan.get_volume_accounting(hdf_path)
    bco_path = ras.project_folder / f"{ras.project_name}.bco{plan_number}"
    assert bco_path.exists()
    bco_text = bco_path.read_text(encoding="utf-8", errors="replace")
    problem_lines = [line.strip() for line in bco_text.splitlines() if diagnostic_pattern.search(line)]
    lower = messages.lower()
    return {
        "condition": condition,
        "plan": f"p{plan_number}",
        "complete_process": "complete process" in lower,
        "precipitation_coverage_warning": "out-of-bounds" in lower and ("2d cells" in lower or "2d faces" in lower),
        "BCO diagnostic lines": len(problem_lines),
        "maximum cell WSEL error (ft)": float(ws_error[error_column].max()),
        "cells WSEL error >0.1 ft": int(np.count_nonzero(ws_error[error_column] > 0.1)),
        "volume error (%)": float(volume.iloc[0]["Error Percent"]) if volume is not None and not volume.empty else np.nan,
    }

compute_audit = pd.DataFrame([
    runtime_audit("No-rain baseline", baseline_plan, baseline_hdf),
    runtime_audit("WPC QPF forecast", forecast_plan, forecast_hdf),
])
display(compute_audit.round(6))
assert compute_audit["complete_process"].all()
assert not compute_audit["precipitation_coverage_warning"].any()
assert (compute_audit["BCO diagnostic lines"] == 0).all()
assert float(compute_audit["maximum cell WSEL error (ft)"].max()) <= MAX_WSEL_ERROR_FT
assert int(compute_audit["cells WSEL error >0.1 ft"].sum()) <= MAX_CELLS_WSEL_ERROR_GT_0_1_FT
if int(compute_audit["cells WSEL error >0.1 ft"].sum()):
    print(
        "Localized WSEL convergence outliers remain in this instructional model; "
        "they are mapped below and require manual review."
    )

Text Only

Computing p07 no-rain baseline in HEC-RAS 7.0...


Baseline p07: complete
Computing p08 WPC QPF forecast in HEC-RAS 7.0...


Forecast p08: complete

	condition	plan	complete_process	precipitation_coverage_warning	BCO diagnostic lines	maximum cell WSEL error (ft)	cells WSEL error >0.1 ft	volume error (%)
0	No-rain baseline	p07	True	False	0	0.00999	0	0.000030
1	WPC QPF forecast	p08	True	False	0	0.87922	3	0.000332

Text Only

Localized WSEL convergence outliers remain in this instructional model; they are mapped below and require manual review.

Completed-plan precipitation validation

The completed plan HDF is the post-compute forcing record. A nearest-DSS-cell comparison is retained alongside the temporary-HDF exact check so interpolation differences remain visible rather than being hidden.

Python

def precipitation_series(plan_hdf: Path):
    rate = HdfResultsMesh.get_mesh_timeseries(
        plan_hdf, MESH_NAME, "Cell Precipitation Rate", truncate=False
    )
    cumulative = HdfResultsMesh.get_mesh_timeseries(
        plan_hdf, MESH_NAME, "Cell Cumulative Precipitation Depth", truncate=False
    )
    return rate, cumulative

baseline_rate, baseline_cumulative = precipitation_series(baseline_hdf)
forecast_rate, forecast_cumulative = precipitation_series(forecast_hdf)
baseline_rate_values = np.asarray(baseline_rate.values, dtype=float)
baseline_cumulative_values = np.asarray(baseline_cumulative.values, dtype=float)
forecast_rate_values = np.asarray(forecast_rate.values, dtype=float)
forecast_cumulative_values = np.asarray(forecast_cumulative.values, dtype=float)
assert np.isclose(np.nanmax(baseline_rate_values), 0.0, atol=1e-9)
assert np.isclose(np.nanmax(baseline_cumulative_values), 0.0, atol=1e-9)
assert float(np.nanmax(forecast_rate_values)) > 0.0
assert float(np.nanmax(forecast_cumulative_values)) > 0.0

final_cumulative = forecast_cumulative_values[-1]
mesh_points = HdfMesh.get_mesh_cell_points(forecast_hdf)
mesh_points = mesh_points.loc[mesh_points["mesh_name"] == MESH_NAME].sort_values("cell_id").reset_index(drop=True)
assert len(mesh_points) == final_cumulative.size

origin_x, origin_y = first_grid["metadata"]["origin"]
cell_size = float(first_grid["cell_size"])
rows, columns = first_grid["shape"]
mesh_in_dss = mesh_points.to_crs(first_grid["crs"])
point_x = mesh_in_dss.geometry.x.to_numpy()
point_y = mesh_in_dss.geometry.y.to_numpy()
sample_columns = np.floor((point_x - origin_x) / cell_size).astype(int)
sample_rows = np.floor((point_y - origin_y) / cell_size).astype(int)
inside_dss = (
    (sample_columns >= 0) & (sample_columns < columns)
    & (sample_rows >= 0) & (sample_rows < rows)
)
sampled_dss_mm = np.full(len(mesh_in_dss), np.nan, dtype=float)
sampled_dss_mm[inside_dss] = selected_dss_total[sample_rows[inside_dss], sample_columns[inside_dss]]
sampled_dss_in = sampled_dss_mm / 25.4
active_mesh_mask = final_cumulative > 1e-6
# HEC-RAS retains result-center slots for 1,531 inactive cells in this fixture.
# They remain exactly zero and are not computational cells, so forcing
# comparisons use the positive-rain active positions and verify that their
# count matches the geometry's Cell Count attribute.
assert int(np.count_nonzero(active_mesh_mask)) == mesh_cell_count
comparison_mask = (
    inside_dss & active_mesh_mask
    & np.isfinite(sampled_dss_in) & np.isfinite(final_cumulative)
)
comparison_error = final_cumulative[comparison_mask] - sampled_dss_in[comparison_mask]
assert comparison_mask.any()
assert abs(float(np.nanmean(comparison_error))) <= 0.02
assert float(np.sqrt(np.nanmean(comparison_error ** 2))) <= 0.08

forcing_summary = pd.DataFrame([{
    "DSS sampled mean (in)": float(np.nanmean(sampled_dss_in[comparison_mask])),
    "temporary-HDF raster mean (in)": float(np.nanmean(solver_total)),
    "final-HDF mesh mean (in)": float(np.nanmean(final_cumulative[comparison_mask])),
    "DSS-to-final mean error (in)": float(np.nanmean(comparison_error)),
    "DSS-to-final RMSE (in)": float(np.sqrt(np.nanmean(comparison_error ** 2))),
    "active positive-rain mesh cells": int(np.count_nonzero(active_mesh_mask)),
}])
display(forcing_summary.round(6))

fig, (ax_scatter, ax_hist) = plt.subplots(1, 2, figsize=(13, 5), constrained_layout=True)
ax_scatter.scatter(sampled_dss_in[comparison_mask], final_cumulative[comparison_mask], s=8, alpha=0.4)
limits = [float(np.nanmin(sampled_dss_in[comparison_mask])), float(np.nanmax(sampled_dss_in[comparison_mask]))]
ax_scatter.plot(limits, limits, color="black", linestyle="--", label="1:1")
ax_scatter.set(title="WPC DSS vs Final-HDF Cumulative Depth", xlabel="Nearest DSS grid depth (in)", ylabel="HEC-RAS final depth (in)")
ax_scatter.grid(alpha=0.25)
ax_scatter.legend()
ax_hist.hist(comparison_error, bins=35, color="#2b6f9f", edgecolor="white")
ax_hist.axvline(0.0, color="black", linestyle="--")
ax_hist.set(title="Interpolation Difference", xlabel="Final HDF minus nearest DSS cell (in)", ylabel="Mesh cells")
plt.show()

	DSS sampled mean (in)	temporary-HDF raster mean (in)	final-HDF mesh mean (in)	DSS-to-final mean error (in)	DSS-to-final RMSE (in)	active positive-rain mesh cells
0	0.163637	0.149107	0.163409	-0.000229	0.003105	18066

Hydraulic response and convergence review

The rainfall case is compared with the no-rain control at every mesh cell. The map includes cells whose maximum water-surface error exceeds 0.1 ft so a reviewer can distinguish forcing response from local numerical concerns.

Python

baseline_depth_ts = HdfResultsMesh.get_mesh_timeseries(
    baseline_hdf, MESH_NAME, "Cell Hydraulic Depth", truncate=False
)
forecast_depth_ts = HdfResultsMesh.get_mesh_timeseries(
    forecast_hdf, MESH_NAME, "Cell Hydraulic Depth", truncate=False
)
baseline_depth_values = np.asarray(baseline_depth_ts.values, dtype=float)
forecast_depth_values = np.asarray(forecast_depth_ts.values, dtype=float)
assert baseline_depth_values.shape == forecast_depth_values.shape
baseline_depth_values[baseline_depth_values < 0] = np.nan
forecast_depth_values[forecast_depth_values < 0] = np.nan
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", message="All-NaN slice encountered", category=RuntimeWarning)
    baseline_depth_max = np.nanmax(baseline_depth_values, axis=0)
    forecast_depth_max = np.nanmax(forecast_depth_values, axis=0)

depth_comparison = mesh_points.copy()
depth_comparison["baseline_max_depth_ft"] = baseline_depth_max
depth_comparison["forecast_max_depth_ft"] = forecast_depth_max
depth_comparison["maximum_depth_envelope_change_ft"] = forecast_depth_max - baseline_depth_max
depth_delta = depth_comparison["maximum_depth_envelope_change_ft"].replace([np.inf, -np.inf], np.nan)
assert float(depth_delta.max()) > 0.0

representative_row = depth_comparison.loc[depth_delta.idxmax()]
representative_cell = int(representative_row["cell_id"])
baseline_times = pd.DatetimeIndex(pd.to_datetime(baseline_depth_ts.coords["time"].values))
forecast_times = pd.DatetimeIndex(pd.to_datetime(forecast_depth_ts.coords["time"].values))

hydraulic_summary = pd.DataFrame([{
    "baseline cells ever wet >0.1 ft": int(np.count_nonzero(baseline_depth_max > 0.1)),
    "forecast cells ever wet >0.1 ft": int(np.count_nonzero(forecast_depth_max > 0.1)),
    "newly wet cells": int(np.count_nonzero((forecast_depth_max > 0.1) & ~(baseline_depth_max > 0.1))),
    "maximum depth-envelope increase (ft)": float(depth_delta.max()),
    "cells with depth increase >0.1 ft": int(np.count_nonzero(depth_delta > 0.1)),
}])
display(hydraulic_summary.round(4).T.rename(columns={0: "value"}))

fig, (ax_ts, ax_delta) = plt.subplots(1, 2, figsize=(15, 5.8), constrained_layout=True)
ax_ts.plot(baseline_times, baseline_depth_values[:, representative_cell], color="#2b6f9f", linewidth=2, label="No-rain baseline")
ax_ts.plot(forecast_times, forecast_depth_values[:, representative_cell], color="#b5453c", linewidth=2, label="WPC QPF")
ax_ts.set(title=f"Hydraulic Depth at Cell {representative_cell}", xlabel="Model time (UTC)", ylabel="Depth (ft)")
ax_ts.tick_params(axis="x", rotation=30)
ax_ts.grid(alpha=0.25)
ax_ts.legend()

finite_abs = np.abs(depth_delta.dropna())
color_limit = max(float(finite_abs.max()), 0.01)
depth_comparison.plot(
    ax=ax_delta, column="maximum_depth_envelope_change_ft", cmap="RdBu_r",
    vmin=-color_limit, vmax=color_limit, markersize=5, legend=True,
    legend_kwds={"label": "WPC minus baseline max-depth envelope (ft)"},
)
mesh_areas.boundary.plot(ax=ax_delta, color="#263238", linewidth=0.8)
forecast_ws_error = HdfResultsMesh.get_mesh_max_ws_err(forecast_hdf)
high_error_points = forecast_ws_error.loc[
    forecast_ws_error["cell_maximum_water_surface_error"] > 0.1
]
if not high_error_points.empty:
    high_error_points.plot(
        ax=ax_delta, color="#ffcc00", edgecolor="black", marker="o",
        markersize=45, label="WSEL error >0.1 ft",
    )
ax_delta.scatter(representative_row.geometry.x, representative_row.geometry.y,
                 color="black", marker="x", s=55, label=f"Cell {representative_cell}")
ax_delta.set(title="Calculated Hydraulic Response to WPC QPF", xlabel="Easting (ft)", ylabel="Northing (ft)")
ax_delta.set_aspect("equal")
ax_delta.legend(loc="best")
plt.show()

	value
baseline cells ever wet >0.1 ft	2631.0000
forecast cells ever wet >0.1 ft	4609.0000
newly wet cells	1978.0000
maximum depth-envelope increase (ft)	2.3557
cells with depth increase >0.1 ft	2174.0000

Completion record

The checks below establish the executed evidence retained in this notebook. They do not replace engineering review of the model or the current forecast.

Python

completion_checks = pd.Series({
    "complete current-day WPC cycle downloaded": downloaded_hours == EXPECTED_QPF_FORECAST_HOURS,
    "all 28 DSS grids read back": len(dss_audit) == len(EXPECTED_QPF_FORECAST_HOURS),
    "selected 24-hour window has positive rain": float(np.nanmean(selected_dss_total)) > 0.0,
    "source and DSS selected records inspected": len(source_dss_audit) == WINDOW_RECORDS,
    "source-to-DSS regional means agree": float(source_dss_audit["difference_mm"].abs().max()) <= MAX_SOURCE_DSS_MEAN_DIFFERENCE_MM,
    "temporary plan HDF contains exact selected DSS rain": float(np.nanmax(np.abs(dss_to_precompute_error))) <= 5e-6,
    "baseline and forecast computes completed": bool(baseline_result) and bool(forecast_result),
    "runtime/BCO precipitation gates pass": (
        compute_audit["complete_process"].all()
        and not compute_audit["precipitation_coverage_warning"].any()
        and (compute_audit["BCO diagnostic lines"] == 0).all()
    ),
    "localized convergence exceptions are bounded and mapped": (
        float(compute_audit["maximum cell WSEL error (ft)"].max()) <= MAX_WSEL_ERROR_FT
        and int(compute_audit["cells WSEL error >0.1 ft"].sum()) <= MAX_CELLS_WSEL_ERROR_GT_0_1_FT
    ),
    "baseline final HDF is dry": np.isclose(np.nanmax(baseline_cumulative_values), 0.0, atol=1e-9),
    "forecast final HDF is wet": float(np.nanmax(final_cumulative)) > 0.0,
    "DSS and final-HDF spatial depths agree": (
        abs(float(np.nanmean(comparison_error))) <= 0.02
        and float(np.sqrt(np.nanmean(comparison_error ** 2))) <= 0.08
    ),
    "forecast creates a positive hydraulic response": float(depth_delta.max()) > 0.0,
}, name="passed").to_frame()
display(completion_checks)
assert completion_checks["passed"].all()

print("WPC QPF DSS-to-HEC-RAS rain-on-grid qualification complete.")
print(f"Cycle: {current_utc:%Y-%m-%d} {selected_cycle_hour:02d}Z; HEC-RAS {RAS_VERSION}; plans p{baseline_plan}/p{forecast_plan}")
print(f"Selected window: {simulation_start:%Y-%m-%d %H:%M} to {simulation_end:%Y-%m-%d %H:%M} UTC")
print("Localized WSEL convergence exceptions are mapped; this notebook does not certify uniform convergence.")
print("Review all HEC-RAS runtime and BCO messages, precipitation maps, timing, units, terrain coverage, mass balance, convergence, and hydraulic plausibility manually.")

	passed
complete current-day WPC cycle downloaded	True
all 28 DSS grids read back	True
selected 24-hour window has positive rain	True
source and DSS selected records inspected	True
source-to-DSS regional means agree	True
temporary plan HDF contains exact selected DSS rain	True
baseline and forecast computes completed	True
runtime/BCO precipitation gates pass	True
localized convergence exceptions are bounded and mapped	True
baseline final HDF is dry	True
forecast final HDF is wet	True
DSS and final-HDF spatial depths agree	True
forecast creates a positive hydraulic response	True

Text Only

WPC QPF DSS-to-HEC-RAS rain-on-grid qualification complete.
Cycle: 2026-09-26 00Z; HEC-RAS 7.0; plans p07/p08
Selected window: 2026-09-27 00:00 to 2026-09-28 00:00 UTC
Localized WSEL convergence exceptions are mapped; this notebook does not certify uniform convergence.
Review all HEC-RAS runtime and BCO messages, precipitation maps, timing, units, terrain coverage, mass balance, convergence, and hydraulic plausibility manually.