917 - MRMS QPE Boundary-Hyetograph Hydraulic Comparison¶
This notebook acquires two archived MRMS QPE events, converts the source grids to native HEC-DSS for independent spatial inspection, derives a spatial arithmetic mean over each configured bounding-box clip, and compares zero-rain and precipitation-boundary runs in HEC-RAS 7.0. It samples the five-minute hydraulic results and creates review figures and MP4 animations.
Absolute-time qualification: the saved outputs freshly inspect and visualize four native runs completed on 2026-09-27 using the corrected precipitation writer. The original producer executed Davis baseline/event and New Orleans baseline; a separate archived continuation executed New Orleans event after resolving the geometry-support gate. No solve was repeated to repair figures. Both boundaries were verified to use simulation-start timing. MRMS valid times are UTC; these examples map them directly to naive HEC-RAS model-clock labels without converting through the Windows local timezone. The simulation dates remain unchanged: Davis 2022-12-31 00:00 through 2023-01-01 18:00, and NewOrleansMetro 2024-04-10 00:00 through 2024-04-11 06:00.
The source/DSS audit retains the frame ending at each simulation start, but that frame covers the preceding hour and is excluded from forcing. Retained rainfall ends are rebased to the actual boundary origin before writing. Source precision is preserved; the writer inserts one zero-time anchor and rounds cumulative mass to hundredths of an inch. Readback checks every cumulative rounding difference against 0.005 inches. After the last observed QPE interval, an explicit zero-rain tail extends to simulation end; that tail is a modeling assumption, not an MRMS observation. Source-only rainfall maps/animations use the observed interval window; the combined display explicitly marks the model dry tail.
Every selected-mesh native cumulative-depth value is checked against the serialized mass curve at the actual absolute five-minute output timestamps and every physical cell identified independently from geometry polygons and validated native endpoint cycles. The remaining one-face boundary cells are verified through mesh topology and must stay zero. Baseline rate and cumulative arrays must be zero in every selected-mesh cell. The native floating-point tolerance is separate from input rounding and from the source precision check (four source-type epsilons for float32, eight for float64, scaled by depth, with exact dry zeros). Rate arrays are inspected for finite, nonnegative values; the event qualification is based on the full cumulative curve.
The GRIB files use local product identifiers that ecCodes labels with unknown units. Each raw identifier and valid time is checked against the filename and loaded stack. NOAA's operational MRMS table identifies product 209/6/37 as MultiSensor_QPE_01H_Pass2, one-hour accumulation in millimeters. Every source frame has the same finite spatial support. Original files, hashes, aligned/readback CSVs, full native cumulative diagnostics, baseline HDF copies, input snapshots, and public-API runtime receipts remain under the printed run root. Earlier exact notebook bytes are archived there separately.
New Orleans cell 16058 is absent from the detailed polygon API because its detailed face perimeter self-crosses. Its closed native endpoint cycle is independently valid, positive-area, and contains its native center; all 19,711 multi-face cells and 775 one-face boundary cells pass the topology checks. This supports rainfall classification only: it does not repair or replace the detailed geometry. See issue #438. Davis's non-target DS Channel mesh has no precipitation arrays in either HDF; absence does not prove zero rainfall there, and that mesh remains outside this forcing qualification.
All final depth graphics use consistent nearest-cell HDF rasterization without network basemaps. Static maps use their actual maximum; videos explicitly label the public API's positive-value 98th-percentile color cap. The combined panel uses the observed hourly interval covering each flood time, with the preceding interval at an exact endpoint and explicit zero at the origin and after the observed window. Those repeated display frames are never used as new observations or fed to the source-only cumulative animation. See the display mapping CSV and issue #440.
This is deliberately a uniform precipitation-boundary workflow. It does not apply the DSS grids through global gridded meteorology, nor does the rectangular grid mean constitute an area-weighted watershed average. Use notebook 924 for spatial MRMS rain-on-grid and notebook 916 for native gridded-DSS qualification. These two examples qualify the demonstrated forcing/time integration; they do not establish calibrated flood predictions, universal model correctness, or a qualified release across HEC-RAS versions. Follow-up context is recorded in issue #398.
Always inspect the source grids, boundary hyetograph, HEC-RAS runtime messages, precipitation outputs, and hydraulic results manually before engineering use.
The New Orleans event HDF also reports a pipe-network iteration row at 10APR2024 14:32:28, Base Node 426 (Mine Blvd - Node 1), with raw ERROR column 0.450. Its reported volume-accounting error is about 8.040 acre-ft (0.1223%). Complete retained HDF messages are displayed below. These numerical diagnostics remain open for engineering review; successful forcing integration does not establish hydraulic acceptance.
Required source version: use an installation or checkout containing the unreleased precipitation-writer repair in PR #425, or a later successor. PyPI/package metadata 0.103.0 alone does not establish that the repair is present. Imported module hashes are recorded, and boundary readback must pass before compute.
The inspection host also emitted a Java tzmappings format diagnostic, retained in final_execution.log. Source, DSS, and native absolute-time checks passed despite that host diagnostic; it is separate from HEC-RAS numerical messages and remains a follow-up limit.
Depth definitions: both event HDFs lack a native Depth dataset, so the public raster exporter uses max(Water Surface − Cells Minimum Elevation, 0): invert depth in feet. Raster maps and flood/combined videos show that quantity; their generic API colorbar “Depth (ft)” refers to invert depth. Reference-cell hydrographs use the separate native Cell Hydraulic Depth variable. The summary labels both quantities explicitly and also reports native invert maxima; nearest-cell sampling can omit a peak cell, so raster/native maxima need not be equal. The wet/dry comparison checks raster invert depth against native hydraulic depth at a 0.01-ft threshold; it does not establish equality of the two depth definitions.
Series position: Start with 915_realtime_forecast_workflow.ipynb for the
operational architecture. Use this notebook when a spatial MRMS product must be
reduced intentionally to a uniform precipitation boundary; use
924_mrms_netcdf_rain_on_grid.ipynb when spatial variability must reach the
HEC-RAS solver. The two pathways are not interchangeable.
Python
from contextlib import contextmanager
from datetime import datetime, timedelta, timezone
from pathlib import Path
import logging
import hashlib
import inspect
import json
import sys
import time
import os
import re
import shutil
import warnings
# tqdm.auto emits this during import when Jupyter widgets are unavailable.
warnings.filterwarnings(
"ignore", message="IProgress not found.*", module="tqdm.auto"
)
import geopandas as gpd
from IPython.display import display
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from rasterio.errors import NotGeoreferencedWarning
from shapely.geometry import box
from ras_commander import (
BcoMonitor,
RasPrj,
RasCmdr,
RasExamples,
RasPlan,
RasProcess,
RasUnsteady,
init_ras_project,
)
from ras_commander.dss import RasDss
from ras_commander.hdf import (
HdfMesh, HdfProject, HdfPump, HdfResultsMesh, HdfResultsPlan,
)
from ras_commander.precip import PrecipMrms
ipython = get_ipython()
if ipython is not None:
ipython.run_line_magic("matplotlib", "inline")
REPO_ROOT = Path.cwd()
if REPO_ROOT.name.lower() == "examples":
REPO_ROOT = REPO_ROOT.parent
RUN_ROOT = Path(os.environ.get(
"RAS_COMMANDER_EXAMPLE_RUN_ROOT",
REPO_ROOT / "working" / f"917_{datetime.now():%Y%m%dT%H%M%S}",
))
ARTIFACT_ROOT = RUN_ROOT / "videos_hdf"
ARTIFACT_ROOT.mkdir(parents=True, exist_ok=True)
OUTPUT_INTERVAL = "5MIN"
FPS = 4
VIDEO_DPI = 110
MIN_CONSOLIDATION_CELL_SIZE_FT = 10.0
ANIMATION_MAX_FRAMES = 30
# Preserve genuine runtime diagnostics in the executed notebook.
DAVIS_BOUNDS = (
-121.78689244974507,
38.523870898941006,
-121.70713010751233,
38.591248584093314,
)
CASES = [
{
"case_id": "davis_ar",
"display_name": "Davis atmospheric river",
"project": "Davis",
"expected_physical_cells": 2705,
"expected_boundary_cells": 236,
"max_runtime": 1200,
"ras_version": "7.0",
"plan_number": "02",
"precip_boundary": "area2",
"event_start": datetime(2022, 12, 31, 0),
"event_last_qpe": datetime(2023, 1, 1, 12),
"sim_end": datetime(2023, 1, 1, 18),
"animation_start": datetime(2022, 12, 31, 0),
"animation_end": datetime(2023, 1, 1, 18),
"mrms_bounds": DAVIS_BOUNDS,
"dss_b": "DAVIS_AR",
"dss_f": "MRMS_AR_20221231",
"event_note": (
"29 Dec 2022 - 1 Jan 2023 Northern California atmospheric river; "
"CLB-672 selected this Davis-local MRMS window."
),
},
{
"case_id": "nola_apr24",
"display_name": "NewOrleansMetro April 2024 flash flood",
"project": "NewOrleansMetro",
"expected_physical_cells": 19711,
"expected_boundary_cells": 775,
"max_runtime": 2400,
"ras_version": "7.0",
"plan_number": "01",
"precip_boundary": "NewOrleans Metro",
"event_start": datetime(2024, 4, 10, 0),
"event_last_qpe": datetime(2024, 4, 10, 23),
"sim_end": datetime(2024, 4, 11, 6),
"animation_start": datetime(2024, 4, 10, 0),
"animation_end": datetime(2024, 4, 11, 6),
"mrms_bounds": None,
"dss_b": "NOLA_APR10",
"dss_f": "MRMS_20240410",
"event_note": (
"NWS New Orleans/Baton Rouge documented the 10 Apr 2024 New Orleans "
"Metro severe thunderstorm and flash-flood event as a non-hurricane "
"rainfall case."
),
},
]
def display_path(path: Path) -> str:
path = Path(path).resolve()
try:
return path.relative_to(REPO_ROOT.resolve()).as_posix()
except ValueError:
return path.name
print(f"Run root: {display_path(RUN_ROOT)}")
print(f"Animation output root: {display_path(ARTIFACT_ROOT)}")
print(f"RAS output/mapping interval target: {OUTPUT_INTERVAL}")
print(
f"Consolidation policy: preserve dominant terrain cells >=10 ft; "
f"otherwise use {MIN_CONSOLIDATION_CELL_SIZE_FT:g} ft; maximum frames: {ANIMATION_MAX_FRAMES}"
)
# Record actual imported files; repository identity is an optional development guard.
import ras_commander
SOURCE_MODULES = {
"ras_commander/RasUnsteady.py": RasUnsteady,
"ras_commander/precip/PrecipMrms.py": PrecipMrms,
"ras_commander/RasCmdr.py": RasCmdr,
"ras_commander/RasPrj.py": RasPrj,
"ras_commander/hdf/HdfMesh.py": HdfMesh,
"ras_commander/hdf/HdfResultsMesh.py": HdfResultsMesh,
}
SOURCE_PATHS = {name: Path(inspect.getfile(klass)).resolve() for name, klass in SOURCE_MODULES.items()}
SOURCE_SHA256 = {name: hashlib.sha256(path.read_bytes()).hexdigest() for name, path in SOURCE_PATHS.items()}
if os.environ.get("RAS_COMMANDER_VERIFY_SOURCE_CHECKOUT") == "1":
for name, path in SOURCE_PATHS.items():
assert path == (REPO_ROOT / name).resolve(), "Imported source differs from requested checkout"
EXECUTION_MODE = os.environ.get("RAS_COMMANDER_917_MODE", "compute")
assert EXECUTION_MODE in {"compute", "inspect_retained"}, "Choose compute or inspect_retained explicitly"
print(f"EXECUTION MODE: {EXECUTION_MODE}")
print(f"Python: {sys.version}; executable: {sys.executable}; ras-commander: {ras_commander.__version__}")
print("MRMS UTC timestamps are used directly as naive HEC-RAS model-clock labels; no local-time conversion.")
print("Imported source paths:", json.dumps({name: str(path) for name, path in SOURCE_PATHS.items()}, indent=2))
print("Source SHA256:", json.dumps(SOURCE_SHA256, indent=2))
plt.rcParams["animation.ffmpeg_args"] = ["-threads", "2"]
plt.rcParams["axes.formatter.useoffset"] = False
Text Only
Run root: working/917_alignment_20260927
Animation output root: working/917_alignment_20260927/videos_hdf
RAS output/mapping interval target: 5MIN
Consolidation policy: preserve dominant terrain cells >=10 ft; otherwise use 10 ft; maximum frames: 30
EXECUTION MODE: inspect_retained
Python: 3.11.13 (main, Jul 11 2025, 22:36:59) [MSC v.1944 64 bit (AMD64)]; executable: C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\.venv\Scripts\python.exe; ras-commander: 0.103.0
MRMS UTC timestamps are used directly as naive HEC-RAS model-clock labels; no local-time conversion.
Imported source paths: {
"ras_commander/RasUnsteady.py": "C:\\Users\\billk_clb\\.codex\\workspaces\\ras-docs-priority-20260926\\ras_commander\\RasUnsteady.py",
"ras_commander/precip/PrecipMrms.py": "C:\\Users\\billk_clb\\.codex\\workspaces\\ras-docs-priority-20260926\\ras_commander\\precip\\PrecipMrms.py",
"ras_commander/RasCmdr.py": "C:\\Users\\billk_clb\\.codex\\workspaces\\ras-docs-priority-20260926\\ras_commander\\RasCmdr.py",
"ras_commander/RasPrj.py": "C:\\Users\\billk_clb\\.codex\\workspaces\\ras-docs-priority-20260926\\ras_commander\\RasPrj.py",
"ras_commander/hdf/HdfMesh.py": "C:\\Users\\billk_clb\\.codex\\workspaces\\ras-docs-priority-20260926\\ras_commander\\hdf\\HdfMesh.py",
"ras_commander/hdf/HdfResultsMesh.py": "C:\\Users\\billk_clb\\.codex\\workspaces\\ras-docs-priority-20260926\\ras_commander\\hdf\\HdfResultsMesh.py"
}
Source SHA256: {
"ras_commander/RasUnsteady.py": "9d20063dd04e08256a7c480f24342551e61d90d30e1b5c44bcfd457662645f5d",
"ras_commander/precip/PrecipMrms.py": "2fb0df862ed121d3f31f765893cf710222b6a6b0c6195b863bf16733e54c0be4",
"ras_commander/RasCmdr.py": "9fd96c94a5165dd69342605b5cf1d9a975d28fd83ddc2d3162ee48cde0fbe50e",
"ras_commander/RasPrj.py": "adb984fd94bbfe0bf8cb9ac76816aac8c74a505ecadde1c97364188049cb4176",
"ras_commander/hdf/HdfMesh.py": "cb5c56b49b2f5f3a4e039ddd50851a065a5979a595e54cadeffba84d6f1ca93a",
"ras_commander/hdf/HdfResultsMesh.py": "44a8cdf45140f05278dc028f9a2489d5db7563973e317bceab941714cf46896d"
}
Development Mode¶
For local source testing, run this notebook from the repository root using the repository .venv and an editable ras-commander install. Generated projects remain under the ignored working/ directory.
Use a new RAS_COMMANDER_EXAMPLE_RUN_ROOT for every compute execution; inspect_retained requires the exact preserved root. This notebook refuses to overwrite a root containing prior compute receipts. The original native sources are identified by the receipt hashes; setup and the emitted summary record the actual inspection Python/library identities and imported module hashes. Package version metadata alone is not a release qualification. Each native run is sequential and limited to two cores: 1,200 seconds per Davis plan and 2,400 seconds per NewOrleans plan.
Workflow Scope¶
The native phase computes paired zero-rain and MRMS uniform precipitation-boundary runs. Retained inspection independently checks source timestamps/depths, exact native input snapshots, source and HDF hashes, complete absolute output windows, geometry-derived rainfall support, and all native cumulative values. The saved run uses explicit inspect_retained mode after the four separately recorded computations.
The visualization phase samples up to 30 timestamps across each five-minute result window, exports every frame through the same public HDF nearest-cell API, and preserves the stated resolution. These diagnostic rasters are not native Mapper interpolation equivalents or engineering flood maps. Raw DSS grids remain spatial products for independent inspection; the solver consumes a uniform arithmetic-mean hyetograph over each configured rectangular clip.
Shared Helpers¶
Python
def hdf_invert_depth_maxima(plan_hdf: Path, geom_hdf: Path, mesh_name: str) -> pd.Series:
"""Mirror the verified exporter fallback through public HDF readers."""
try:
HdfResultsMesh.get_mesh_timeseries(plan_hdf, mesh_name, "Depth", truncate=False, time_selection=0)
except ValueError as error:
if "not found in HDF file" not in str(error): raise
else:
raise ValueError("A native Depth dataset exists; verify its definition before using invert-depth labels")
wse = HdfResultsMesh.get_mesh_timeseries(plan_hdf, mesh_name, "Water Surface", truncate=False)
topology = HdfMesh.get_mesh_sloped_topology(geom_hdf, mesh_name)
minimum = np.asarray(topology["cell_min_elev"], dtype=np.float32)
values = np.maximum(np.asarray(wse.values, dtype=np.float32) - minimum, 0.0)
values = np.where(np.isfinite(values), values, 0.0).astype(np.float32)
print("Verified raster variable: no native Depth dataset; exporter uses max(Water Surface - Cells Minimum Elevation, 0), in feet.")
return pd.Series(np.max(values, axis=1), index=pd.to_datetime(wse.coords["time"].values))
def build_qualification_summary(runs, *, mode, source_module_sha256, library_version, inspection_time=None):
inspected = inspection_time if inspection_time is not None else datetime.now(timezone.utc)
assert inspected.tzinfo is not None, "inspection_time must be timezone-aware"
inspected = inspected.astimezone(timezone.utc)
return {"mode": mode, "inspection_utc": inspected.isoformat().replace("+00:00", "Z"),
"execution_date": inspected.date().isoformat(),
"source_module_sha256": dict(source_module_sha256), "library_version": library_version,
"clock": "UTC-equivalent naive HEC-RAS model labels", "runs": runs,
"native_runtime_seconds_total": sum(row["native_runtime_seconds"] for row in runs)}
def catalog_dss(dss_file: Path) -> pd.DataFrame:
catalog = RasDss.get_catalog(dss_file)
if catalog.empty:
raise RuntimeError(f"DSS catalog is empty: {Path(dss_file).name}")
return catalog
def validate_dss_precipitation_grids(
dss_file: Path,
catalog: pd.DataFrame,
case: dict,
expected_count: int,
) -> pd.DataFrame:
if len(catalog) != expected_count:
raise ValueError(f"Expected {expected_count} DSS grids, found {len(catalog)}")
rows = []
for pathname in catalog["pathname"].astype(str):
grid = RasDss.read_grid(dss_file, pathname)
parts = grid["metadata"]["pathname_parts"]
assert parts["A"] == "SHG"
assert parts["B"] == case["dss_b"]
assert parts["C"] == "PRECIPITATION"
assert parts["F"] == case["dss_f"]
assert grid["units"].upper() == "MM"
assert grid["data_type"].upper() == "PER-CUM"
assert grid["grid_type"] == "albers"
assert "North_American_1983" in grid["crs"]
assert np.isclose(grid["cell_size"], 2000.0)
values = np.asarray(grid["data"], dtype=float)
finite = values[np.isfinite(values)]
if not len(finite) or np.any(finite < 0):
raise ValueError(f"Invalid DSS grid values for {pathname}")
rows.append(
{
"pathname": pathname,
"start_time": grid["start_time"],
"end_time": grid["end_time"],
"rows": grid["shape"][0],
"columns": grid["shape"][1],
"minimum_mm": float(finite.min()),
"maximum_mm": float(finite.max()),
"mean_mm": float(finite.mean()),
}
)
audit = pd.DataFrame(rows).sort_values("start_time").reset_index(drop=True)
expected_ends = pd.date_range(case["event_start"], periods=expected_count, freq="h")
expected_starts = expected_ends - pd.Timedelta(hours=1)
assert audit["start_time"].tolist() == expected_starts.tolist()
assert audit["end_time"].tolist() == expected_ends.tolist()
assert audit["maximum_mm"].max() > 0
return audit
def get_plan_row(ras, plan_number: str) -> pd.Series:
plan = ras.plan_df[ras.plan_df["plan_number"].astype(str).str.zfill(2) == plan_number.zfill(2)]
if plan.empty:
raise ValueError(f"Plan {plan_number} not found. Available: {ras.plan_df['plan_number'].tolist()}")
return plan.iloc[0]
def get_geometry_hdf(ras, plan_number: str) -> Path:
plan_row = get_plan_row(ras, plan_number)
geom_number = str(plan_row["geometry_number"]).zfill(2)
geom_hdf = ras.project_folder / f"{ras.project_name}.g{geom_number}.hdf"
if not geom_hdf.exists():
raise FileNotFoundError(f"Geometry HDF not found: {geom_hdf}")
return geom_hdf
def get_unsteady_number(ras, plan_number: str) -> str:
plan_row = get_plan_row(ras, plan_number)
return str(plan_row["unsteady_number"]).zfill(2)
def display_preview(frame: pd.DataFrame, rows: int = 3) -> None:
if len(frame) <= rows * 2:
display(frame)
return
print(f"Showing first and last {rows} of {len(frame)} rows")
display(pd.concat([frame.head(rows), frame.tail(rows)]))
@contextmanager
def capture_logger_records(logger_name: str, level: int = logging.WARNING):
logger = logging.getLogger(logger_name)
records = []
class RecordHandler(logging.Handler):
def emit(self, record):
records.append(record)
handler = RecordHandler(level=level)
old_handlers = logger.handlers[:]
old_level = logger.level
old_propagate = logger.propagate
logger.handlers = [handler]
logger.setLevel(level)
logger.propagate = False
try:
yield records
except Exception:
print(f"{logger_name} failed; captured diagnostics:")
for record in records:
print(f"- {record.levelname}: {record.getMessage()}")
raise
finally:
for record in records:
print(f"{logger_name}: {record.levelname}: {record.getMessage()}")
logger.handlers = old_handlers
logger.setLevel(old_level)
logger.propagate = old_propagate
def zero_hyetograph_like(hyetograph):
baseline = hyetograph.copy(deep=True)
baseline["incremental_depth"] = 0.0
baseline["cumulative_depth"] = 0.0
return baseline
def sha256_file(path):
digest = hashlib.sha256()
with Path(path).open("rb") as stream:
for block in iter(lambda: stream.read(1024 * 1024), b""):
digest.update(block)
return digest.hexdigest()
def audit_mrms_source(files, stack):
"""Identify the local GRIB product using NOAA's operational product table."""
import eccodes
rows = []
for index, path in enumerate(files):
with Path(path).open("rb") as stream:
message = eccodes.codes_grib_new_from_file(stream)
try:
metadata = {key: eccodes.codes_get(message, key) for key in (
"discipline", "parameterCategory", "parameterNumber",
"units", "stepType", "validityDate", "validityTime",
)}
finally:
eccodes.codes_release(message)
assert tuple(metadata[key] for key in (
"discipline", "parameterCategory", "parameterNumber"
)) == (209, 6, 37), metadata
assert "MultiSensor_QPE_01H_Pass2" in Path(path).name
raw_time = pd.to_datetime(f"{metadata['validityDate']}{metadata['validityTime']:04d}", format="%Y%m%d%H%M")
filename_time = pd.to_datetime(re.search(r"(\d{8}-\d{6})", Path(path).name).group(1), format="%Y%m%d-%H%M%S")
assert raw_time == filename_time == pd.Timestamp(stack.time.values[index]), "GRIB/filename/stack clock mismatch"
rows.append({"file": str(path), "sha256": sha256_file(path), **metadata,
"resolved_units": "mm", "accumulation_hours": 1})
values = np.asarray(stack.values, dtype=float)
assert values.ndim == 3
finite_mask = np.isfinite(values)
assert np.all(finite_mask == finite_mask[0]), "Changing finite spatial support in MRMS frames"
for row, count in zip(rows, finite_mask.sum(axis=(1, 2))):
row["finite_spatial_cells"] = int(count)
row["total_spatial_cells"] = int(values.shape[1] * values.shape[2])
assert np.all(np.isfinite(values).any(axis=(1, 2))), "All-missing MRMS frame"
assert not np.isinf(values).any(), "Infinite MRMS depth"
assert np.all(values[np.isfinite(values)] >= 0), "Negative/missing-sentinel MRMS depth"
assert stack.attrs["units"] == "mm"
print("GRIB local product 209/6/37: NOAA MultiSensor_QPE_01H_Pass2, 1-hour mm.")
print("ecCodes units may be unknown; units/period resolved from NOAA's product table, not inferred from cadence.")
print("Source: https://www.nssl.noaa.gov/projects/mrms/operational/tables.php")
return pd.DataFrame(rows)
def align_mrms_hyetograph(hyetograph, origin, last_qpe, simulation_end):
"""Keep complete hourly accumulations after origin; append an explicit dry tail."""
source = hyetograph.copy(deep=True)
times = pd.DatetimeIndex(pd.to_datetime(source["time"])).as_unit("ns")
origin, last_qpe, simulation_end = map(pd.Timestamp, (origin, last_qpe, simulation_end))
if times.tz is not None or any(t.tz is not None for t in (origin, last_qpe, simulation_end)):
raise ValueError("Use explicit UTC-equivalent naive model-clock labels")
if times.has_duplicates or not times.is_monotonic_increasing:
raise ValueError("MRMS times must be unique and increasing")
if len(times) < 2 or not np.all(np.diff(times.asi8) == pd.Timedelta(hours=1).value):
raise ValueError("MRMS hourly source has a missing/irregular interval")
if np.any((times.asi8 - origin.value) % pd.Timedelta(hours=1).value):
raise ValueError("MRMS interval straddles the boundary origin")
depths = source["incremental_depth"].to_numpy(dtype=float)
if not np.isfinite(depths).all() or (depths < 0).any():
raise ValueError("Invalid source precipitation")
kept = source.loc[(times > origin) & (times <= last_qpe)].copy()
expected = pd.date_range(origin + pd.Timedelta(hours=1), last_qpe, freq="h")
if not pd.DatetimeIndex(kept["time"]).equals(expected) or not len(expected):
raise ValueError("Source does not completely cover the requested rainfall window")
tail = pd.date_range(last_qpe + pd.Timedelta(hours=1), simulation_end, freq="h")
if simulation_end < last_qpe or (simulation_end - origin) % pd.Timedelta(hours=1):
raise ValueError("Simulation end must be on the hourly lattice after the QPE window")
if len(tail):
kept = pd.concat([kept, pd.DataFrame({"time": tail, "incremental_depth": 0.0})], ignore_index=True)
kept["hour"] = (pd.to_datetime(kept["time"]) - origin) / pd.Timedelta(hours=1)
kept["cumulative_depth"] = kept["incremental_depth"].cumsum()
return kept[["time", "hour", "incremental_depth", "cumulative_depth"]].reset_index(drop=True)
def precipitation_boundary(ras, unsteady_number, name):
boundaries = ras.get_boundary_conditions()
selected = boundaries.loc[
boundaries["unsteady_number"].astype(str).str.zfill(2).eq(unsteady_number)
& boundaries["area_2d"].astype(str).str.casefold().eq(name.casefold())
& boundaries["hydrograph_type"].eq("Precipitation Hydrograph")
]
if len(selected) != 1:
raise ValueError("Expected exactly one selected precipitation boundary")
boundary = selected.iloc[0]
if str(boundary["Use Fixed Start Time"]).strip().lower() != "false":
raise ValueError("This recipe requires the selected boundary to follow simulation start")
return boundary
def readback_hyetograph(ras, unsteady_number, case, intended, condition, output):
boundary = precipitation_boundary(ras, unsteady_number, case["precip_boundary"])
assert boundary["Interval"] == "1HOUR", "Expected 1HOUR boundary interval"
values = np.asarray(boundary["hydrograph_values"], dtype=float)
assert len(values) == len(intended) + 1, "Serialized ordinate count must equal intended plus one"
assert values[0] == 0, "Time-zero anchor must be zero"
times = pd.date_range(case["event_start"], case["sim_end"], freq="h")
assert len(times) == len(values), "Serialized window ordinate count mismatch"
assert times[1:].equals(pd.DatetimeIndex(intended["time"])), "Intended and serialized absolute times differ"
serialized = pd.DataFrame({"time": times, "hour": np.arange(len(values)),
"incremental_depth": values, "cumulative_depth": values.cumsum()})
errors = np.abs(serialized["cumulative_depth"].to_numpy()[1:] - intended["cumulative_depth"].to_numpy())
assert errors.max() <= 0.005 + 1e-12, "Cumulative serialization error exceeds 0.005 in"
serialized.to_csv(output / f"{condition}_serialized.csv", index=False)
print(f"{condition}: readback cumulative rounding error <= {errors.max():.9f} in; total {values.sum():.9f} in")
return serialized
def rainfall_support(geometry_hdf, mesh_name):
polygons = HdfMesh.get_mesh_cell_polygons(geometry_hdf)
physical = polygons.loc[polygons.mesh_name.eq(mesh_name)]
assert len(physical) and physical.cell_id.is_unique and (physical.geometry.area > 0).all()
ids = np.sort(physical.cell_id.to_numpy(dtype=int))
points = HdfMesh.get_mesh_cell_points(geometry_hdf)
all_ids = np.sort(points.loc[points.mesh_name.eq(mesh_name), "cell_id"].to_numpy(dtype=int))
complement = np.setdiff1d(all_ids, ids)
print(f"Detailed polygon API support: {len(ids)} cell polygons; {len(complement)} other cells require topology classification.")
return ids, complement, all_ids
def verify_boundary_topology(plan_hdf, mesh_name, support):
ids, complement, all_ids = support
topology = HdfMesh.get_mesh_sloped_topology(plan_hdf, mesh_name)
assert topology["n_cells"] == len(all_ids), "Topology cell count mismatch"
for cell in complement:
offset, count = topology["cell_face_info"][cell, :2]
assert count == 1, "Complement cell must have one face"
face = int(topology["cell_face_values"][offset, 0])
neighbors = topology["face_cells"][face]
assert cell in neighbors and len(set(neighbors).intersection(ids)) == 1, "Complement face must join one physical cell"
print(f"Geometry-derived rainfall support: {len(ids)} physical cells; {len(complement)} one-face boundary cells.")
def compare_native_cumulative(times, cell_ids, cumulative, rates, expected, support, origin, end):
"""Compare all native values at absolute timestamps; never rebase a result clock."""
times = pd.DatetimeIndex(times).as_unit("ns")
required_times = pd.date_range(origin, end, freq="5min")
if not times.equals(required_times):
raise ValueError("Native output does not have the complete absolute simulation clock")
physical, boundary, all_ids = support
if not np.array_equal(np.sort(cell_ids), all_ids) or len(np.unique(cell_ids)) != len(all_ids):
raise ValueError("Native cell coordinates do not cover the geometry exactly")
if cumulative.shape != rates.shape or cumulative.shape != (len(times), len(cell_ids)):
raise ValueError("Native array dimensions disagree")
if not np.isfinite(cumulative).all() or not np.isfinite(rates).all():
raise ValueError("Nonfinite native precipitation")
if np.any(cumulative < 0) or np.any(rates < 0):
raise ValueError("Negative native precipitation")
source_times = pd.DatetimeIndex(expected["time"]).as_unit("ns")
if source_times[0] != pd.Timestamp(origin) or source_times[-1] != pd.Timestamp(end):
raise ValueError("Serialized boundary does not cover the simulation window")
expected_curve = np.interp(times.asi8, source_times.asi8, expected["cumulative_depth"])
tolerance = 8 * np.finfo(np.float32).eps * max(1.0, float(expected_curve.max()))
columns = np.flatnonzero(np.isin(cell_ids, physical))
errors = np.abs(cumulative[:, columns] - expected_curve[:, None])
max_error = float(errors.max())
if max_error > tolerance:
raise ValueError(f"Absolute-time native cumulative mismatch: {max_error} > {tolerance} in")
outside = np.flatnonzero(np.isin(cell_ids, boundary))
if np.any(cumulative[:, outside] != 0) or np.any(rates[:, outside] != 0):
raise ValueError("Boundary-only cells contain precipitation")
if expected_curve.max() == 0 and (np.any(cumulative != 0) or np.any(rates != 0)):
raise ValueError("No-rain baseline is not zero everywhere")
diagnostics = pd.DataFrame({"time": times, "expected_cumulative_in": expected_curve,
"native_physical_min_in": cumulative[:, columns].min(axis=1),
"native_physical_max_in": cumulative[:, columns].max(axis=1),
"max_abs_error_in": errors.max(axis=1)})
return diagnostics, max_error, tolerance
def validate_hdf_precipitation(plan_hdf, mesh_name, expected_hyetograph, label, case, support, output):
support = reconstruct_physical_support(plan_hdf, mesh_name, support, output)
assert len(support[0]) == case["expected_physical_cells"], "Native physical cell count differs from qualified example"
assert len(support[1]) == case["expected_boundary_cells"], "Native boundary cell count differs from qualified example"
verify_boundary_topology(plan_hdf, mesh_name, support)
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)
assert rate.attrs.get("units") == "in/hr" and cumulative.attrs.get("units") == "in"
assert rate.dims == cumulative.dims and len(rate.dims) == 2
for dim in rate.dims:
assert np.array_equal(rate.coords[dim].values, cumulative.coords[dim].values)
assert rate.dims[0] == "time"
times = pd.DatetimeIndex(cumulative.coords["time"].values)
cell_ids = np.asarray(cumulative.coords[cumulative.dims[1]].values, dtype=int)
diagnostics, error, tolerance = compare_native_cumulative(
times, cell_ids, np.asarray(cumulative.values, dtype=float), np.asarray(rate.values, dtype=float),
expected_hyetograph, support, case["event_start"], case["sim_end"])
diagnostics.to_csv(output / f"{label}_native_cumulative.csv", index=False)
physical, boundary, _ = support
return pd.DataFrame([{"condition": label, "timestamps": len(times), "physical_cells_checked": len(physical),
"boundary_cells_zero": len(boundary), "absolute_start_utc": times[0], "absolute_end_utc": times[-1],
"expected_total_in": expected_hyetograph["cumulative_depth"].iloc[-1],
"hdf_final_total_in": float(cumulative.values[-1].max()),
"max_all_time_cell_error_in": error, "native_tolerance_in": tolerance,
"hdf_peak_rate_in_hr": float(rate.values.max()), "absolute_time_qualified": True}])
def compute_with_receipt(ras, case, variables, condition, output):
plan = case["plan_number"]
row = get_plan_row(ras, plan)
inputs = [Path(row[key]) for key in ("full_path", "Geom Path", "Flow Path")]
before = {str(path): sha256_file(path) for path in inputs}
snapshot = output / f"{condition}_inputs"
snapshot.mkdir(exist_ok=False)
for path in inputs:
shutil.copy2(path, snapshot / path.name)
started = time.time()
print(f"Launching {case['project']} {condition}: 2 cores, deadline {case['max_runtime']} s", flush=True)
result = RasCmdr.compute_plan(plan, ras_object=ras, force_rerun=True, verify=True,
num_cores=2, max_runtime=case["max_runtime"], use_optimal_hdf_settings=True,
hdf_output_variables=variables)
inventory = RasCmdr.inspect_plan_processes(plan, ras_object=ras)
hdf = ras.project_folder / f"{ras.project_name}.p{plan}.hdf"
after_snapshot = output / f"{condition}_inputs_after"
after_snapshot.mkdir(exist_ok=True)
for path in inputs:
shutil.copy2(path, after_snapshot / path.name)
receipt = {"condition": condition, "started_epoch": started, "finished_epoch": time.time(),
"max_runtime_seconds": case["max_runtime"], "num_cores": 2,
"success": bool(result), "completion_verified": result.completion_verified,
"execution_details": result.execution_details,
"results_df_row": None if result.results_df_row is None else result.results_df_row.to_dict(),
"process_inventory": inventory.to_dict(), "input_sha256_before": before,
"input_sha256_after": {str(path): sha256_file(path) for path in inputs},
"hdf_sha256": sha256_file(hdf) if hdf.exists() else None,
"source_sha256": SOURCE_SHA256, "python": sys.version, "executable": sys.executable}
(output / f"{condition}_compute_receipt.json").write_text(json.dumps(receipt, default=str, indent=2))
print(f"Terminal {condition}: {result}; elapsed {receipt['finished_epoch'] - started:.2f}s; process matches={len(inventory.matched)}", flush=True)
assert bool(result) and result.completion_verified
assert inventory.complete and not inventory.matched and hdf.stat().st_mtime >= started
return result
def hdf_hydraulic_depth_maxima(plan_hdf: Path, mesh_name: str) -> pd.Series:
depth_da = HdfResultsMesh.get_mesh_timeseries(
plan_hdf, mesh_name, "Cell Hydraulic Depth", truncate=False
)
values = np.asarray(depth_da.values, dtype=float)
maxima = np.nanmax(values, axis=1)
return pd.Series(maxima, index=pd.to_datetime(depth_da.coords["time"].values))
def get_terrain_inventory(ras) -> pd.DataFrame:
import rasterio
terrain_dir = ras.project_folder / "Terrain"
if not terrain_dir.exists() or "terrain_hdf_path" not in ras.rasmap_df.columns:
return pd.DataFrame()
active_hdfs = []
for value in ras.rasmap_df["terrain_hdf_path"].dropna():
values = value if isinstance(value, (list, tuple, set)) else [value]
active_hdfs.extend(Path(path) for path in values if path)
candidates: dict[Path, Path] = {}
for terrain_hdf in active_hdfs:
for terrain_path in terrain_dir.glob(f"{terrain_hdf.stem}*.tif"):
candidates[terrain_path.resolve()] = terrain_hdf
rows = []
for terrain_path, terrain_hdf in sorted(candidates.items()):
with rasterio.open(terrain_path) as src:
x_resolution = abs(float(src.res[0]))
y_resolution = abs(float(src.res[1]))
valid_cells = sum(
int(np.count_nonzero(src.read_masks(1, window=window)))
for _, window in src.block_windows(1)
)
rows.append(
{
"terrain": terrain_path.name,
"terrain_hdf": terrain_hdf.name,
"path": terrain_path,
"x_resolution": x_resolution,
"y_resolution": y_resolution,
"linear_units": src.crs.linear_units if src.crs else None,
"valid_cells": valid_cells,
"valid_coverage_area": valid_cells * x_resolution * y_resolution,
}
)
inventory = pd.DataFrame(rows)
if inventory.empty:
return inventory
inventory = inventory.sort_values("valid_coverage_area", ascending=False).reset_index(drop=True)
inventory["nominal_resolution"] = inventory["x_resolution"].round(2)
precedence = inventory.sort_values(
["nominal_resolution", "valid_coverage_area"],
ascending=[False, False],
).index
inventory["mosaic_order"] = 0
inventory.loc[precedence, "mosaic_order"] = np.arange(1, len(inventory) + 1)
return inventory
def get_pump_stations(geom_hdf: Path, mesh_areas: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
try:
pumps = HdfPump.get_pump_stations(geom_hdf)
except Exception:
return gpd.GeoDataFrame({"Name": []}, geometry=[], crs=mesh_areas.crs)
if pumps.empty and mesh_areas.crs is not None:
pumps = pumps.set_crs(mesh_areas.crs)
return pumps
def make_clip_shp(bounds: tuple[float, float, float, float], output_path: Path, label: str) -> Path:
output_path.parent.mkdir(parents=True, exist_ok=True)
gdf = gpd.GeoDataFrame({"name": [label]}, geometry=[box(*bounds)], crs="EPSG:4326")
gdf.to_file(output_path)
return output_path
def parse_ras_timestamp(timestamp: str) -> datetime:
return datetime.strptime(timestamp, "%d%b%Y %H:%M:%S")
def select_animation_timestamps(
timestamps: list[str],
start: datetime,
end: datetime,
) -> list[str]:
selected = [
timestamp
for timestamp in timestamps
if start <= parse_ras_timestamp(timestamp) <= end
]
if len(selected) < 2:
raise ValueError(f"Animation window {start} to {end} selected {len(selected)} timestamp(s)")
selected_start = parse_ras_timestamp(selected[0])
selected_end = parse_ras_timestamp(selected[-1])
if selected_start != start or selected_end != end:
raise ValueError(
f"Result coverage {selected_start} to {selected_end} does not match "
f"the requested event window {start} to {end}"
)
deltas = [
parse_ras_timestamp(selected[idx + 1]) - parse_ras_timestamp(selected[idx])
for idx in range(len(selected) - 1)
]
if any(delta != timedelta(minutes=5) for delta in deltas):
raise ValueError("Selected animation timestamps are not consecutive 5-minute frames")
return selected
def get_depth_tiles(
depth_results: dict[str, dict[str, list[Path]]],
timestamp: str,
terrain_inventory: pd.DataFrame | None = None,
) -> list[Path]:
priority = {}
if terrain_inventory is not None and not terrain_inventory.empty:
priority = {
str(row.terrain).casefold(): int(row.mosaic_order)
for row in terrain_inventory.itertuples()
}
def tile_sort_key(path: Path) -> tuple[int, str]:
path_name = path.name.casefold()
matches = [order for terrain, order in priority.items() if terrain in path_name]
return (matches[0] if matches else 0, path_name)
tiles = sorted(
(Path(path) for path in depth_results.get(timestamp, {}).get("depth", [])),
key=tile_sort_key,
)
return [path for path in tiles if path.exists()]
def audit_mapper_depth_frames(
depth_results: dict[str, dict[str, list[Path]]],
timestamps: list[str],
terrain_inventory: pd.DataFrame,
hdf_depth_max: pd.Series,
) -> pd.DataFrame:
import rasterio
dominant_name = str(terrain_inventory.iloc[0]["terrain"]).casefold()
rows = []
for timestamp in timestamps:
tiles = get_depth_tiles(depth_results, timestamp, terrain_inventory)
dominant_present = any(dominant_name in path.name.casefold() for path in tiles)
readable = bool(tiles)
georeferenced = bool(tiles)
all_have_data = bool(tiles)
mapper_max_depth = 0.0
for path in tiles:
try:
with rasterio.open(path) as src:
georeferenced = georeferenced and src.crs is not None and not src.transform.is_identity
data = src.read(1, masked=True)
has_data = bool(data.count())
all_have_data = all_have_data and has_data
if has_data:
mapper_max_depth = max(mapper_max_depth, float(data.max()))
except Exception:
readable = False
result_time = pd.Timestamp(parse_ras_timestamp(timestamp))
hdf_max_depth = float(hdf_depth_max.get(result_time, np.nan))
wetness_matches_hdf = not (hdf_max_depth > 0.01 and mapper_max_depth <= 0.01)
valid = (
bool(tiles)
and dominant_present
and readable
and georeferenced
and all_have_data
and wetness_matches_hdf
)
reasons = []
if not tiles:
reasons.append("no depth raster")
if tiles and not dominant_present:
reasons.append("dominant terrain missing")
if not readable:
reasons.append("unreadable raster")
if tiles and not georeferenced:
reasons.append("missing georeferencing")
if tiles and not all_have_data:
reasons.append("all-nodata raster")
if not wetness_matches_hdf:
reasons.append("blank while HDF is wet")
rows.append(
{
"timestamp": timestamp,
"terrain_tiles": len(tiles),
"dominant_present": dominant_present,
"mapper_max_depth_ft": mapper_max_depth,
"hdf_max_depth_ft": hdf_max_depth,
"valid": valid,
"reason": "; ".join(reasons),
}
)
return pd.DataFrame(rows)
def collect_depth_raster_frames(
depth_results: dict[str, dict[str, list[Path]]],
timestamps: list[str],
terrain_inventory: pd.DataFrame,
) -> tuple[list[list[Path]], pd.DataFrame]:
frames: list[list[Path]] = []
audit_rows: list[dict] = []
missing: list[str] = []
for timestamp in timestamps:
tiles = get_depth_tiles(depth_results, timestamp, terrain_inventory)
if not tiles:
missing.append(timestamp)
continue
frames.append(tiles)
audit_rows.append(
{
"timestamp": parse_ras_timestamp(timestamp),
"terrain_tiles": len(tiles),
"total_raster_mb": sum(path.stat().st_size for path in tiles) / 1_000_000,
"rasters": "; ".join(path.name for path in tiles),
}
)
if missing:
preview = ", ".join(missing[:5])
raise RuntimeError(
f"Stored-map depth rasters are missing for {len(missing)} timestep(s): {preview}"
)
return frames, pd.DataFrame(audit_rows)
def summarize_depth_stack(depth_stack) -> pd.DataFrame:
values = np.asarray(depth_stack.values, dtype=np.float32)
x = np.asarray(depth_stack.coords["x"].values, dtype=float)
y = np.asarray(depth_stack.coords["y"].values, dtype=float)
dx = float(np.median(np.abs(np.diff(x)))) if len(x) > 1 else MIN_CONSOLIDATION_CELL_SIZE_FT
dy = float(np.median(np.abs(np.diff(y)))) if len(y) > 1 else MIN_CONSOLIDATION_CELL_SIZE_FT
cell_area_acres = dx * dy / 43560.0
rows = []
for frame_index, timestamp in enumerate(pd.to_datetime(depth_stack.coords["time"].values)):
data = values[frame_index]
wet = np.isfinite(data) & (data > 0)
rows.append(
{
"timestamp": timestamp,
"wet_cells": int(np.count_nonzero(wet)),
"wet_area_acres": float(np.count_nonzero(wet) * cell_area_acres),
"raster_invert_max_ft": float(np.nanmax(data)) if np.isfinite(data).any() else 0.0,
"raster_invert_mean_wet_ft": float(np.mean(data[wet])) if wet.any() else 0.0,
}
)
return pd.DataFrame(rows)
def max_depth_grid_from_stack(depth_stack):
values = np.asarray(depth_stack.values, dtype=np.float32)
finite = np.isfinite(values)
max_depth = np.max(np.where(finite, values, -np.inf), axis=0)
max_depth[~finite.any(axis=0)] = np.nan
max_extent, _ = PrecipMrms._data_extent(depth_stack)
max_crs = depth_stack.attrs.get("crs")
return max_depth, max_extent, max_crs
def plot_accumulation_map(
case: dict,
mrms_stack,
hyetograph: pd.DataFrame,
mesh_areas: gpd.GeoDataFrame,
pump_stations: gpd.GeoDataFrame,
):
accum_in = (mrms_stack.sum(dim="time") / 25.4).rename("MRMS accumulation")
extent, origin = PrecipMrms._data_extent(mrms_stack)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), constrained_layout=True)
axes[0].bar(hyetograph["time"], hyetograph["incremental_depth"], width=-1/24, align="edge", color="#2b8cbe")
axes[0].plot(hyetograph["time"], hyetograph["cumulative_depth"], color="#08306b", marker="o")
axes[0].set_ylabel("Depth (in)")
axes[0].set_title(f"{case['display_name']} MRMS hyetograph")
axes[0].tick_params(axis="x", rotation=30)
im = axes[1].imshow(
accum_in.values,
extent=extent,
origin=origin,
cmap="turbo",
vmin=0,
alpha=0.78,
zorder=2,
)
axes[1].set_title("MRMS accumulation with model context")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")
from matplotlib.ticker import MaxNLocator
axes[1].xaxis.set_major_locator(MaxNLocator(nbins=4))
PrecipMrms._plot_spatial_overlays(
axes[1],
data_crs="EPSG:4326",
mesh_boundary=mesh_areas,
pump_stations=pump_stations,
add_basemap=False,
)
fig.colorbar(im, ax=axes[1], shrink=0.8, label="Accumulation (in)")
for axis in fig.axes:
axis.grid(True, alpha=0.2)
plt.show()
def plot_depth_response(
case: dict,
depth_summary: pd.DataFrame,
depth_stack,
mesh_areas: gpd.GeoDataFrame,
pump_stations: gpd.GeoDataFrame,
):
max_depth, max_extent, max_crs = max_depth_grid_from_stack(depth_stack)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), constrained_layout=True)
axes[0].plot(
depth_summary["timestamp"],
depth_summary["raster_invert_max_ft"],
color="#08519c",
marker="o",
)
axes[0].set_title(f"{case['project']}: selected-frame invert depth")
axes[0].set_ylabel("Maximum raster invert depth (ft)")
axes[0].tick_params(axis="x", rotation=30)
vmax = max(float(np.nanmax(max_depth)), 0.1)
im = axes[1].imshow(
max_depth,
extent=max_extent,
origin="upper",
cmap="Blues",
vmin=0,
vmax=vmax,
alpha=0.78,
zorder=2,
)
axes[1].set_title("Maximum selected-frame invert depth")
axes[1].set_xlabel("Easting (ft)")
axes[1].set_ylabel("Northing (ft)")
axes[1].set_aspect("equal", adjustable="box")
PrecipMrms._plot_spatial_overlays(
axes[1],
data_crs=max_crs,
mesh_boundary=mesh_areas,
pump_stations=pump_stations,
add_basemap=False,
)
fig.colorbar(im, ax=axes[1], shrink=0.8, label="Invert depth (ft)")
for axis in fig.axes:
axis.grid(True, alpha=0.2)
plt.show()
def fresh_artifact_path(path: Path) -> Path:
path = Path(path)
path.parent.mkdir(parents=True, exist_ok=True)
path.unlink(missing_ok=True)
return path
def write_consolidated_depth_rasters(
depth_stack,
output_dir: Path,
) -> tuple[list[Path], pd.DataFrame]:
import rasterio
from rasterio.transform import from_origin
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
x = np.asarray(depth_stack.coords["x"].values, dtype=float)
y = np.asarray(depth_stack.coords["y"].values, dtype=float)
if len(x) < 2 or len(y) < 2 or not np.all(np.diff(x) > 0) or not np.all(np.diff(y) < 0):
raise ValueError("Consolidated raster coordinates must be x-ascending and y-descending")
dx = float(np.median(np.diff(x)))
dy = float(np.median(np.abs(np.diff(y))))
transform = from_origin(x[0] - dx / 2, y[0] + dy / 2, dx, dy)
nodata = np.float32(-9999.0)
output_paths = []
audit_rows = []
expected_crs = rasterio.crs.CRS.from_user_input(depth_stack.attrs.get("crs"))
for frame_index, time_value in enumerate(depth_stack.coords["time"].values):
timestamp = pd.Timestamp(time_value)
output_path = output_dir / f"depth_{timestamp:%Y%m%d_%H%M%S}.tif"
values = np.asarray(depth_stack.isel(time=frame_index).values, dtype=np.float32)
write_values = np.where(np.isfinite(values), values, nodata).astype(np.float32)
with rasterio.open(
output_path,
"w",
driver="GTiff",
height=write_values.shape[0],
width=write_values.shape[1],
count=1,
dtype="float32",
crs=depth_stack.attrs.get("crs"),
transform=transform,
nodata=float(nodata),
compress="deflate",
predictor=3,
BIGTIFF="IF_SAFER",
) as dst:
dst.write(write_values, 1)
dst.update_tags(timestamp=timestamp.isoformat(), units=depth_stack.attrs.get("units", ""))
with rasterio.open(output_path) as check:
assert check.crs == expected_crs
assert check.transform.almost_equals(transform)
assert np.allclose(check.res, (dx, dy))
assert (check.height, check.width) == write_values.shape
assert np.isclose(check.nodata, float(nodata))
assert check.tags().get("timestamp") == timestamp.isoformat()
audit_rows.append(
{
"timestamp": timestamp,
"raster": output_path.name,
"width": check.width,
"height": check.height,
"cell_size_ft": check.res[0],
"size_mb": output_path.stat().st_size / 1_000_000,
}
)
output_paths.append(output_path)
return output_paths, pd.DataFrame(audit_rows)
def select_reference_cell_by_depth_delta(
event_plan_hdf: Path,
baseline_plan_hdf: Path,
mesh_name: str,
) -> tuple[int, pd.DataFrame]:
"""Select the cell with the largest concurrent event-minus-baseline depth."""
event_da = HdfResultsMesh.get_mesh_timeseries(
event_plan_hdf, mesh_name, "Cell Hydraulic Depth", truncate=False
)
baseline_da = HdfResultsMesh.get_mesh_timeseries(
baseline_plan_hdf, mesh_name, "Cell Hydraulic Depth", truncate=False
)
event_cells = np.asarray(event_da.coords["cell_id"].values, dtype=int)
baseline_cells = np.asarray(baseline_da.coords["cell_id"].values, dtype=int)
event_times = pd.to_datetime(event_da.coords["time"].values)
baseline_times = pd.to_datetime(baseline_da.coords["time"].values)
if not np.array_equal(event_cells, baseline_cells):
raise ValueError(f"Baseline and event cell IDs differ for mesh {mesh_name}")
if not np.array_equal(event_times, baseline_times):
raise ValueError(f"Baseline and event output times differ for mesh {mesh_name}")
event_depth = np.asarray(event_da.values, dtype=float)
baseline_depth = np.asarray(baseline_da.values, dtype=float)
depth_delta = event_depth - baseline_depth
finite = np.isfinite(depth_delta)
peak_delta_by_cell = np.max(np.where(finite, depth_delta, -np.inf), axis=0)
peak_delta_by_cell[~finite.any(axis=0)] = np.nan
if not np.isfinite(peak_delta_by_cell).any():
raise RuntimeError(f"No finite event-minus-baseline depths found for mesh {mesh_name}")
reference_idx = int(np.nanargmax(peak_delta_by_cell))
peak_time_idx = int(np.nanargmax(depth_delta[:, reference_idx]))
reference_cell_id = int(event_cells[reference_idx])
selection = pd.DataFrame(
[
{
"mesh_name": mesh_name,
"reference_cell_id": reference_cell_id,
"selection_metric": "maximum concurrent event-minus-baseline depth",
"peak_response_time": event_times[peak_time_idx],
"baseline_hydraulic_depth_ft": float(baseline_depth[peak_time_idx, reference_idx]),
"event_hydraulic_depth_ft": float(event_depth[peak_time_idx, reference_idx]),
"event_minus_baseline_hydraulic_depth_ft": float(
depth_delta[peak_time_idx, reference_idx]
),
}
]
)
return reference_cell_id, selection
def build_reference_response(
plan_hdf: Path,
geom_hdf: Path,
mesh_name: str,
reference_cell_id: int | None = None,
):
depth_da = HdfResultsMesh.get_mesh_timeseries(plan_hdf, mesh_name, "Cell Hydraulic Depth", truncate=False)
wse_da = HdfResultsMesh.get_mesh_timeseries(plan_hdf, mesh_name, "Water Surface", truncate=False)
precip_da = HdfResultsMesh.get_mesh_timeseries(plan_hdf, mesh_name, "Cell Precipitation Rate", truncate=False)
cumulative_precip_da = HdfResultsMesh.get_mesh_timeseries(
plan_hdf, mesh_name, "Cell Cumulative Precipitation Depth", truncate=False
)
depth_values = np.asarray(depth_da.values, dtype=float)
peak_by_cell = np.nanmax(depth_values, axis=0)
if not np.isfinite(peak_by_cell).any():
raise RuntimeError(f"No finite depth values found for mesh {mesh_name}")
if reference_cell_id is None:
reference_idx = int(np.nanargmax(peak_by_cell))
reference_cell_id = int(depth_da.coords["cell_id"].values[reference_idx])
else:
cell_ids = np.asarray(depth_da.coords["cell_id"].values, dtype=int)
matches = np.flatnonzero(cell_ids == int(reference_cell_id))
if not len(matches):
raise ValueError(f"Cell {reference_cell_id} not found in mesh {mesh_name}")
reference_idx = int(matches[0])
reference_response = pd.DataFrame(
{
"time": pd.to_datetime(depth_da.coords["time"].values),
"hydraulic_depth_ft": depth_values[:, reference_idx],
"wse_ft": np.asarray(wse_da.sel(cell_id=reference_cell_id).values, dtype=float),
"precip_in_hr": np.asarray(precip_da.sel(cell_id=reference_cell_id).values, dtype=float),
"cumulative_precip_in": np.asarray(
cumulative_precip_da.sel(cell_id=reference_cell_id).values, dtype=float
),
}
)
reference_response.attrs["mesh_name"] = mesh_name
reference_response.attrs["cell_id"] = reference_cell_id
cell_points = HdfMesh.get_mesh_cell_points(geom_hdf)
reference_point = cell_points[
(cell_points["mesh_name"] == mesh_name) & (cell_points["cell_id"] == reference_cell_id)
].copy()
return reference_response, reference_point, reference_cell_id
def plot_reference_response(
case: dict,
hyetograph: pd.DataFrame,
event_response: pd.DataFrame,
baseline_response: pd.DataFrame,
reference_point: gpd.GeoDataFrame,
depth_stack,
mesh_areas: gpd.GeoDataFrame,
pump_stations: gpd.GeoDataFrame,
):
max_depth, max_extent, max_crs = max_depth_grid_from_stack(depth_stack)
import matplotlib.dates as mdates
cell_id = event_response.attrs["cell_id"]
mesh_name = event_response.attrs["mesh_name"]
fig = plt.figure(figsize=(13, 7), constrained_layout=True)
grid = fig.add_gridspec(2, 2, width_ratios=[1.25, 1.0])
ax_precip = fig.add_subplot(grid[0, 0])
ax_hydro = fig.add_subplot(grid[1, 0], sharex=ax_precip)
ax_map = fig.add_subplot(grid[:, 1])
ax_precip.bar(
hyetograph["time"],
hyetograph["incremental_depth"],
width=-1/24, align="edge",
color="#2b8cbe",
label="Spatial mean MRMS",
)
ax_precip.set_ylabel("Incremental depth (in)")
ax_precip.set_title(f"{case['display_name']} rainfall input")
ax_precip.legend(loc="upper right")
ax_hydro.plot(
baseline_response["time"], baseline_response["hydraulic_depth_ft"],
color="#636363", linestyle="--", label="No-rain hydraulic depth"
)
ax_hydro.plot(
event_response["time"], event_response["hydraulic_depth_ft"],
color="#08519c", label="MRMS-event hydraulic depth"
)
ax_hydro.set_ylabel("Hydraulic depth (ft)")
ax_hydro_wse = ax_hydro.twinx()
ax_hydro_wse.plot(
baseline_response["time"], baseline_response["wse_ft"],
color="#969696", linestyle=":", label="No-rain WSE"
)
ax_hydro_wse.plot(
event_response["time"], event_response["wse_ft"],
color="#238b45", label="MRMS-event WSE"
)
ax_hydro_wse.set_ylabel("WSE (ft)")
ax_hydro.set_title(f"Maximum rainfall-response cell {cell_id} hydrograph ({mesh_name})")
ax_hydro.xaxis.set_major_locator(mdates.AutoDateLocator(minticks=4, maxticks=7))
ax_hydro.xaxis.set_major_formatter(mdates.ConciseDateFormatter(ax_hydro.xaxis.get_major_locator()))
lines = ax_hydro.get_lines() + ax_hydro_wse.get_lines()
ax_hydro.legend(lines, [line.get_label() for line in lines], loc="upper left", fontsize=8)
vmax = max(float(np.nanmax(max_depth)), 0.1)
im = ax_map.imshow(
max_depth,
extent=max_extent,
origin="upper",
cmap="Blues",
vmin=0,
vmax=vmax,
alpha=0.78,
zorder=2,
)
ax_map.set_title("Reference cell on maximum invert-depth map")
ax_map.set_xlabel("Easting (ft)")
ax_map.set_ylabel("Northing (ft)")
ax_map.set_aspect("equal", adjustable="box")
PrecipMrms._plot_spatial_overlays(
ax_map,
data_crs=max_crs,
mesh_boundary=mesh_areas,
pump_stations=pump_stations,
add_basemap=False,
)
if not reference_point.empty:
plot_point = reference_point.to_crs(max_crs) if max_crs and reference_point.crs else reference_point
plot_point.plot(ax=ax_map, marker="*", color="#ffea00", edgecolor="black", markersize=160, zorder=8)
fig.colorbar(im, ax=ax_map, shrink=0.8, label="Invert depth (ft)")
for axis in fig.axes:
axis.grid(True, alpha=0.2)
plt.show()
def get_pump_operation(plan_hdf: Path, pump_stations: gpd.GeoDataFrame) -> dict[str, pd.DataFrame]:
operations: dict[str, pd.DataFrame] = {}
if pump_stations.empty or "Name" not in pump_stations.columns:
return operations
for station_name in pump_stations["Name"].dropna().astype(str):
try:
pump_da = HdfPump.get_pump_station_timeseries(plan_hdf, pump_station=station_name)
except Exception as exc:
print(f"Pump operation unavailable for {station_name}: {exc}")
continue
frame = pd.DataFrame(
np.asarray(pump_da.values, dtype=float),
columns=[str(value) for value in pump_da.coords["variable"].values],
)
frame.insert(0, "time", pd.to_datetime(pump_da.coords["time"].values))
frame.attrs["unit_by_variable"] = pump_da.attrs.get("unit_by_variable", {})
operations[station_name] = frame
return operations
def plot_pump_operation(
case: dict,
event_operation: dict[str, pd.DataFrame],
baseline_operation: dict[str, pd.DataFrame],
) -> pd.DataFrame:
if not event_operation:
print(f"No pump operation time series available for {case['display_name']}")
return pd.DataFrame()
import matplotlib.dates as mdates
fig, axes = plt.subplots(
len(event_operation),
2,
figsize=(13, 3.5 * len(event_operation)),
squeeze=False,
constrained_layout=True,
)
summary_rows = []
for row_idx, (station_name, event_frame) in enumerate(event_operation.items()):
baseline_frame = baseline_operation.get(station_name, pd.DataFrame())
flow_col = "Flow" if "Flow" in event_frame else next(
(column for column in event_frame.columns if column.endswith(" Flow")), None
)
on_col = "Pumps on" if "Pumps on" in event_frame else next(
(column for column in event_frame.columns if column.endswith("Pumps on")), None
)
ax_flow = axes[row_idx, 0]
if flow_col is not None and not baseline_frame.empty and flow_col in baseline_frame:
ax_flow.plot(
baseline_frame["time"], baseline_frame[flow_col],
color="#636363", linestyle="--", label="No-rain"
)
if flow_col is not None:
ax_flow.plot(
event_frame["time"], event_frame[flow_col],
color="#08519c", label="MRMS event"
)
ax_flow.set_title(f"{station_name} total pump flow")
ax_flow.set_ylabel("Flow (cfs)")
ax_flow.legend(loc="upper left", fontsize=8)
ax_on = axes[row_idx, 1]
if on_col is not None and not baseline_frame.empty and on_col in baseline_frame:
ax_on.step(
baseline_frame["time"], baseline_frame[on_col], where="post",
color="#636363", linestyle="--", label="No-rain"
)
if on_col is not None:
ax_on.step(
event_frame["time"], event_frame[on_col], where="post",
color="#238b45", label="MRMS event"
)
ax_on.set_title(f"{station_name} pump activation")
ax_on.set_ylabel("Pumps on")
ax_on.legend(loc="upper left", fontsize=8)
for axis in (ax_flow, ax_on):
locator = mdates.AutoDateLocator(minticks=3, maxticks=5)
axis.xaxis.set_major_locator(locator)
axis.xaxis.set_major_formatter(mdates.ConciseDateFormatter(locator))
event_peak_flow = float(event_frame[flow_col].max()) if flow_col else np.nan
baseline_peak_flow = (
float(baseline_frame[flow_col].max())
if flow_col and not baseline_frame.empty and flow_col in baseline_frame
else np.nan
)
event_on = event_frame[on_col] if on_col else pd.Series(0, index=event_frame.index)
active = event_frame.loc[event_on > 0, "time"]
summary_rows.append(
{
"station": station_name,
"baseline_peak_flow_cfs": baseline_peak_flow,
"event_peak_flow_cfs": event_peak_flow,
"peak_flow_delta_cfs": event_peak_flow - baseline_peak_flow,
"event_max_pumps_on": float(event_on.max()) if len(event_on) else 0.0,
"first_on": active.iloc[0] if not active.empty else pd.NaT,
"last_on": active.iloc[-1] if not active.empty else pd.NaT,
}
)
fig.suptitle(f"No-rain baseline versus MRMS-event pump operation: {case['display_name']}")
for axis in fig.axes:
axis.grid(True, alpha=0.2)
plt.show()
return pd.DataFrame(summary_rows)
def display_case_video_artifacts(case: dict, videos: list[tuple[str, Path]]) -> None:
rows = []
for label, path in videos:
path = Path(path)
assert path.exists() and path.stat().st_size > 0
rows.append(
{
"animation": label,
"file": path.name,
"size_mb": path.stat().st_size / 1_000_000,
}
)
print(f"MP4 animations generated for {case['display_name']}:")
display(pd.DataFrame(rows).round({"size_mb": 2}))
def audit_runtime_messages(ras, plan_number: str, plan_hdf: Path, label: str) -> dict:
messages = HdfResultsPlan.get_compute_messages_hdf_only(plan_hdf)
lower = messages.lower()
bco_path = ras.project_folder / f"{ras.project_name}.bco{plan_number}"
bco_text = BcoMonitor(project_path=ras.project_folder, plan_number=plan_number, project_name=ras.project_name).get_final_messages() or ""
(ras.project_folder / "P" / f"{label.replace(' ', '_')}_runtime_messages.json").write_text(
json.dumps({"hdf_compute_messages": messages, "bco_messages": bco_text}, indent=2))
return summarize_runtime_messages(messages, bco_text, label)
def validate_source_hyetograph(source_hyetograph, stack):
assert stack.attrs.get("units") == "mm", "Source GRIB stack must use millimeters"
assert pd.DatetimeIndex(source_hyetograph["time"]).equals(pd.DatetimeIndex(stack.time.values)), "Source hyetograph timestamps differ from GRIB stack"
expected_depths = np.nanmean(np.asarray(stack.values, dtype=float), axis=(1, 2)) / 25.4
dtype = np.dtype(np.result_type(stack.dtype, np.float32))
factor = 4 if dtype.itemsize <= 4 else 8
# Unit conversion occurs per cell before the API takes its float64 mean.
# Bound source floating-point roundoff separately from writer quantization.
bound = factor * np.finfo(dtype).eps * np.abs(expected_depths)
error = np.abs(source_hyetograph["incremental_depth"].to_numpy() - expected_depths)
assert np.all(error <= bound), "Source hyetograph depths differ from unrounded spatial mean"
print(f"Source precision: dtype={dtype}; max difference={error.max():.12g} in; max {factor}-epsilon bound={bound.max():.12g} in; exact zero required for dry frames.")
return pd.DataFrame({"time": source_hyetograph["time"], "reference_float64_mean_in": expected_depths,
"api_mean_in": source_hyetograph["incremental_depth"], "abs_error_in": error, "source_roundoff_bound_in": bound})
def summarize_runtime_messages(messages, bco_text, label):
pattern = re.compile(r"\b(?:warning|fatal|failed|missing)\b|\berror(?::|\s+(?:reading|opening|processing|computing|occurred|encountered))", re.IGNORECASE)
lines = [f"{source}: {line.strip()}" for source, text in (("HDF", messages), ("BCO", bco_text))
for line in text.splitlines() if pattern.search(line)]
iteration_lines, volume_lines = [], []
in_iterations = False
for line in messages.splitlines():
if "Maximum iteration location" in line or "Pipe Network Iter" in line:
in_iterations = True
if "Overall Volume Accounting" in line:
in_iterations = False
volume_lines.append(line)
if in_iterations and line.strip(): iteration_lines.append(line)
print(f"{label}: warning/error-pattern lines from HDF and BCO: {len(lines)}; iteration-table lines: {len(iteration_lines)}")
for line in lines: print(line)
print(f"{label}: complete retained HDF compute messages (numerical diagnostics require engineering review):")
print(messages)
return {"condition": label, "complete_process": "complete process" in messages.lower(),
"bco_exists": bool(bco_text), "diagnostic_line_count": len(lines), "diagnostic_lines": lines,
"iteration_diagnostic_lines": iteration_lines, "volume_accounting_lines": volume_lines,
"hydraulic_acceptance": "not_established"}
def preserve_after_input_snapshots(receipt, evidence, condition):
"""Verify existing after-images first; resolve all recoveries before copying."""
destination = evidence / f"{condition}_inputs_after"
manifest_path = destination / "capture_manifest.json"
if manifest_path.exists():
manifest = json.loads(manifest_path.read_text())
recorded = {Path(item["snapshot"]).name: item for item in manifest["files"]}
assert set(recorded) == {Path(value).name for value in receipt["input_sha256_after"]}
for original, expected_hash in receipt["input_sha256_after"].items():
item = recorded[Path(original).name]
assert item["sha256"] == expected_hash
assert sha256_file(Path(item["snapshot"])) == expected_hash, "Postcompute snapshot hash mismatch"
return
captures = []
for original, expected_hash in receipt["input_sha256_after"].items():
target = destination / Path(original).name
if target.exists():
assert sha256_file(target) == expected_hash, "Postcompute snapshot hash mismatch"
source = target
else:
candidates = [evidence / f"{condition}_inputs" / Path(original).name, Path(original)]
source = next((candidate for candidate in candidates if candidate.exists() and sha256_file(candidate) == expected_hash), None)
if source is None:
raise ValueError(f"Cannot recover exact postcompute input bytes: {original}")
captures.append({"source": str(source), "snapshot": str(target), "sha256": expected_hash})
# Missing/mismatched sources above leave all existing evidence untouched.
destination.mkdir(exist_ok=True)
for item in captures:
target = Path(item["snapshot"])
if not target.exists(): shutil.copy2(Path(item["source"]), target)
assert sha256_file(target) == item["sha256"], "Postcompute snapshot hash mismatch"
manifest = {"captured_epoch": time.time(), "capture_scope": "Existing after-images verified; missing images recovered by exact original postcompute hash", "files": captures}
manifest_path.write_text(json.dumps(manifest, indent=2))
def inspect_retained_case(case):
print(f"INSPECT_RETAINED: {case['project']}; no native solver launched by this cell.")
project_path = RUN_ROOT / "m" / case["project"]
ras = init_ras_project(project_path, ras_version=case["ras_version"], ras_object=RasPrj(), load_results_summary=False, hide_intro=True)
assert RasPrj.get_project_units(RasPrj.find_ras_prj(ras.project_folder)) == "ft", "This example requires a US-customary project"
plan_number = case["plan_number"]
unsteady_number = get_unsteady_number(ras, plan_number)
geom_hdf = get_geometry_hdf(ras, plan_number)
selected = precipitation_boundary(ras, unsteady_number, case["precip_boundary"])
print("Verified selected native clock:", selected[["Use Fixed Start Time", "Fixed Start Date/Time"]].to_dict())
support = rainfall_support(geom_hdf, case["precip_boundary"])
terrain_inventory = get_terrain_inventory(ras)
dominant_terrain = terrain_inventory.iloc[0]
terrain_tif = Path(dominant_terrain["path"])
dominant_resolution_ft = float(dominant_terrain["x_resolution"])
consolidation_cell_size_ft = max(dominant_resolution_ft, MIN_CONSOLIDATION_CELL_SIZE_FT)
mesh_areas = HdfMesh.get_mesh_areas(geom_hdf)
pump_stations = get_pump_stations(geom_hdf, mesh_areas)
bounds = case["mrms_bounds"] or tuple(HdfProject.get_project_bounds_latlon(geom_hdf, buffer_percent=25.0))
precip_root = project_path / "P"
original_source_audit = pd.read_csv(precip_root / "source_grib_audit.csv")
mrms_files = [Path(value) for value in original_source_audit["file"]]
for path, expected_hash in zip(mrms_files, original_source_audit["sha256"]):
assert sha256_file(path) == expected_hash
mrms_stack = PrecipMrms.load_grib2_stack(mrms_files, bounds=bounds)
source_audit = audit_mrms_source(mrms_files, mrms_stack)
source_hyetograph = PrecipMrms.to_hyetograph(mrms_stack)
precision_audit = validate_source_hyetograph(source_hyetograph, mrms_stack)
precision_audit.to_csv(precip_root / "source_precision_final.csv", index=False)
hyetograph = align_mrms_hyetograph(source_hyetograph, case["event_start"], case["event_last_qpe"], case["sim_end"])
retained = pd.read_csv(precip_root / "aligned_hyetograph.csv", parse_dates=["time"])
pd.testing.assert_frame_equal(hyetograph, retained, check_dtype=False, atol=1e-12, rtol=0)
print(f"Source support: {source_audit.finite_spatial_cells.min()} of {source_audit.total_spatial_cells.iloc[0]} cells, unchanged over {len(source_audit)} frames; bounds={bounds}")
print(f"Raw aligned storm total: {hyetograph.cumulative_depth.iloc[-1]:.9f} in; rectangular arithmetic mean, not polygon weighting.")
display_preview(hyetograph)
aligned_stack = mrms_stack.sel(time=slice(case["event_start"] + timedelta(hours=1), case["event_last_qpe"]))
mrms_dss = precip_root / f"{case['case_id']}_mrms_qpe.dss"
dss_grid_audit = validate_dss_precipitation_grids(mrms_dss, catalog_dss(mrms_dss), case, len(mrms_files))
rows, runtime_rows, qualifications = [], [], []
hdfs = {}
for condition, label in (("baseline", "No-rain baseline"), ("event", "MRMS boundary event")):
receipt_path = precip_root / f"{condition}_compute_receipt.json"
receipt = json.loads(receipt_path.read_text())
assert receipt["success"] and receipt["completion_verified"]
assert receipt["process_inventory"]["complete"] and not receipt["process_inventory"]["matched"]
assert receipt["execution_details"]["solver_quiescence_confirmed"]
assert receipt["execution_details"]["actual_engine_provenance_confirmed"]
for relative, expected_hash in receipt["source_sha256"].items():
assert SOURCE_SHA256[relative] == expected_hash, f"Native producer source changed: {relative}"
for original, expected_hash in receipt["input_sha256_before"].items():
assert sha256_file(precip_root / f"{condition}_inputs" / Path(original).name) == expected_hash
preserve_after_input_snapshots(receipt, precip_root, condition)
if condition == "baseline":
hdf = precip_root / f"{case['case_id']}_no_rain_baseline.p{plan_number}.hdf"
else:
current = ras.project_folder / f"{ras.project_name}.p{plan_number}.hdf"
assert sha256_file(current) == receipt["hdf_sha256"]
hdf = precip_root / f"event_retained.p{plan_number}.hdf"
if not hdf.exists(): shutil.copy2(current, hdf)
assert sha256_file(hdf) == receipt["hdf_sha256"]
hdfs[condition] = hdf
expected = pd.read_csv(precip_root / f"{condition}_serialized.csv", parse_dates=["time"])
intended = zero_hyetograph_like(hyetograph) if condition == "baseline" else hyetograph
assert pd.DatetimeIndex(expected.time.iloc[1:]).equals(pd.DatetimeIndex(intended.time)), "Intended and serialized absolute times differ"
assert np.max(np.abs(expected.cumulative_depth.to_numpy()[1:] - intended.cumulative_depth.to_numpy())) <= .005 + 1e-12
# Independently re-read the immutable native input snapshot via the public table API.
u_path = precip_root / f"{condition}_inputs_after" / f"{ras.project_name}.u{unsteady_number}"
tables = RasUnsteady.extract_boundary_and_tables(u_path, ras_object=ras)
selected_tables = tables.loc[tables["Storage Area Name"].eq(case["precip_boundary"]) & tables["Tables"].map(lambda item: "Precipitation Hydrograph" in item)]
assert len(selected_tables) == 1
values = selected_tables.iloc[0]["Tables"]["Precipitation Hydrograph"]["Value"].to_numpy(dtype=float)
assert np.array_equal(values, expected.incremental_depth.to_numpy())
audit = validate_hdf_precipitation(hdf, case["precip_boundary"], expected, condition, case, support, precip_root)
rows.append(audit)
saved_messages = json.loads((precip_root / f"{label.replace(' ', '_')}_runtime_messages.json").read_text())
messages = HdfResultsPlan.get_compute_messages_hdf_only(hdf)
assert messages == saved_messages["hdf_compute_messages"]
runtime_rows.append(summarize_runtime_messages(messages, saved_messages["bco_messages"], condition))
qualifications.append({"project": case["project"],
"hydraulic_acceptance": "not_established",
"iteration_diagnostic_lines": runtime_rows[-1]["iteration_diagnostic_lines"],
"volume_accounting_lines": runtime_rows[-1]["volume_accounting_lines"],
"native_started_epoch": receipt["started_epoch"],
"native_source_module_sha256": dict(receipt["source_sha256"]),
"native_runtime_seconds": receipt["finished_epoch"]-receipt["started_epoch"],
"native_hdf_sha256": receipt["hdf_sha256"], "native_executable_sha256": receipt["execution_details"]["selected_executable_sha256"],
**audit.iloc[0].to_dict()})
plan_hdf, baseline_plan_hdf = hdfs["event"], hdfs["baseline"]
precipitation_audit = pd.concat(rows, ignore_index=True)
runtime_audit = pd.DataFrame([{k:v for k,v in row.items() if k != "diagnostic_lines"} for row in runtime_rows])
inventory = RasCmdr.inspect_plan_processes(plan_number, ras_object=ras)
assert inventory.complete and not inventory.matched
display(precipitation_audit)
display(runtime_audit)
# Missing variables are an evidence limit, never an assertion of zero rain.
for mesh in mesh_areas.mesh_name:
if mesh == case["precip_boundary"]: continue
for condition, hdf in hdfs.items():
try:
rate = HdfResultsMesh.get_mesh_timeseries(hdf, mesh, "Cell Precipitation Rate", truncate=False)
cumulative = HdfResultsMesh.get_mesh_timeseries(hdf, mesh, "Cell Cumulative Precipitation Depth", truncate=False)
except (KeyError, ValueError) as error:
print(f"Non-target mesh {mesh}, {condition}: precipitation variables unavailable; not qualified as zero ({error}).")
else:
assert np.isfinite(rate.values).all() and np.isfinite(cumulative.values).all()
assert np.all(rate.values == 0) and np.all(cumulative.values == 0)
print(f"Non-target mesh {mesh}, {condition}: precipitation arrays are zero.")
return locals()
def physical_support_from_topology(topology, polygon_ids):
"""Validate native endpoint cycles independently of precipitation values."""
from collections import Counter
from shapely.geometry import LineString, Point
from shapely.ops import polygonize
count = int(topology["n_cells"])
all_ids = np.arange(count, dtype=int)
face_counts = np.asarray(topology["cell_face_info"])[:, 1]
assert len(face_counts) == count and np.all(face_counts >= 1), "Native topology contains a cell with no faces"
physical = all_ids[face_counts > 1]
boundary = all_ids[face_counts == 1]
polygon_ids = set(int(i) for i in polygon_ids)
assert polygon_ids.issubset(set(physical)), "Polygon API cell is not a native multi-face cell"
rows = []
for cell in physical:
offset, n_faces = topology["cell_face_info"][cell, :2]
assert n_faces >= 3, "Physical endpoint cycle must have at least three faces"
faces = np.asarray(topology["cell_face_values"][offset:offset+n_faces, 0], dtype=int)
assert len(np.unique(faces)) == n_faces, "Physical endpoint cycle contains duplicate faces"
endpoints = np.asarray(topology["face_facepoints"])[faces]
degree = Counter(int(vertex) for pair in endpoints for vertex in pair)
assert degree and all(value == 2 for value in degree.values()), "Physical endpoint cycle is open or branched"
vertices = np.asarray(topology["facepoint_coords"])[endpoints]
assert np.isfinite(vertices).all(), "Physical endpoint cycle contains nonfinite coordinates"
rings = list(polygonize([LineString(pair) for pair in vertices]))
assert len(rings) == 1 and rings[0].is_valid and rings[0].area > 0, "Physical endpoint cycle must form one valid positive-area polygon"
assert np.isfinite(topology["cell_min_elev"][cell]), "Physical cell has nonfinite minimum elevation"
center_covered = bool(rings[0].covers(Point(topology["cell_centers"][cell])))
omitted = int(cell) not in polygon_ids
if omitted:
assert center_covered, "Reconstructed omitted cell does not contain its native center"
rows.append({"cell_id": int(cell), "faces": int(n_faces), "endpoint_area": float(rings[0].area),
"center_covered": center_covered, "polygon_api_omission": omitted})
physical_set = set(physical)
for cell in boundary:
offset, n_faces = topology["cell_face_info"][cell, :2]
face = int(topology["cell_face_values"][offset, 0])
neighbors = topology["face_cells"][face]
assert cell in neighbors and len(set(neighbors).intersection(physical_set)) == 1, "Complement face must join one physical cell"
return physical, boundary, all_ids, pd.DataFrame(rows)
def reconstruct_physical_support(plan_hdf, mesh_name, polygon_support, output):
topology = HdfMesh.get_mesh_sloped_topology(plan_hdf, mesh_name)
polygons = HdfMesh.get_mesh_cell_polygons(plan_hdf)
same_hdf_ids = np.sort(polygons.loc[polygons.mesh_name.eq(mesh_name), "cell_id"].to_numpy(dtype=int))
assert np.array_equal(same_hdf_ids, polygon_support[0]), "Geometry polygon IDs changed between source and preserved plan HDF"
physical, boundary, all_ids, audit = physical_support_from_topology(topology, same_hdf_ids)
omissions = audit.loc[audit.polygon_api_omission]
if not omissions.empty:
print("Detailed polygon API omissions reconstructed only for rainfall support from native endpoint cycles:")
display(omissions)
print("The detailed face perimeter remains invalid; model geometry is unchanged. This is an explicit geometry-reader limitation.")
audit.to_csv(output / f"{Path(plan_hdf).stem}_endpoint_geometry_audit.csv", index=False)
print(f"Native endpoint-cycle support: {len(physical)} physical cells; {len(boundary)} one-face boundary cells.")
return physical, boundary, all_ids
def align_animation_precipitation(aligned_stack, flood_times, origin, last_qpe):
"""Align display panels only; never use repeated display frames as new rainfall."""
import xarray as xr
times = pd.DatetimeIndex(flood_times).as_unit("ns")
source_times = pd.DatetimeIndex(aligned_stack.time.values).as_unit("ns")
origin, last_qpe = pd.Timestamp(origin), pd.Timestamp(last_qpe)
assert times.is_monotonic_increasing and not times.has_duplicates, "Display times must be unique and increasing"
assert len(times) and times[0] >= origin, "Display times precede boundary origin"
assert source_times.equals(pd.date_range(origin + pd.Timedelta(hours=1), last_qpe, freq="h")), "Observed source intervals do not cover the display window"
frames, rows = [], []
for timestamp in times:
if timestamp == origin or timestamp > last_qpe:
template = aligned_stack.isel(time=0, drop=True)
frame = xr.where(np.isfinite(template), 0.0, np.nan)
source_index, interval_start, interval_end = None, pd.NaT, pd.NaT
mode = "initial_zero" if timestamp == origin else "dry_tail"
else:
source_index = int(np.searchsorted(source_times.asi8, timestamp.value, side="left"))
frame = aligned_stack.isel(time=source_index, drop=True).copy(deep=True)
interval_end = source_times[source_index]
interval_start = interval_end - pd.Timedelta(hours=1)
mode = "observed_interval"
frames.append(frame)
rows.append({"flood_time": timestamp, "source_interval_start": interval_start,
"source_interval_end": interval_end, "display_mode": mode, "source_index": source_index})
display_stack = xr.concat(frames, dim=pd.DatetimeIndex(times, name="time"))
display_stack.attrs.update(aligned_stack.attrs)
display_stack.attrs["display_role"] = "spatial source illustration; uniform model boundary"
display_stack.attrs["display_alignment"] = "hourly interval covering flood time; exact interval end uses preceding interval"
return display_stack, pd.DataFrame(rows)
Case Runner¶
Python
def run_native_case(case: dict) -> dict:
print("=" * 88)
print(case["display_name"])
print(case["event_note"])
project_path = RUN_ROOT / "m" / case["project"]
if not project_path.exists():
project_path = RasExamples.extract_project(case["project"], output_path=RUN_ROOT / "m")
# A new run root is required once either baseline/event compute has been attempted.
if list((project_path / "P").glob("*_compute_receipt.json")):
raise RuntimeError("Retain prior evidence and select a new RAS_COMMANDER_EXAMPLE_RUN_ROOT")
assert max(len(str(p.resolve())) for p in project_path.rglob("*")) <= 240
ras = init_ras_project(
project_path,
ras_version=case["ras_version"],
ras_object=RasPrj(),
load_results_summary=False,
hide_intro=True,
)
assert RasPrj.get_project_units(RasPrj.find_ras_prj(ras.project_folder)) == "ft", "This example requires a US-customary project"
plan_number = case["plan_number"]
unsteady_number = get_unsteady_number(ras, plan_number)
geom_hdf = get_geometry_hdf(ras, plan_number)
selected_boundary = precipitation_boundary(ras, unsteady_number, case["precip_boundary"])
print("Selected boundary clock:", selected_boundary[["Use Fixed Start Time", "Fixed Start Date/Time"]].to_dict())
support = rainfall_support(geom_hdf, case["precip_boundary"])
terrain_inventory = get_terrain_inventory(ras)
if terrain_inventory.empty:
raise FileNotFoundError(f"No active RasMapper terrain TIFFs found for {case['display_name']}")
dominant_terrain = terrain_inventory.iloc[0]
terrain_tif = Path(dominant_terrain["path"])
terrain_units = str(dominant_terrain["linear_units"] or "").lower()
if "foot" not in terrain_units:
raise ValueError(f"Expected a foot-based terrain CRS, found {terrain_units!r}")
x_resolution_ft = float(dominant_terrain["x_resolution"])
y_resolution_ft = float(dominant_terrain["y_resolution"])
if not np.isclose(x_resolution_ft, y_resolution_ft):
raise ValueError(
f"Expected square dominant terrain cells, found "
f"{x_resolution_ft:.6f} x {y_resolution_ft:.6f} ft"
)
dominant_resolution_ft = x_resolution_ft
consolidation_cell_size_ft = max(
dominant_resolution_ft,
MIN_CONSOLIDATION_CELL_SIZE_FT,
)
mesh_areas = HdfMesh.get_mesh_areas(geom_hdf)
mesh_cell_counts = (
HdfMesh.get_mesh_cell_points(geom_hdf)
.groupby("mesh_name")
.size()
.rename("mesh_cells")
.reset_index()
)
pump_stations = get_pump_stations(geom_hdf, mesh_areas)
project_bounds = tuple(
float(value)
for value in HdfProject.get_project_bounds_latlon(
geom_hdf,
buffer_percent=25.0,
)
)
bounds = case["mrms_bounds"] or project_bounds
print(f"Project folder: {display_path(ras.project_folder)}")
print(f"Plan: {plan_number}; unsteady file: {unsteady_number}; geometry HDF: {geom_hdf.name}")
print(f"Mesh areas: {mesh_areas['mesh_name'].tolist()}")
print(f"HEC-RAS execution version: {case['ras_version']}")
display(mesh_cell_counts)
print(f"Pump stations: {pump_stations['Name'].tolist() if 'Name' in pump_stations else len(pump_stations)}")
boundary_columns = [
column
for column in [
"boundary_condition_number", "area_2d", "bc_line_name",
"bc_type", "hydrograph_type", "Interval",
"hydrograph_num_values", "Use DSS",
]
if column in ras.boundaries_df.columns
]
boundary_context = ras.boundaries_df.loc[
ras.boundaries_df["unsteady_number"].astype(str).str.zfill(2) == unsteady_number,
boundary_columns,
].copy()
print("Boundary conditions retained in both the no-rain baseline and MRMS event:")
display(boundary_context)
terrain_display = terrain_inventory.drop(columns="path").copy()
terrain_display["valid_coverage_sq_mi"] = (
terrain_display["valid_coverage_area"] / 5280.0**2
)
print("Terrain mosaic order (1 is written first; later local detail takes precedence):")
display(terrain_display.drop(columns="valid_coverage_area").round(3))
print(
f"Dominant terrain: {terrain_tif.name}; native cell: {dominant_resolution_ft:.3f} ft; "
f"consolidated cell: {consolidation_cell_size_ft:.3f} ft"
)
print(f"Project WGS84 bounds: {project_bounds}")
print(f"MRMS download bounds: {bounds}")
precip_root = ras.project_folder / "P"
precip_root.mkdir(parents=True, exist_ok=True)
clip_shp = make_clip_shp(bounds, precip_root / f"{case['case_id']}_bounds.shp", case["display_name"])
catalog = PrecipMrms.catalog(
bounds=bounds,
start_date=case["event_start"],
end_date=case["event_last_qpe"],
source="noaa_s3",
)
expected_files = int((case["event_last_qpe"] - case["event_start"]).total_seconds() // 3600) + 1
assert len(catalog) >= expected_files, f"Expected at least {expected_files} MRMS files, found {len(catalog)}"
display_preview(catalog[["valid_time", "product", "archive_product", "filename", "size_bytes", "compressed"]])
mrms_files = PrecipMrms.download(
bounds=bounds,
start_time=case["event_start"],
end_time=case["event_last_qpe"],
output_dir=precip_root / "grib2",
source="noaa_s3",
)
assert len(mrms_files) >= expected_files
print(f"Downloaded {len(mrms_files)} MRMS GRIB2 files")
mrms_dss_path = precip_root / f"{case['case_id']}_mrms_qpe.dss"
mrms_dss_path.unlink(missing_ok=True)
mrms_dss = PrecipMrms.to_dss(
mrms_files,
mrms_dss_path,
clip_shp=clip_shp,
dss_parts={"B": case["dss_b"], "F": case["dss_f"]},
timeout=1800,
)
assert mrms_dss.exists() and mrms_dss.stat().st_size > 0
dss_catalog = catalog_dss(mrms_dss)
display_preview(dss_catalog)
dss_grid_audit = validate_dss_precipitation_grids(
mrms_dss, dss_catalog, case, expected_files
)
display_preview(dss_grid_audit)
mrms_stack = PrecipMrms.load_grib2_stack(mrms_files, bounds=bounds)
source_audit = audit_mrms_source(mrms_files, mrms_stack)
source_audit.to_csv(precip_root / "source_grib_audit.csv", index=False)
dss_grid_audit.to_csv(precip_root / "source_dss_audit.csv", index=False)
source_hyetograph = PrecipMrms.to_hyetograph(mrms_stack)
validate_source_hyetograph(source_hyetograph, mrms_stack)
source_hyetograph.to_csv(precip_root / "source_hyetograph.csv", index=False)
hyetograph = align_mrms_hyetograph(source_hyetograph, case["event_start"], case["event_last_qpe"], case["sim_end"])
hyetograph.to_csv(precip_root / "aligned_hyetograph.csv", index=False)
baseline_hyetograph = zero_hyetograph_like(hyetograph)
print("Pre-start source accumulations excluded from forcing:")
display(source_hyetograph.loc[source_hyetograph["time"] <= case["event_start"]])
print(f"Explicit zero-rain tail: after {case['event_last_qpe']} through {case['sim_end']} UTC-equivalent model time")
display_preview(hyetograph)
print(f"MRMS stack shape: {mrms_stack.shape}")
print(f"Aligned spatial-mean event depth before serialization: {hyetograph['cumulative_depth'].iloc[-1]:.9f} inches")
# The map shows the source accumulation over the retained rainfall window, excluding pre-start data.
aligned_stack = mrms_stack.sel(time=slice(case["event_start"] + timedelta(hours=1), case["event_last_qpe"]))
# Visualization is a separate, rerunnable phase.
RasPlan.update_simulation_date(
plan_number,
case["event_start"],
case["sim_end"],
ras_object=ras,
)
RasPlan.update_plan_intervals(
plan_number,
output_interval=OUTPUT_INTERVAL,
instantaneous_interval=OUTPUT_INTERVAL,
mapping_interval=OUTPUT_INTERVAL,
ras_object=ras,
)
print(f"Updated precipitation boundary: {case['precip_boundary']}")
print(f"Simulation window: {case['event_start']} to {case['sim_end']}")
for key in ["Output Interval", "Instantaneous Interval", "Mapping Interval"]:
print(f"{key}: {RasPlan.get_plan_value(plan_number, key, ras_object=ras)}")
hdf_output_variables = [
"Cell Hydraulic Depth",
"Cell Invert Depth (WSE - Cell Min Elev)",
"Water Surface",
"Cell Precipitation Rate",
"Cell Cumulative Precipitation Depth",
]
plan_hdf = ras.project_folder / f"{ras.project_name}.p{plan_number}.hdf"
RasUnsteady.set_precipitation_hyetograph(
unsteady_number,
baseline_hyetograph,
boundary_name=case["precip_boundary"],
ras_object=ras,
)
baseline_serialized = readback_hyetograph(ras, unsteady_number, case, baseline_hyetograph, "baseline", precip_root)
baseline_compute = compute_with_receipt(ras, case, hdf_output_variables, "baseline", precip_root)
assert bool(baseline_compute), f"No-rain baseline failed for {case['display_name']}"
baseline_runtime_audit = audit_runtime_messages(
ras, plan_number, plan_hdf, "No-rain baseline"
)
assert baseline_runtime_audit["complete_process"]
baseline_plan_hdf = precip_root / f"{case['case_id']}_no_rain_baseline.p{plan_number}.hdf"
baseline_plan_hdf.unlink(missing_ok=True)
shutil.copy2(plan_hdf, baseline_plan_hdf)
assert sha256_file(baseline_plan_hdf) == json.loads((precip_root / "baseline_compute_receipt.json").read_text())["hdf_sha256"]
baseline_precip_audit = validate_hdf_precipitation(
baseline_plan_hdf, case["precip_boundary"], baseline_serialized, "baseline", case, support, precip_root
)
RasUnsteady.set_precipitation_hyetograph(
unsteady_number,
hyetograph,
boundary_name=case["precip_boundary"],
ras_object=ras,
)
serialized_hyetograph = readback_hyetograph(ras, unsteady_number, case, hyetograph, "event", precip_root)
compute_result = compute_with_receipt(ras, case, hdf_output_variables, "event", precip_root)
assert bool(compute_result), f"MRMS-event compute failed for {case['display_name']}"
event_runtime_audit = audit_runtime_messages(
ras, plan_number, plan_hdf, "MRMS boundary event"
)
assert event_runtime_audit["complete_process"]
runtime_audit = pd.DataFrame([
{key: value for key, value in row.items() if key != "diagnostic_lines"}
for row in (baseline_runtime_audit, event_runtime_audit)
])
display(runtime_audit)
for row in (baseline_runtime_audit, event_runtime_audit):
if row["diagnostic_lines"]:
print(f"{row['condition']} HDF/BCO diagnostic lines requiring manual review:")
for line in row["diagnostic_lines"]:
print(f"- {line}")
event_precip_audit = validate_hdf_precipitation(
plan_hdf, case["precip_boundary"], serialized_hyetograph, "event", case, support, precip_root
)
precipitation_audit = pd.concat(
[baseline_precip_audit, event_precip_audit], ignore_index=True
)
precipitation_audit.to_csv(precip_root / "precipitation_qualification.csv", index=False)
display(precipitation_audit)
return locals()
def render_case(native):
case = native["case"]
ras = native["ras"]
plan_number = native["plan_number"]
geom_hdf = native["geom_hdf"]
terrain_inventory = native["terrain_inventory"]
terrain_tif = native["terrain_tif"]
dominant_resolution_ft = native["dominant_resolution_ft"]
consolidation_cell_size_ft = native["consolidation_cell_size_ft"]
mesh_areas = native["mesh_areas"]
pump_stations = native["pump_stations"]
bounds = native["bounds"]
precip_root = native["precip_root"]
mrms_files = native["mrms_files"]
mrms_dss = native["mrms_dss"]
dss_grid_audit = native["dss_grid_audit"]
hyetograph = native["hyetograph"]
aligned_stack = native["aligned_stack"]
plan_hdf = native["plan_hdf"]
baseline_plan_hdf = native["baseline_plan_hdf"]
precipitation_audit = native["precipitation_audit"]
runtime_audit = native["runtime_audit"]
plot_accumulation_map(case, aligned_stack, hyetograph, mesh_areas, pump_stations)
timestamps = RasProcess.get_plan_timestamps(plan_number, ras_object=ras)
assert timestamps, "No HEC-RAS output timesteps found"
available_animation_timestamps = select_animation_timestamps(
timestamps,
case["animation_start"],
case["animation_end"],
)
if len(available_animation_timestamps) > ANIMATION_MAX_FRAMES:
frame_indices = np.linspace(
0,
len(available_animation_timestamps) - 1,
ANIMATION_MAX_FRAMES,
dtype=int,
)
animation_timestamps = [
available_animation_timestamps[index] for index in frame_indices
]
else:
animation_timestamps = available_animation_timestamps
print("Rendering previously qualified native result; no solver execution in this cell.")
print(f"Plan HDF: {display_path(plan_hdf)}")
print(f"Total output timesteps: {len(timestamps)}")
print(f"Available consecutive 5-minute timesteps: {len(available_animation_timestamps)}")
print(f"Selected map and animation frames: {len(animation_timestamps)}")
print(f"Selected frame coverage: {animation_timestamps[0]} to {animation_timestamps[-1]}")
event_hdf_depth_max = hdf_hydraulic_depth_maxima(plan_hdf, case["precip_boundary"])
event_hdf_invert_max = hdf_invert_depth_maxima(plan_hdf, geom_hdf, case["precip_boundary"])
depth_frame_dir = (ras.project_folder / "HdfFrames").resolve()
depth_results = HdfResultsMesh.export_depth_rasters_at_times(
plan_hdf, animation_timestamps, depth_frame_dir,
mesh_name=case["precip_boundary"], geom_hdf_path=geom_hdf,
resolution=consolidation_cell_size_ft, method="nearest",
)
depth_tile_frames = [[Path(depth_results[timestamp])] for timestamp in animation_timestamps]
depth_tile_audit = pd.DataFrame({"timestamp": [parse_ras_timestamp(t) for t in animation_timestamps],
"source_rasters": 1, "source": "HDF nearest-cell invert-depth raster"})
mapper_frame_audit = pd.DataFrame()
hdf_gap_fill_count = 0
raster_source = "Consistent HDF nearest-cell invert-depth rasters (not native Mapper equivalence)"
print(raster_source)
depth_times = [parse_ras_timestamp(timestamp) for timestamp in animation_timestamps]
depth_stack = PrecipMrms.load_stored_map_stack(
depth_tile_frames,
times=depth_times,
max_frames=ANIMATION_MAX_FRAMES,
name="flood_inundation",
units="ft",
cell_size=consolidation_cell_size_ft,
resampling="nearest",
)
consolidated_frame_dir = (
ras.project_folder
/ "HdfConsolidated"
).resolve()
consolidated_depth_rasters, consolidated_raster_audit = write_consolidated_depth_rasters(
depth_stack,
consolidated_frame_dir,
)
assert len(consolidated_depth_rasters) == depth_stack.sizes["time"]
depth_summary = summarize_depth_stack(depth_stack)
depth_summary["native_hydraulic_max_ft"] = depth_summary["timestamp"].map(
lambda value: float(event_hdf_depth_max.get(pd.Timestamp(value), np.nan))
)
depth_summary["native_invert_max_ft"] = depth_summary["timestamp"].map(event_hdf_invert_max)
if depth_summary["native_hydraulic_max_ft"].isna().any():
raise RuntimeError("Selected raster timestamps do not align with plan-HDF depth output")
depth_summary["raster_wetness_matches_hydraulic"] = (
(depth_summary["raster_invert_max_ft"] > 0.01)
== (depth_summary["native_hydraulic_max_ft"] > 0.01)
)
if not depth_summary["raster_wetness_matches_hydraulic"].all():
failed = depth_summary.loc[
~depth_summary["raster_wetness_matches_hydraulic"],
["timestamp", "raster_invert_max_ft", "native_hydraulic_max_ft"],
]
raise RuntimeError(f"Raster invert-depth/native hydraulic-depth wetness mismatch:\n{failed}")
total_depth_tiles = sum(len(frame) for frame in depth_tile_frames)
# Preserve full frame-level diagnostics alongside the original Mapper attempt.
depth_tile_audit.to_csv(precip_root / "hdf_frame_source_audit.csv", index=False)
consolidated_raster_audit.to_csv(precip_root / "hdf_consolidated_raster_audit.csv", index=False)
depth_summary.to_csv(precip_root / "hdf_depth_summary.csv", index=False)
display_preview(depth_tile_audit)
display_preview(consolidated_raster_audit)
display_preview(depth_summary)
print(
f"{raster_source}: {len(depth_tile_frames)} source timesteps, "
f"{total_depth_tiles} HDF source frames written to "
f"{len(consolidated_depth_rasters)} single {consolidation_cell_size_ft:.3f}-ft rasters"
)
print(f"Diagnostic raster folder: {display_path(consolidated_frame_dir)}")
print(
"HDF rasters per selected timestep: "
f"min={depth_tile_audit['source_rasters'].min()}, "
f"median={depth_tile_audit['source_rasters'].median():.0f}, "
f"max={depth_tile_audit['source_rasters'].max()}"
)
print("Static invert-depth maps use the actual raster maximum. Animation APIs cap colors at the 98th percentile of positive values (minimum 0.1); this is a display cap, not maximum invert depth.")
precip_video = PrecipMrms.animate_precipitation(
aligned_stack,
fresh_artifact_path(ARTIFACT_ROOT / f"mrms_precipitation_{case['case_id']}.mp4"),
bounds=bounds,
mesh_boundary=mesh_areas,
pump_stations=pump_stations,
add_basemap=False,
crs="EPSG:4326",
units="in/hr",
fps=FPS,
dpi=VIDEO_DPI,
title=f"{case['project']}: hourly MRMS QPE\nColor capped at positive P98 (minimum 0.1 in/hr)",
)
flood_video = PrecipMrms.animate_flood_inundation(
depth_stack,
fresh_artifact_path(ARTIFACT_ROOT / f"flood_inundation_depth_{case['case_id']}_5min.mp4"),
mesh_boundary=mesh_areas,
pump_stations=pump_stations,
add_basemap=False,
units="ft",
fps=FPS,
dpi=VIDEO_DPI,
title=f"{case['project']}: HDF nearest-cell invert depth\nColor capped at positive P98 (minimum 0.1 ft)",
)
combined_display_stack, combined_alignment = align_animation_precipitation(
aligned_stack, depth_stack.time.values, case["event_start"], case["event_last_qpe"])
combined_alignment.to_csv(precip_root / "combined_display_interval_audit.csv", index=False)
display_preview(combined_alignment)
print("Combined panel: hourly interval covering flood time; exact endpoint uses the preceding interval; initial and dry-tail frames zero. Spatial source illustration, uniform model boundary.")
# The animation API creates its axes internally; use scoped automatic offsets
# to keep geographic tick labels readable without mutating pyplot or the API.
with plt.rc_context({"axes.formatter.useoffset": True, "axes.formatter.offset_threshold": 2}):
combined_video = PrecipMrms.animate_combined(
combined_display_stack,
depth_stack,
fresh_artifact_path(ARTIFACT_ROOT / f"combined_precip_flood_{case['case_id']}_5min.mp4"),
precip_bounds=bounds,
mesh_boundary=mesh_areas,
pump_stations=pump_stations,
add_basemap=False,
precip_crs="EPSG:4326",
fps=FPS,
dpi=VIDEO_DPI,
title=f"{case['project']}: hourly interval covering flood time\nSpatial source illustration; uniform model boundary | HDF invert depth\nColors capped at positive P98 (minimum 0.1)",
)
videos = [
("Precipitation only", precip_video),
("Flood inundation only", flood_video),
("Combined precipitation and flood", combined_video),
]
for _, path in videos:
assert path.exists() and path.stat().st_size > 0
print(f"Video: {display_path(path)} ({path.stat().st_size:,} bytes)")
plot_depth_response(case, depth_summary, depth_stack, mesh_areas, pump_stations)
reference_cell_id, reference_selection = select_reference_cell_by_depth_delta(
plan_hdf,
baseline_plan_hdf,
case["precip_boundary"],
)
display(reference_selection)
event_reference_response, reference_point, _ = build_reference_response(
plan_hdf,
geom_hdf,
case["precip_boundary"],
reference_cell_id=reference_cell_id,
)
baseline_reference_response, _, _ = build_reference_response(
baseline_plan_hdf,
geom_hdf,
case["precip_boundary"],
reference_cell_id=reference_cell_id,
)
reference_summary = pd.DataFrame(
[
{
"case": case["case_id"],
"mesh_name": case["precip_boundary"],
"reference_cell_id": reference_cell_id,
"baseline_peak_hydraulic_depth_ft": float(baseline_reference_response["hydraulic_depth_ft"].max()),
"event_peak_hydraulic_depth_ft": float(event_reference_response["hydraulic_depth_ft"].max()),
"peak_hydraulic_depth_delta_ft": float(
np.nanmax(
event_reference_response["hydraulic_depth_ft"].to_numpy()
- baseline_reference_response["hydraulic_depth_ft"].to_numpy()
)
),
"baseline_peak_wse_ft": float(baseline_reference_response["wse_ft"].max()),
"event_peak_wse_ft": float(event_reference_response["wse_ft"].max()),
"peak_wse_delta_ft": float(
np.nanmax(
event_reference_response["wse_ft"].to_numpy()
- baseline_reference_response["wse_ft"].to_numpy()
)
),
"event_peak_precip_in_hr": float(event_reference_response["precip_in_hr"].max()),
}
]
)
display(reference_summary)
plot_reference_response(
case,
hyetograph,
event_reference_response,
baseline_reference_response,
reference_point,
depth_stack,
mesh_areas,
pump_stations,
)
event_pump_operation = get_pump_operation(plan_hdf, pump_stations)
baseline_pump_operation = get_pump_operation(baseline_plan_hdf, pump_stations)
pump_summary = plot_pump_operation(
case, event_pump_operation, baseline_pump_operation
)
if not pump_summary.empty:
display(pump_summary)
display_case_video_artifacts(case, videos)
return {
"case": case["case_id"],
"project_folder": ras.project_folder,
"plan_hdf": plan_hdf,
"baseline_plan_hdf": baseline_plan_hdf,
"mrms_dss": mrms_dss,
"dss_grid_audit": dss_grid_audit,
"mrms_files": mrms_files,
"depth_tile_frames": depth_tile_frames,
"consolidated_depth_rasters": consolidated_depth_rasters,
"depth_tile_audit": depth_tile_audit,
"consolidated_raster_audit": consolidated_raster_audit,
"depth_stack": depth_stack,
"depth_summary": depth_summary,
"reference_cell_id": reference_cell_id,
"reference_selection": reference_selection,
"reference_summary": reference_summary,
"event_reference_response": event_reference_response,
"baseline_reference_response": baseline_reference_response,
"event_pump_operation": event_pump_operation,
"baseline_pump_operation": baseline_pump_operation,
"precipitation_audit": precipitation_audit,
"runtime_audit": runtime_audit,
"mapper_frame_audit": mapper_frame_audit,
"pump_summary": pump_summary,
"videos": [path for _, path in videos],
"mesh_count": len(mesh_areas),
"pump_count": len(pump_stations),
"total_output_timesteps": len(timestamps),
"available_animation_timesteps": len(available_animation_timestamps),
"animation_timesteps": len(animation_timestamps),
"animation_frames": int(depth_stack.sizes["time"]),
"depth_tile_count": total_depth_tiles,
"consolidated_raster_count": len(consolidated_depth_rasters),
"dominant_terrain_resolution_ft": dominant_resolution_ft,
"consolidation_cell_size_ft": consolidation_cell_size_ft,
"hdf_gap_fill_count": hdf_gap_fill_count,
}
Native execution or explicit retained-result inspection¶
The default mode is compute: prepare fresh model copies, execute four sequential
native plans, then inspect the results. Set RAS_COMMANDER_917_MODE=inspect_retained
and the exact existing RAS_COMMANDER_EXAMPLE_RUN_ROOT only to inspect completed
runs. Missing receipts, changed sources/inputs/output hashes, failed native gates,
or active processes stop inspection. Existing files never silently skip compute.
The saved execution below uses inspect_retained for the four native runs completed earlier on 2026-09-27. The original producer source and exact output notebook are archived under that run root. This refactored default compute recipe has not been rerun through the solvers; final inspection and visualization were executed afresh. Original source/DSS preparation was pre-staged before those native runs. No solve was repeated solely to replace blocked basemap graphics.
The saved notebook runtime measures this inspection/render execution only. The four original native run durations and their total are separately emitted in the qualification summary. The refactored default compute recipe is executable but was not natively rerun after separating visualization; its original producer sources and continuation are retained.
The final retained inspection/render was executed again after visual corrections and explicit per-case project objects. Earlier complete inspection/render outputs and videos are archived; no native solve was repeated. The refactored compute recipe remains unexecuted by this saved inspect_retained run.
Python
native_cases = []
if EXECUTION_MODE == "compute":
run_native_case(CASES[0])
native_cases.append(inspect_retained_case(CASES[0]))
Text Only
INSPECT_RETAINED: Davis; no native solver launched by this cell.
Verified selected native clock: {'Use Fixed Start Time': 'False', 'Fixed Start Date/Time': ','}
Detailed polygon API support: 2705 cell polygons; 236 other cells require topology classification.
GRIB local product 209/6/37: NOAA MultiSensor_QPE_01H_Pass2, 1-hour mm.
ecCodes units may be unknown; units/period resolved from NOAA's product table, not inferred from cadence.
Source: https://www.nssl.noaa.gov/projects/mrms/operational/tables.php
Source precision: dtype=float32; max difference=4.30775973226e-09 in; max 4-epsilon bound=1.03587375098e-07 in; exact zero required for dry frames.
Source support: 56 of 56 cells, unchanged over 37 frames; bounds=(-121.78689244974507, 38.523870898941006, -121.70713010751233, 38.591248584093314)
Raw aligned storm total: 1.741071451 in; rectangular arithmetic mean, not polygon weighting.
Showing first and last 3 of 42 rows
| time | hour | incremental_depth | cumulative_depth | |
|---|---|---|---|---|
| 0 | 2022-12-31 01:00:00 | 1.0 | 0.022497 | 0.022497 |
| 1 | 2022-12-31 02:00:00 | 2.0 | 0.003375 | 0.025872 |
| 2 | 2022-12-31 03:00:00 | 3.0 | 0.039159 | 0.065031 |
| 39 | 2023-01-01 16:00:00 | 40.0 | 0.000000 | 1.741071 |
| 40 | 2023-01-01 17:00:00 | 41.0 | 0.000000 | 1.741071 |
| 41 | 2023-01-01 18:00:00 | 42.0 | 0.000000 | 1.741071 |
Text Only
2026-09-27 06:12:58 - ras_commander.RasUnsteady - INFO - Successfully extracted boundaries and tables from C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\m\Davis\P\baseline_inputs_after\DavisStormSystem.u01
2026-09-27 06:12:58 - ras_commander.RasUnsteady - INFO - Successfully extracted boundaries and tables from C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\m\Davis\P\event_inputs_after\DavisStormSystem.u01
Native endpoint-cycle support: 2705 physical cells; 236 one-face boundary cells.
Geometry-derived rainfall support: 2705 physical cells; 236 one-face boundary cells.
baseline: warning/error-pattern lines from HDF and BCO: 0; iteration-table lines: 0
baseline: complete retained HDF compute messages (numerical diagnostics require engineering review):
Plan: 'Full System ROM with Pump' (DavisStormSystem.p02)
Simulation started at: 27Sep2026 04:49:16 AM
Writing Plan GIS Data...
Completed Writing Plan GIS Data
Writing Geometry...
Terrain: Terrain (1)
Computing 2D Flow Area 'area2' tables: Property tables do not exist.
2D Flow Area 'area2' tables complete 1.14 sec
Computing 2D Flow Area 'DS Channel' tables: Property tables do not exist.
2D Flow Area 'DS Channel' tables complete 0.03 sec
Computing Pipe Network 'Davis' tables: Property tables do not exist.
Pipe Network 'Davis' tables complete 0.14 sec
Davis: Pipe Network property tables are current.
Completed Writing Geometry
Writing Event Conditions ...
Completed Writing Event Condition Data
Geometric Preprocessor HEC-RAS 7.0 April 2026
Finished Processing Geometry
Parallelization CPU Affinity ON (Performance Cores 0,1)
Performing Unsteady Flow Simulation HEC-RAS 7.0 April 2026
Unsteady Input Summary:
2D Area Cells Cores Solvers
area2 2705 2 Diffusion Wave Equation Set (fastest)
DS Channel 11 1 Diffusion Wave Equation Set (fastest)
Pipe Network Nodes Cores Courant Target Solvers
Davis 133 1 Diffusion Wave Equation Set (fastest)
Overall Volume Accounting Error in Acre Feet: 0.0000000000
Overall Volume Accounting Error as percentage: 0.0000000000
Please review "Computational Log File" output for volume accounting details
Writing Results to DSS
Finished Unsteady Flow Simulation
Computations Summary
Computation Task Time(hh:mm:ss)
Completing Geometry, Flow and Plan 3
Preprocessing Geometry <1
Unsteady Flow Computations 51
Complete Process 55
Computation Speed Simulation/Runtime
Unsteady Flow Computations 2940x
Complete Process 2734x
Native endpoint-cycle support: 2705 physical cells; 236 one-face boundary cells.
Geometry-derived rainfall support: 2705 physical cells; 236 one-face boundary cells.
event: warning/error-pattern lines from HDF and BCO: 0; iteration-table lines: 0
event: complete retained HDF compute messages (numerical diagnostics require engineering review):
Plan: 'Full System ROM with Pump' (DavisStormSystem.p02)
Simulation started at: 27Sep2026 04:50:14 AM
Writing Plan GIS Data...
Completed Writing Plan GIS Data
Writing Geometry...
Terrain: Terrain (1)
Computing 2D Flow Area 'area2' tables: Property tables do not exist.
2D Flow Area 'area2' tables complete 1.07 sec
Computing 2D Flow Area 'DS Channel' tables: Property tables do not exist.
2D Flow Area 'DS Channel' tables complete 0.03 sec
Davis: Pipe Network property tables are current.
Completed Writing Geometry
Writing Event Conditions ...
Completed Writing Event Condition Data
Geometric Preprocessor HEC-RAS 7.0 April 2026
Finished Processing Geometry
Parallelization CPU Affinity ON (Performance Cores 0,1)
Performing Unsteady Flow Simulation HEC-RAS 7.0 April 2026
Unsteady Input Summary:
2D Area Cells Cores Solvers
area2 2705 2 Diffusion Wave Equation Set (fastest)
DS Channel 11 1 Diffusion Wave Equation Set (fastest)
Pipe Network Nodes Cores Courant Target Solvers
Davis 133 1 Diffusion Wave Equation Set (fastest)
Overall Volume Accounting Error in Acre Feet: 0.07065
Overall Volume Accounting Error as percentage: 0.03423
Please review "Computational Log File" output for volume accounting details
Writing Results to DSS
Finished Unsteady Flow Simulation
Computations Summary
Computation Task Time(hh:mm:ss)
Completing Geometry, Flow and Plan 2
Preprocessing Geometry <1
Unsteady Flow Computations 3:48
Complete Process 3:51
Computation Speed Simulation/Runtime
Unsteady Flow Computations 662x
Complete Process 654x
| condition | timestamps | physical_cells_checked | boundary_cells_zero | absolute_start_utc | absolute_end_utc | expected_total_in | hdf_final_total_in | max_all_time_cell_error_in | native_tolerance_in | hdf_peak_rate_in_hr | absolute_time_qualified | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | baseline | 505 | 2705 | 236 | 2022-12-31 | 2023-01-01 18:00:00 | 0.00 | 0.00 | 0.000000e+00 | 9.536743e-07 | 0.00 | True |
| 1 | event | 505 | 2705 | 236 | 2022-12-31 | 2023-01-01 18:00:00 | 1.74 | 1.74 | 2.098084e-07 | 1.659393e-06 | 0.21 | True |
| condition | complete_process | bco_exists | diagnostic_line_count | iteration_diagnostic_lines | volume_accounting_lines | hydraulic_acceptance | |
|---|---|---|---|---|---|---|---|
| 0 | baseline | True | True | 0 | [] | [Overall Volume Accounting Error in Acre Feet:... | not_established |
| 1 | event | True | True | 0 | [] | [Overall Volume Accounting Error in Acre Feet:... | not_established |
Text Only
Non-target mesh DS Channel, baseline: precipitation variables unavailable; not qualified as zero (Path Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series/2D Flow Areas/DS Channel/Cell Precipitation Rate not found in HDF file).
Non-target mesh DS Channel, event: precipitation variables unavailable; not qualified as zero (Path Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series/2D Flow Areas/DS Channel/Cell Precipitation Rate not found in HDF file).
Python
if EXECUTION_MODE == "compute":
run_native_case(CASES[1])
native_cases.append(inspect_retained_case(CASES[1]))
Text Only
INSPECT_RETAINED: NewOrleansMetro; no native solver launched by this cell.
Verified selected native clock: {'Use Fixed Start Time': 'False', 'Fixed Start Date/Time': ','}
Detailed polygon API support: 19710 cell polygons; 776 other cells require topology classification.
GRIB local product 209/6/37: NOAA MultiSensor_QPE_01H_Pass2, 1-hour mm.
ecCodes units may be unknown; units/period resolved from NOAA's product table, not inferred from cadence.
Source: https://www.nssl.noaa.gov/projects/mrms/operational/tables.php
Source precision: dtype=float32; max difference=1.52041612633e-08 in; max 4-epsilon bound=4.53999972612e-07 in; exact zero required for dry frames.
Source support: 182 of 182 cells, unchanged over 24 frames; bounds=(np.float64(-90.18290001411526), np.float64(29.900389160096417), np.float64(-90.04639080882717), np.float64(30.04210617844562))
Raw aligned storm total: 4.095613113 in; rectangular arithmetic mean, not polygon weighting.
Showing first and last 3 of 30 rows
| time | hour | incremental_depth | cumulative_depth | |
|---|---|---|---|---|
| 0 | 2024-04-10 01:00:00 | 1.0 | 0.0 | 0.000000 |
| 1 | 2024-04-10 02:00:00 | 2.0 | 0.0 | 0.000000 |
| 2 | 2024-04-10 03:00:00 | 3.0 | 0.0 | 0.000000 |
| 27 | 2024-04-11 04:00:00 | 28.0 | 0.0 | 4.095613 |
| 28 | 2024-04-11 05:00:00 | 29.0 | 0.0 | 4.095613 |
| 29 | 2024-04-11 06:00:00 | 30.0 | 0.0 | 4.095613 |
Text Only
2026-09-27 06:13:12 - ras_commander.RasUnsteady - INFO - Successfully extracted boundaries and tables from C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\m\NewOrleansMetro\P\baseline_inputs_after\NewOrleansMetro.u01
Detailed polygon API omissions reconstructed only for rainfall support from native endpoint cycles:
| cell_id | faces | endpoint_area | center_covered | polygon_api_omission | |
|---|---|---|---|---|---|
| 16058 | 16058 | 5 | 10730.610567 | True | True |
Text Only
The detailed face perimeter remains invalid; model geometry is unchanged. This is an explicit geometry-reader limitation.
Native endpoint-cycle support: 19711 physical cells; 775 one-face boundary cells.
Geometry-derived rainfall support: 19711 physical cells; 775 one-face boundary cells.
baseline: warning/error-pattern lines from HDF and BCO: 0; iteration-table lines: 0
baseline: complete retained HDF compute messages (numerical diagnostics require engineering review):
Plan: '2D Surface Model w Pipes' (NewOrleansMetro.p01)
Simulation started at: 27Sep2026 04:56:59 AM
Writing Plan GIS Data...
Completed Writing Plan GIS Data
Writing Geometry...
Terrain: Terrain
Land Cover: LandCover
Computing 2D Flow Area 'NewOrleans Metro' tables: Property tables do not exist.
2D Flow Area 'NewOrleans Metro' tables complete 16.36 sec
Computing Pipe Network 'Base' tables: Property tables do not exist.
Pipe Network 'Base' tables complete 0.13 sec
Base: Pipe Network property tables are current.
Completed Writing Geometry
Writing Event Conditions ...
Computing IC for 2D Flow Area: NewOrleans Metro
Completed IC for 2D Flow Area: NewOrleans Metro
Completed Writing Event Condition Data
Geometric Preprocessor HEC-RAS 7.0 April 2026
Finished Processing Geometry
Parallelization CPU Affinity ON (Performance Cores 0,1)
Performing Unsteady Flow Simulation HEC-RAS 7.0 April 2026
Unsteady Input Summary:
2D Area Cells Cores Solvers
NewOrleans Metro 19711 2 SWE-ELM Equation Set (original)
Pipe Network Nodes Cores Courant Target Solvers
Base 64 2 0.90 SWE-ELM Equation Set (original)
Overall Volume Accounting Error in Acre Feet: 0.000103
Overall Volume Accounting Error as percentage: 0.0000095671
Please review "Computational Log File" output for volume accounting details
Writing Results to DSS
Finished Unsteady Flow Simulation
Computations Summary
Computation Task Time(hh:mm:ss)
Completing Geometry, Flow and Plan 19
Preprocessing Geometry <1
Unsteady Flow Computations 8:33
Complete Process 8:52
Computation Speed Simulation/Runtime
Unsteady Flow Computations 210x
Complete Process 203x
2026-09-27 06:13:17 - ras_commander.RasUnsteady - INFO - Successfully extracted boundaries and tables from C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\m\NewOrleansMetro\P\event_inputs_after\NewOrleansMetro.u01
Detailed polygon API omissions reconstructed only for rainfall support from native endpoint cycles:
| cell_id | faces | endpoint_area | center_covered | polygon_api_omission | |
|---|---|---|---|---|---|
| 16058 | 16058 | 5 | 10730.610567 | True | True |
Text Only
The detailed face perimeter remains invalid; model geometry is unchanged. This is an explicit geometry-reader limitation.
Native endpoint-cycle support: 19711 physical cells; 775 one-face boundary cells.
Geometry-derived rainfall support: 19711 physical cells; 775 one-face boundary cells.
event: warning/error-pattern lines from HDF and BCO: 0; iteration-table lines: 3
event: complete retained HDF compute messages (numerical diagnostics require engineering review):
Plan: '2D Surface Model w Pipes' (NewOrleansMetro.p01)
Simulation started at: 27Sep2026 05:16:51 AM
Writing Plan GIS Data...
Completed Writing Plan GIS Data
Writing Geometry...
Terrain: Terrain
Land Cover: LandCover
Computing 2D Flow Area 'NewOrleans Metro' tables: Property tables do not exist.
2D Flow Area 'NewOrleans Metro' tables complete 16.52 sec
Base: Pipe Network property tables are current.
Completed Writing Geometry
Writing Event Conditions ...
Computing IC for 2D Flow Area: NewOrleans Metro
Completed IC for 2D Flow Area: NewOrleans Metro
Completed Writing Event Condition Data
Geometric Preprocessor HEC-RAS 7.0 April 2026
Finished Processing Geometry
Parallelization CPU Affinity ON (Performance Cores 0,1)
Performing Unsteady Flow Simulation HEC-RAS 7.0 April 2026
Unsteady Input Summary:
2D Area Cells Cores Solvers
NewOrleans Metro 19711 2 SWE-ELM Equation Set (original)
Pipe Network Nodes Cores Courant Target Solvers
Base 64 2 0.90 SWE-ELM Equation Set (original)
Maximum iteration location Cell WSEL ERROR ITERATIONS
Pipe Network Iter Type Cell ERROR Node or Conduit
10APR2024 14:32:28 Base Node 426 0.450 Mine Blvd - Node 1
Overall Volume Accounting Error in Acre Feet: 8.040
Overall Volume Accounting Error as percentage: 0.1223
Please review "Computational Log File" output for volume accounting details
Writing Results to DSS
Finished Unsteady Flow Simulation
Computations Summary
Computation Task Time(hh:mm:ss)
Completing Geometry, Flow and Plan 19
Preprocessing Geometry <1
Unsteady Flow Computations 10:28
Complete Process 10:48
Computation Speed Simulation/Runtime
Unsteady Flow Computations 172x
Complete Process 167x
| condition | timestamps | physical_cells_checked | boundary_cells_zero | absolute_start_utc | absolute_end_utc | expected_total_in | hdf_final_total_in | max_all_time_cell_error_in | native_tolerance_in | hdf_peak_rate_in_hr | absolute_time_qualified | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | baseline | 361 | 19711 | 775 | 2024-04-10 | 2024-04-11 06:00:00 | 0.0 | 0.000000 | 0.000000e+00 | 9.536743e-07 | 0.00 | True |
| 1 | event | 361 | 19711 | 775 | 2024-04-10 | 2024-04-11 06:00:00 | 4.1 | 4.099999 | 5.849202e-07 | 3.910065e-06 | 0.95 | True |
| condition | complete_process | bco_exists | diagnostic_line_count | iteration_diagnostic_lines | volume_accounting_lines | hydraulic_acceptance | |
|---|---|---|---|---|---|---|---|
| 0 | baseline | True | True | 0 | [] | [Overall Volume Accounting Error in Acre Feet:... | not_established |
| 1 | event | True | True | 0 | [Maximum iteration location\tCell\t WSEL\tERRO... | [Overall Volume Accounting Error in Acre Feet:... | not_established |
Python
qualification_summary = build_qualification_summary(
[row for case in native_cases for row in case["qualifications"]],
mode=EXECUTION_MODE, source_module_sha256=SOURCE_SHA256,
library_version=ras_commander.__version__,
)
print("MRMS_ABSOLUTE_TIME_QUALIFICATION=" + json.dumps(qualification_summary, default=str))
(RUN_ROOT / "final_qualification_summary.json").write_text(json.dumps(qualification_summary, default=str, indent=2))
display(pd.DataFrame(qualification_summary["runs"]))
Text Only
MRMS_ABSOLUTE_TIME_QUALIFICATION={"mode": "inspect_retained", "inspection_utc": "2026-09-27T10:13:21.956507Z", "execution_date": "2026-09-27", "source_module_sha256": {"ras_commander/RasUnsteady.py": "9d20063dd04e08256a7c480f24342551e61d90d30e1b5c44bcfd457662645f5d", "ras_commander/precip/PrecipMrms.py": "2fb0df862ed121d3f31f765893cf710222b6a6b0c6195b863bf16733e54c0be4", "ras_commander/RasCmdr.py": "9fd96c94a5165dd69342605b5cf1d9a975d28fd83ddc2d3162ee48cde0fbe50e", "ras_commander/RasPrj.py": "adb984fd94bbfe0bf8cb9ac76816aac8c74a505ecadde1c97364188049cb4176", "ras_commander/hdf/HdfMesh.py": "cb5c56b49b2f5f3a4e039ddd50851a065a5979a595e54cadeffba84d6f1ca93a", "ras_commander/hdf/HdfResultsMesh.py": "44a8cdf45140f05278dc028f9a2489d5db7563973e317bceab941714cf46896d"}, "library_version": "0.103.0", "clock": "UTC-equivalent naive HEC-RAS model labels", "runs": [{"project": "Davis", "hydraulic_acceptance": "not_established", "iteration_diagnostic_lines": [], "volume_accounting_lines": ["Overall Volume Accounting Error in Acre Feet: 0.0000000000", "Overall Volume Accounting Error as percentage: 0.0000000000"], "native_started_epoch": 1790498953.5059812, "native_source_module_sha256": {"ras_commander/RasUnsteady.py": "9d20063dd04e08256a7c480f24342551e61d90d30e1b5c44bcfd457662645f5d", "ras_commander/precip/PrecipMrms.py": "2fb0df862ed121d3f31f765893cf710222b6a6b0c6195b863bf16733e54c0be4", "ras_commander/RasCmdr.py": "9fd96c94a5165dd69342605b5cf1d9a975d28fd83ddc2d3162ee48cde0fbe50e"}, "native_runtime_seconds": 59.33219337463379, "native_hdf_sha256": "2d0211bbd6bfdf845e86cc26972b4118c958bdda770f95f2ee17021bd63c87f2", "native_executable_sha256": "9990c10531221469bf51a9a62ae91f36ec01e651bed83aadf632e678130ae797", "condition": "baseline", "timestamps": 505, "physical_cells_checked": 2705, "boundary_cells_zero": 236, "absolute_start_utc": "2022-12-31 00:00:00", "absolute_end_utc": "2023-01-01 18:00:00", "expected_total_in": 0.0, "hdf_final_total_in": 0.0, "max_all_time_cell_error_in": 0.0, "native_tolerance_in": 9.5367431640625e-07, "hdf_peak_rate_in_hr": 0.0, "absolute_time_qualified": true}, {"project": "Davis", "hydraulic_acceptance": "not_established", "iteration_diagnostic_lines": [], "volume_accounting_lines": ["Overall Volume Accounting Error in Acre Feet: 0.07065", "Overall Volume Accounting Error as percentage: 0.03423"], "native_started_epoch": 1790499012.9724948, "native_source_module_sha256": {"ras_commander/RasUnsteady.py": "9d20063dd04e08256a7c480f24342551e61d90d30e1b5c44bcfd457662645f5d", "ras_commander/precip/PrecipMrms.py": "2fb0df862ed121d3f31f765893cf710222b6a6b0c6195b863bf16733e54c0be4", "ras_commander/RasCmdr.py": "9fd96c94a5165dd69342605b5cf1d9a975d28fd83ddc2d3162ee48cde0fbe50e"}, "native_runtime_seconds": 233.41464376449585, "native_hdf_sha256": "a0fc20df27cd023ac2e879d979ada9178f4e961d7a738c509ff665ea41f622cd", "native_executable_sha256": "9990c10531221469bf51a9a62ae91f36ec01e651bed83aadf632e678130ae797", "condition": "event", "timestamps": 505, "physical_cells_checked": 2705, "boundary_cells_zero": 236, "absolute_start_utc": "2022-12-31 00:00:00", "absolute_end_utc": "2023-01-01 18:00:00", "expected_total_in": 1.7400000000000004, "hdf_final_total_in": 1.7400000095367432, "max_all_time_cell_error_in": 2.0980835002681886e-07, "native_tolerance_in": 1.659393319641822e-06, "hdf_peak_rate_in_hr": 0.20999997854232788, "absolute_time_qualified": true}, {"project": "NewOrleansMetro", "hydraulic_acceptance": "not_established", "iteration_diagnostic_lines": [], "volume_accounting_lines": ["Overall Volume Accounting Error in Acre Feet: 0.000103", "Overall Volume Accounting Error as percentage: 0.0000095671"], "native_started_epoch": 1790499416.668967, "native_source_module_sha256": {"ras_commander/RasUnsteady.py": "9d20063dd04e08256a7c480f24342551e61d90d30e1b5c44bcfd457662645f5d", "ras_commander/precip/PrecipMrms.py": "2fb0df862ed121d3f31f765893cf710222b6a6b0c6195b863bf16733e54c0be4", "ras_commander/RasCmdr.py": "9fd96c94a5165dd69342605b5cf1d9a975d28fd83ddc2d3162ee48cde0fbe50e"}, "native_runtime_seconds": 536.6128127574921, "native_hdf_sha256": "2187249cc9e5de9429f83d07eebd5f7606d6fed943dc4506067a225fdd58a3c0", "native_executable_sha256": "9990c10531221469bf51a9a62ae91f36ec01e651bed83aadf632e678130ae797", "condition": "baseline", "timestamps": 361, "physical_cells_checked": 19711, "boundary_cells_zero": 775, "absolute_start_utc": "2024-04-10 00:00:00", "absolute_end_utc": "2024-04-11 06:00:00", "expected_total_in": 0.0, "hdf_final_total_in": 0.0, "max_all_time_cell_error_in": 0.0, "native_tolerance_in": 9.5367431640625e-07, "hdf_peak_rate_in_hr": 0.0, "absolute_time_qualified": true}, {"project": "NewOrleansMetro", "hydraulic_acceptance": "not_established", "iteration_diagnostic_lines": ["Maximum iteration location\tCell\t WSEL\tERROR\tITERATIONS", "Pipe Network Iter\tType\tCell\tERROR\tNode or Conduit", "10APR2024 14:32:28 Base\tNode\t 426\t 0.450\tMine Blvd - Node 1"], "volume_accounting_lines": ["Overall Volume Accounting Error in Acre Feet: 8.040", "Overall Volume Accounting Error as percentage: 0.1223"], "native_started_epoch": 1790500608.2070925, "native_source_module_sha256": {"ras_commander/RasUnsteady.py": "9d20063dd04e08256a7c480f24342551e61d90d30e1b5c44bcfd457662645f5d", "ras_commander/precip/PrecipMrms.py": "2fb0df862ed121d3f31f765893cf710222b6a6b0c6195b863bf16733e54c0be4", "ras_commander/RasCmdr.py": "9fd96c94a5165dd69342605b5cf1d9a975d28fd83ddc2d3162ee48cde0fbe50e"}, "native_runtime_seconds": 651.8917608261108, "native_hdf_sha256": "945e56ec76459dfc3067fb9689e291f2270351778e4ec13ae56f31be27cc8a57", "native_executable_sha256": "9990c10531221469bf51a9a62ae91f36ec01e651bed83aadf632e678130ae797", "condition": "event", "timestamps": 361, "physical_cells_checked": 19711, "boundary_cells_zero": 775, "absolute_start_utc": "2024-04-10 00:00:00", "absolute_end_utc": "2024-04-11 06:00:00", "expected_total_in": 4.1, "hdf_final_total_in": 4.09999942779541, "max_all_time_cell_error_in": 5.849202473484638e-07, "native_tolerance_in": 3.910064606316155e-06, "hdf_peak_rate_in_hr": 0.9499999284744263, "absolute_time_qualified": true}], "native_runtime_seconds_total": 1481.2514107227325}
| project | hydraulic_acceptance | iteration_diagnostic_lines | volume_accounting_lines | native_started_epoch | native_source_module_sha256 | native_runtime_seconds | native_hdf_sha256 | native_executable_sha256 | condition | ... | physical_cells_checked | boundary_cells_zero | absolute_start_utc | absolute_end_utc | expected_total_in | hdf_final_total_in | max_all_time_cell_error_in | native_tolerance_in | hdf_peak_rate_in_hr | absolute_time_qualified | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Davis | not_established | [] | [Overall Volume Accounting Error in Acre Feet:... | 1.790499e+09 | {'ras_commander/RasUnsteady.py': '9d20063dd04e... | 59.332193 | 2d0211bbd6bfdf845e86cc26972b4118c958bdda770f95... | 9990c10531221469bf51a9a62ae91f36ec01e651bed83a... | baseline | ... | 2705 | 236 | 2022-12-31 | 2023-01-01 18:00:00 | 0.00 | 0.000000 | 0.000000e+00 | 9.536743e-07 | 0.00 | True |
| 1 | Davis | not_established | [] | [Overall Volume Accounting Error in Acre Feet:... | 1.790499e+09 | {'ras_commander/RasUnsteady.py': '9d20063dd04e... | 233.414644 | a0fc20df27cd023ac2e879d979ada9178f4e961d7a738c... | 9990c10531221469bf51a9a62ae91f36ec01e651bed83a... | event | ... | 2705 | 236 | 2022-12-31 | 2023-01-01 18:00:00 | 1.74 | 1.740000 | 2.098084e-07 | 1.659393e-06 | 0.21 | True |
| 2 | NewOrleansMetro | not_established | [] | [Overall Volume Accounting Error in Acre Feet:... | 1.790499e+09 | {'ras_commander/RasUnsteady.py': '9d20063dd04e... | 536.612813 | 2187249cc9e5de9429f83d07eebd5f7606d6fed943dc45... | 9990c10531221469bf51a9a62ae91f36ec01e651bed83a... | baseline | ... | 19711 | 775 | 2024-04-10 | 2024-04-11 06:00:00 | 0.00 | 0.000000 | 0.000000e+00 | 9.536743e-07 | 0.00 | True |
| 3 | NewOrleansMetro | not_established | [Maximum iteration location\tCell\t WSEL\tERRO... | [Overall Volume Accounting Error in Acre Feet:... | 1.790501e+09 | {'ras_commander/RasUnsteady.py': '9d20063dd04e... | 651.891761 | 945e56ec76459dfc3067fb9689e291f2270351778e4ec1... | 9990c10531221469bf51a9a62ae91f36ec01e651bed83a... | event | ... | 19711 | 775 | 2024-04-10 | 2024-04-11 06:00:00 | 4.10 | 4.099999 | 5.849202e-07 | 3.910065e-06 | 0.95 | True |
4 rows × 21 columns
Consistent diagnostic visualization without network basemaps¶
All final depth frames use the same public HDF nearest-cell export at the stated resolution. They are diagnostic rasterizations, not native Mapper pixel/extent equivalents or calibrated flood maps. The earlier Davis Mapper attempt exposed TIFF decode/nodata warnings and mixed interpolation footprints; its files and outputs are preserved separately. The original online map tiles were blocked, so this execution makes no basemap requests and uses coordinate grids, mesh outlines, and pump overlays. Observed rainfall bars occupy their preceding hourly intervals; dry-tail zeros remain a declared assumption.
Flood videos omit the optional terrain backdrop: the animation API stretches a raw terrain array to the flood extent, and these examples have different terrain/flood extents. Prior renders are archived separately. Final geographic static axes use four tick intervals; combined videos use scoped Matplotlib coordinate offsets for readability. All full-resolution depth rasters/statistics are retained.
Terrain-overlay follow-up: issue #441. The final videos omit terrain; no display downsampling was used.
Text Only
Rendering previously qualified native result; no solver execution in this cell.
Plan HDF: working/917_alignment_20260927/m/Davis/P/event_retained.p02.hdf
Total output timesteps: 505
Available consecutive 5-minute timesteps: 505
Selected map and animation frames: 30
Selected frame coverage: 31DEC2022 00:00:00 to 01JAN2023 18:00:00
Verified raster variable: no native Depth dataset; exporter uses max(Water Surface - Cells Minimum Elevation, 0), in feet.
2026-09-27 06:13:23 - ras_commander.hdf.HdfResultsMesh - INFO - Wrote 30 depth raster(s) for mesh 'area2'
Consistent HDF nearest-cell invert-depth rasters (not native Mapper equivalence)
Showing first and last 3 of 30 rows
| timestamp | source_rasters | source | |
|---|---|---|---|
| 0 | 2022-12-31 00:00:00 | 1 | HDF nearest-cell invert-depth raster |
| 1 | 2022-12-31 01:25:00 | 1 | HDF nearest-cell invert-depth raster |
| 2 | 2022-12-31 02:50:00 | 1 | HDF nearest-cell invert-depth raster |
| 27 | 2023-01-01 15:05:00 | 1 | HDF nearest-cell invert-depth raster |
| 28 | 2023-01-01 16:30:00 | 1 | HDF nearest-cell invert-depth raster |
| 29 | 2023-01-01 18:00:00 | 1 | HDF nearest-cell invert-depth raster |
Text Only
Showing first and last 3 of 30 rows
| timestamp | raster | width | height | cell_size_ft | size_mb | |
|---|---|---|---|---|---|---|
| 0 | 2022-12-31 00:00:00 | depth_20221231_000000.tif | 393 | 302 | 28.845245 | 0.009795 |
| 1 | 2022-12-31 01:25:00 | depth_20221231_012500.tif | 393 | 302 | 28.845245 | 0.036023 |
| 2 | 2022-12-31 02:50:00 | depth_20221231_025000.tif | 393 | 302 | 28.845245 | 0.037631 |
| 27 | 2023-01-01 15:05:00 | depth_20230101_150500.tif | 393 | 302 | 28.845245 | 0.040511 |
| 28 | 2023-01-01 16:30:00 | depth_20230101_163000.tif | 393 | 302 | 28.845245 | 0.040403 |
| 29 | 2023-01-01 18:00:00 | depth_20230101_180000.tif | 393 | 302 | 28.845245 | 0.040363 |
Text Only
Showing first and last 3 of 30 rows
| timestamp | wet_cells | wet_area_acres | raster_invert_max_ft | raster_invert_mean_wet_ft | native_hydraulic_max_ft | native_invert_max_ft | raster_wetness_matches_hydraulic | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2022-12-31 00:00:00 | 0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | True |
| 1 | 2022-12-31 01:25:00 | 72392 | 1382.773872 | 0.520729 | 0.102102 | 0.189000 | 0.520729 | True |
| 2 | 2022-12-31 02:50:00 | 72392 | 1382.773872 | 0.781048 | 0.154556 | 0.348944 | 0.781048 | True |
| 27 | 2023-01-01 15:05:00 | 72392 | 1382.773872 | 3.515953 | 0.228310 | 2.008721 | 3.515953 | True |
| 28 | 2023-01-01 16:30:00 | 72392 | 1382.773872 | 3.516327 | 0.214692 | 2.009027 | 3.516327 | True |
| 29 | 2023-01-01 18:00:00 | 72392 | 1382.773872 | 3.516560 | 0.200669 | 2.009217 | 3.516560 | True |
Text Only
Consistent HDF nearest-cell invert-depth rasters (not native Mapper equivalence): 30 source timesteps, 30 HDF source frames written to 30 single 28.845-ft rasters
Diagnostic raster folder: working/917_alignment_20260927/m/Davis/HdfConsolidated
HDF rasters per selected timestep: min=1, median=1, max=1
Static invert-depth maps use the actual raster maximum. Animation APIs cap colors at the 98th percentile of positive values (minimum 0.1); this is a display cap, not maximum invert depth.
2026-09-27 06:13:24 - matplotlib.animation - INFO - Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2026-09-27 06:13:24 - matplotlib.animation - INFO - MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 880x660 -pix_fmt rgba -framerate 4 -loglevel error -i pipe: -vcodec h264 -pix_fmt yuv420p -b 1800k -threads 2 -y 'C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\videos_hdf\mrms_precipitation_davis_ar.mp4'
2026-09-27 06:13:34 - matplotlib.animation - INFO - Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2026-09-27 06:13:34 - matplotlib.animation - INFO - MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 880x660 -pix_fmt rgba -framerate 4 -loglevel error -i pipe: -vcodec h264 -pix_fmt yuv420p -b 1800k -threads 2 -y 'C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\videos_hdf\flood_inundation_depth_davis_ar_5min.mp4'
Showing first and last 3 of 30 rows
| flood_time | source_interval_start | source_interval_end | display_mode | source_index | |
|---|---|---|---|---|---|
| 0 | 2022-12-31 00:00:00 | NaT | NaT | initial_zero | NaN |
| 1 | 2022-12-31 01:25:00 | 2022-12-31 01:00:00 | 2022-12-31 02:00:00 | observed_interval | 1.0 |
| 2 | 2022-12-31 02:50:00 | 2022-12-31 02:00:00 | 2022-12-31 03:00:00 | observed_interval | 2.0 |
| 27 | 2023-01-01 15:05:00 | NaT | NaT | dry_tail | NaN |
| 28 | 2023-01-01 16:30:00 | NaT | NaT | dry_tail | NaN |
| 29 | 2023-01-01 18:00:00 | NaT | NaT | dry_tail | NaN |
Text Only
Combined panel: hourly interval covering flood time; exact endpoint uses the preceding interval; initial and dry-tail frames zero. Spatial source illustration, uniform model boundary.
2026-09-27 06:13:44 - matplotlib.animation - INFO - Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2026-09-27 06:13:44 - matplotlib.animation - INFO - MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 1320x550 -pix_fmt rgba -framerate 4 -loglevel error -i pipe: -vcodec h264 -pix_fmt yuv420p -b 1800k -threads 2 -y 'C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\videos_hdf\combined_precip_flood_davis_ar_5min.mp4'
Video: working/917_alignment_20260927/videos_hdf/mrms_precipitation_davis_ar.mp4 (416,652 bytes)
Video: working/917_alignment_20260927/videos_hdf/flood_inundation_depth_davis_ar_5min.mp4 (977,482 bytes)
Video: working/917_alignment_20260927/videos_hdf/combined_precip_flood_davis_ar_5min.mp4 (907,127 bytes)
| mesh_name | reference_cell_id | selection_metric | peak_response_time | baseline_hydraulic_depth_ft | event_hydraulic_depth_ft | event_minus_baseline_hydraulic_depth_ft | |
|---|---|---|---|---|---|---|---|
| 0 | area2 | 600 | maximum concurrent event-minus-baseline depth | 2022-12-31 21:15:00 | 0.0 | 2.571257 | 2.571257 |
| case | mesh_name | reference_cell_id | baseline_peak_hydraulic_depth_ft | event_peak_hydraulic_depth_ft | peak_hydraulic_depth_delta_ft | baseline_peak_wse_ft | event_peak_wse_ft | peak_wse_delta_ft | event_peak_precip_in_hr | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | davis_ar | area2 | 600 | 0.0 | 2.571257 | 2.571257 | 37.542934 | 41.224609 | 3.681675 | 0.21 |
| station | baseline_peak_flow_cfs | event_peak_flow_cfs | peak_flow_delta_cfs | event_max_pumps_on | first_on | last_on | |
|---|---|---|---|---|---|---|---|
| 0 | Pump Station #1 | 0.0 | 67.372246 | 67.372246 | 1.0 | 2022-12-31 03:30:00 | 2023-01-01 18:00:00 |
Text Only
MP4 animations generated for Davis atmospheric river:
| animation | file | size_mb | |
|---|---|---|---|
| 0 | Precipitation only | mrms_precipitation_davis_ar.mp4 | 0.42 |
| 1 | Flood inundation only | flood_inundation_depth_davis_ar_5min.mp4 | 0.98 |
| 2 | Combined precipitation and flood | combined_precip_flood_davis_ar_5min.mp4 | 0.91 |
Text Only
Rendering previously qualified native result; no solver execution in this cell.
Plan HDF: working/917_alignment_20260927/m/NewOrleansMetro/P/event_retained.p01.hdf
Total output timesteps: 361
Available consecutive 5-minute timesteps: 361
Selected map and animation frames: 30
Selected frame coverage: 10APR2024 00:00:00 to 11APR2024 06:00:00
Verified raster variable: no native Depth dataset; exporter uses max(Water Surface - Cells Minimum Elevation, 0), in feet.
2026-09-27 06:14:21 - ras_commander.hdf.HdfResultsMesh - INFO - Wrote 30 depth raster(s) for mesh 'NewOrleans Metro'
Consistent HDF nearest-cell invert-depth rasters (not native Mapper equivalence)
Showing first and last 3 of 30 rows
| timestamp | source_rasters | source | |
|---|---|---|---|
| 0 | 2024-04-10 00:00:00 | 1 | HDF nearest-cell invert-depth raster |
| 1 | 2024-04-10 01:00:00 | 1 | HDF nearest-cell invert-depth raster |
| 2 | 2024-04-10 02:00:00 | 1 | HDF nearest-cell invert-depth raster |
| 27 | 2024-04-11 03:55:00 | 1 | HDF nearest-cell invert-depth raster |
| 28 | 2024-04-11 04:55:00 | 1 | HDF nearest-cell invert-depth raster |
| 29 | 2024-04-11 06:00:00 | 1 | HDF nearest-cell invert-depth raster |
Text Only
Showing first and last 3 of 30 rows
| timestamp | raster | width | height | cell_size_ft | size_mb | |
|---|---|---|---|---|---|---|
| 0 | 2024-04-10 00:00:00 | depth_20240410_000000.tif | 3413 | 4086 | 10.0 | 0.492137 |
| 1 | 2024-04-10 01:00:00 | depth_20240410_010000.tif | 3413 | 4086 | 10.0 | 0.493851 |
| 2 | 2024-04-10 02:00:00 | depth_20240410_020000.tif | 3413 | 4086 | 10.0 | 0.490827 |
| 27 | 2024-04-11 03:55:00 | depth_20240411_035500.tif | 3413 | 4086 | 10.0 | 3.498587 |
| 28 | 2024-04-11 04:55:00 | depth_20240411_045500.tif | 3413 | 4086 | 10.0 | 3.500441 |
| 29 | 2024-04-11 06:00:00 | depth_20240411_060000.tif | 3413 | 4086 | 10.0 | 3.505505 |
Text Only
Showing first and last 3 of 30 rows
| timestamp | wet_cells | wet_area_acres | raster_invert_max_ft | raster_invert_mean_wet_ft | native_hydraulic_max_ft | native_invert_max_ft | raster_wetness_matches_hydraulic | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2024-04-10 00:00:00 | 46363 | 106.434803 | 44.527458 | 13.098394 | 22.845560 | 44.527458 | True |
| 1 | 2024-04-10 01:00:00 | 46363 | 106.434803 | 44.527355 | 13.097767 | 22.845457 | 44.527355 | True |
| 2 | 2024-04-10 02:00:00 | 46363 | 106.434803 | 44.526001 | 13.097028 | 22.844105 | 44.526001 | True |
| 27 | 2024-04-11 03:55:00 | 6704897 | 15392.325528 | 45.480717 | 1.046657 | 23.797974 | 45.480717 | True |
| 28 | 2024-04-11 04:55:00 | 6704897 | 15392.325528 | 45.312046 | 1.019962 | 23.629322 | 45.312046 | True |
| 29 | 2024-04-11 06:00:00 | 6704897 | 15392.325528 | 45.174095 | 0.992680 | 23.491385 | 45.174095 | True |
Text Only
Consistent HDF nearest-cell invert-depth rasters (not native Mapper equivalence): 30 source timesteps, 30 HDF source frames written to 30 single 10.000-ft rasters
Diagnostic raster folder: working/917_alignment_20260927/m/NewOrleansMetro/HdfConsolidated
HDF rasters per selected timestep: min=1, median=1, max=1
Static invert-depth maps use the actual raster maximum. Animation APIs cap colors at the 98th percentile of positive values (minimum 0.1); this is a display cap, not maximum invert depth.
2026-09-27 06:14:49 - matplotlib.animation - INFO - Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2026-09-27 06:14:49 - matplotlib.animation - INFO - MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 880x660 -pix_fmt rgba -framerate 4 -loglevel error -i pipe: -vcodec h264 -pix_fmt yuv420p -b 1800k -threads 2 -y 'C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\videos_hdf\mrms_precipitation_nola_apr24.mp4'
2026-09-27 06:15:01 - matplotlib.animation - INFO - Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2026-09-27 06:15:01 - matplotlib.animation - INFO - MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 880x660 -pix_fmt rgba -framerate 4 -loglevel error -i pipe: -vcodec h264 -pix_fmt yuv420p -b 1800k -threads 2 -y 'C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\videos_hdf\flood_inundation_depth_nola_apr24_5min.mp4'
Showing first and last 3 of 30 rows
| flood_time | source_interval_start | source_interval_end | display_mode | source_index | |
|---|---|---|---|---|---|
| 0 | 2024-04-10 00:00:00 | NaT | NaT | initial_zero | NaN |
| 1 | 2024-04-10 01:00:00 | 2024-04-10 00:00:00 | 2024-04-10 01:00:00 | observed_interval | 0.0 |
| 2 | 2024-04-10 02:00:00 | 2024-04-10 01:00:00 | 2024-04-10 02:00:00 | observed_interval | 1.0 |
| 27 | 2024-04-11 03:55:00 | NaT | NaT | dry_tail | NaN |
| 28 | 2024-04-11 04:55:00 | NaT | NaT | dry_tail | NaN |
| 29 | 2024-04-11 06:00:00 | NaT | NaT | dry_tail | NaN |
Text Only
Combined panel: hourly interval covering flood time; exact endpoint uses the preceding interval; initial and dry-tail frames zero. Spatial source illustration, uniform model boundary.
2026-09-27 06:16:19 - matplotlib.animation - INFO - Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2026-09-27 06:16:19 - matplotlib.animation - INFO - MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 1320x550 -pix_fmt rgba -framerate 4 -loglevel error -i pipe: -vcodec h264 -pix_fmt yuv420p -b 1800k -threads 2 -y 'C:\Users\billk_clb\.codex\workspaces\ras-docs-priority-20260926\working\917_alignment_20260927\videos_hdf\combined_precip_flood_nola_apr24_5min.mp4'
Video: working/917_alignment_20260927/videos_hdf/mrms_precipitation_nola_apr24.mp4 (370,914 bytes)
Video: working/917_alignment_20260927/videos_hdf/flood_inundation_depth_nola_apr24_5min.mp4 (867,675 bytes)
Video: working/917_alignment_20260927/videos_hdf/combined_precip_flood_nola_apr24_5min.mp4 (701,362 bytes)
| mesh_name | reference_cell_id | selection_metric | peak_response_time | baseline_hydraulic_depth_ft | event_hydraulic_depth_ft | event_minus_baseline_hydraulic_depth_ft | |
|---|---|---|---|---|---|---|---|
| 0 | NewOrleans Metro | 17519 | maximum concurrent event-minus-baseline depth | 2024-04-10 14:45:00 | 0.0 | 20.142014 | 20.142014 |
| case | mesh_name | reference_cell_id | baseline_peak_hydraulic_depth_ft | event_peak_hydraulic_depth_ft | peak_hydraulic_depth_delta_ft | baseline_peak_wse_ft | event_peak_wse_ft | peak_wse_delta_ft | event_peak_precip_in_hr | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | nola_apr24 | NewOrleans Metro | 17519 | 0.0 | 20.142014 | 20.142014 | -30.0 | -9.029428 | 20.970572 | 0.95 |
| station | baseline_peak_flow_cfs | event_peak_flow_cfs | peak_flow_delta_cfs | event_max_pumps_on | first_on | last_on | |
|---|---|---|---|---|---|---|---|
| 0 | 17th St Pumps | 0.0 | 2594.920410 | 2594.920410 | 4.0 | 2024-04-10 15:40:00 | 2024-04-11 06:00:00 |
| 1 | DPS 12 | 0.0 | 265.534302 | 265.534302 | 1.0 | 2024-04-10 14:45:00 | 2024-04-11 06:00:00 |
Text Only
MP4 animations generated for NewOrleansMetro April 2024 flash flood:
| animation | file | size_mb | |
|---|---|---|---|
| 0 | Precipitation only | mrms_precipitation_nola_apr24.mp4 | 0.37 |
| 1 | Flood inundation only | flood_inundation_depth_nola_apr24_5min.mp4 | 0.87 |
| 2 | Combined precipitation and flood | combined_precip_flood_nola_apr24_5min.mp4 | 0.70 |
Python
summary = pd.DataFrame(
[
{
"case": result["case"],
"mesh_count": result["mesh_count"],
"pump_count": result["pump_count"],
"pump_operation_stations": len(result["event_pump_operation"]),
"reference_cell_id": result["reference_cell_id"],
"peak_hydraulic_depth_delta_ft": result["reference_summary"].iloc[0]["peak_hydraulic_depth_delta_ft"],
"peak_wse_delta_ft": result["reference_summary"].iloc[0]["peak_wse_delta_ft"],
"total_output_timesteps": result["total_output_timesteps"],
"available_animation_timesteps": result["available_animation_timesteps"],
"animation_timesteps": result["animation_timesteps"],
"animation_frames": result["animation_frames"],
"source_raster_frames": result["depth_tile_count"],
"diagnostic_rasters": result["consolidated_raster_count"],
"dominant_resolution_ft": result["dominant_terrain_resolution_ft"],
"consolidated_resolution_ft": result["consolidation_cell_size_ft"],
"dss_file": result["mrms_dss"].name,
"videos": len(result["videos"]),
}
for result in results
]
)
display(summary)
| case | mesh_count | pump_count | pump_operation_stations | reference_cell_id | peak_hydraulic_depth_delta_ft | peak_wse_delta_ft | total_output_timesteps | available_animation_timesteps | animation_timesteps | animation_frames | source_raster_frames | diagnostic_rasters | dominant_resolution_ft | consolidated_resolution_ft | dss_file | videos | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | davis_ar | 2 | 1 | 1 | 600 | 2.571257 | 3.681675 | 505 | 505 | 30 | 30 | 30 | 30 | 28.845245 | 28.845245 | davis_ar_mrms_qpe.dss | 3 |
| 1 | nola_apr24 | 1 | 2 | 2 | 17519 | 20.142014 | 20.970572 | 361 | 361 | 30 | 30 | 30 | 30 | 3.280000 | 10.000000 | nola_apr24_mrms_qpe.dss | 3 |
Reviewer Artifact Summary¶
Python
print("Final reviewer artifacts generated by this notebook run:")
for result in results:
print(f"- {result['case']}")
print(f" Project folder: {display_path(result['project_folder'])}")
print(f" Plan HDF: {display_path(result['plan_hdf'])} ({result['plan_hdf'].stat().st_size:,} bytes)")
print(
f" No-rain baseline HDF: {display_path(result['baseline_plan_hdf'])} "
f"({result['baseline_plan_hdf'].stat().st_size:,} bytes)"
)
print(f" MRMS DSS: {display_path(result['mrms_dss'])} ({result['mrms_dss'].stat().st_size:,} bytes)")
print(
f" Selected HDF time frames: {len(result['depth_tile_frames'])}; "
f"source rasters: {result['depth_tile_count']}; "
f"diagnostic rasters: {result['consolidated_raster_count']}; "
f"cell size: {result['consolidation_cell_size_ft']:.3f} ft"
)
print(f" Diagnostic raster folder: {display_path(result['consolidated_depth_rasters'][0].parent)}")
print(f" Reference cell: {result['reference_cell_id']}")
print(f" Pump operation stations plotted: {len(result['event_pump_operation'])}")
for video in result["videos"]:
print(f" Video artifact: {display_path(video)} ({video.stat().st_size:,} bytes)")
print("\nQA/QC reminder: these runs use an intentionally uniform precipitation boundary.")
print("Review all BCO/runtime diagnostics, precipitation maps and time series, terrain coverage, and hydraulic results manually.")
print("Use notebook 924 when spatial MRMS variability must be consumed by the solver.")
Text Only
Final reviewer artifacts generated by this notebook run:
- davis_ar
Project folder: working/917_alignment_20260927/m/Davis
Plan HDF: working/917_alignment_20260927/m/Davis/P/event_retained.p02.hdf (50,621,942 bytes)
No-rain baseline HDF: working/917_alignment_20260927/m/Davis/P/davis_ar_no_rain_baseline.p02.hdf (3,595,465 bytes)
MRMS DSS: working/917_alignment_20260927/m/Davis/P/davis_ar_mrms_qpe.dss (220,608 bytes)
Selected HDF time frames: 30; source rasters: 30; diagnostic rasters: 30; cell size: 28.845 ft
Diagnostic raster folder: working/917_alignment_20260927/m/Davis/HdfConsolidated
Reference cell: 600
Pump operation stations plotted: 1
Video artifact: working/917_alignment_20260927/videos_hdf/mrms_precipitation_davis_ar.mp4 (416,652 bytes)
Video artifact: working/917_alignment_20260927/videos_hdf/flood_inundation_depth_davis_ar_5min.mp4 (977,482 bytes)
Video artifact: working/917_alignment_20260927/videos_hdf/combined_precip_flood_davis_ar_5min.mp4 (907,127 bytes)
- nola_apr24
Project folder: working/917_alignment_20260927/m/NewOrleansMetro
Plan HDF: working/917_alignment_20260927/m/NewOrleansMetro/P/event_retained.p01.hdf (105,503,626 bytes)
No-rain baseline HDF: working/917_alignment_20260927/m/NewOrleansMetro/P/nola_apr24_no_rain_baseline.p01.hdf (40,749,856 bytes)
MRMS DSS: working/917_alignment_20260927/m/NewOrleansMetro/P/nola_apr24_mrms_qpe.dss (155,648 bytes)
Selected HDF time frames: 30; source rasters: 30; diagnostic rasters: 30; cell size: 10.000 ft
Diagnostic raster folder: working/917_alignment_20260927/m/NewOrleansMetro/HdfConsolidated
Reference cell: 17519
Pump operation stations plotted: 2
Video artifact: working/917_alignment_20260927/videos_hdf/mrms_precipitation_nola_apr24.mp4 (370,914 bytes)
Video artifact: working/917_alignment_20260927/videos_hdf/flood_inundation_depth_nola_apr24_5min.mp4 (867,675 bytes)
Video artifact: working/917_alignment_20260927/videos_hdf/combined_precip_flood_nola_apr24_5min.mp4 (701,362 bytes)
QA/QC reminder: these runs use an intentionally uniform precipitation boundary.
Review all BCO/runtime diagnostics, precipitation maps and time series, terrain coverage, and hydraulic results manually.
Use notebook 924 when spatial MRMS variability must be consumed by the solver.