MRMS NetCDF Rain-on-Grid Validation

This notebook validates the direct MRMS gridded precipitation path in ras-commander: MRMS GRIB2 download, NetCDF export, RasUnsteady.set_gridded_precipitation() HDF import, HEC-RAS compute, and model-output precipitation checks on a RasExamples project.

Python

# =============================================================================
# DEVELOPMENT MODE TOGGLE
# =============================================================================
# True uses the local checkout; False uses the installed ras-commander package.

USE_LOCAL_SOURCE = True

import sys
from pathlib import Path

REPO_ROOT = Path.cwd()
if REPO_ROOT.name.lower() == "examples":
    REPO_ROOT = REPO_ROOT.parent

if USE_LOCAL_SOURCE:
    local_path = str(REPO_ROOT)
    if local_path not in sys.path:
        sys.path.insert(0, local_path)
    print(f"LOCAL SOURCE MODE: Loading from {REPO_ROOT / 'ras_commander'}")
else:
    print("PIP PACKAGE MODE: Loading installed ras-commander")

Python

from __future__ import annotations

from datetime import datetime
import os
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from ras_commander import RasCmdr, RasExamples, RasPlan, RasUnsteady, init_ras_project
from ras_commander.hdf import HdfMesh, HdfProject, HdfResultsMesh
from ras_commander.precip import PrecipMrms

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" / "example_924_mrms_netcdf_rain_on_grid",
    )
)
RUN_ROOT.mkdir(parents=True, exist_ok=True)

PROJECT_NAME = "BaldEagleCrkMulti2D"
PROJECT_SUFFIX = "mrms_netcdf_924"
RAS_VERSION = "7.0"
EVENT_START = datetime(2024, 8, 9, 10)
EVENT_LAST_QPE = datetime(2024, 8, 9, 13)
SIM_END = datetime(2024, 8, 9, 16)
MRMS_RESOLUTION_METERS = 8000.0

print(f"Run root: {RUN_ROOT}")
print(f"MRMS event: {EVENT_START} through {EVENT_LAST_QPE}")

Python

def find_gridded_precipitation_plan(ras):
    """Return the first plan already configured for gridded precipitation."""
    for _, row in ras.plan_df.iterrows():
        unsteady_number = str(row["unsteady_number"]).zfill(2)
        config = RasUnsteady.get_met_precipitation_config(unsteady_number, ras_object=ras)
        if config.get("enabled") and config.get("mode") == "Gridded":
            return row, config
    raise RuntimeError("No enabled gridded precipitation plan found in RasExamples project")


def get_plan_geometry_hdf(ras, plan_row: pd.Series) -> Path:
    geometry_number = str(plan_row["geometry_number"]).zfill(2)
    geom_hdf = ras.project_folder / f"{ras.project_name}.g{geometry_number}.hdf"
    if not geom_hdf.exists():
        raise FileNotFoundError(f"Geometry HDF not found: {geom_hdf}")
    return geom_hdf


def mesh_precip_summary(plan_hdf: Path, variable: str) -> pd.DataFrame:
    mesh_data = HdfResultsMesh.get_mesh_cells_timeseries(plan_hdf, var=variable)
    records = []
    for mesh_name, dataset in mesh_data.items():
        values = dataset[variable].values
        final_values = values[-1]
        records.append(
            {
                "mesh": mesh_name,
                "timesteps": int(values.shape[0]),
                "cells": int(values.shape[1]),
                "min": float(np.nanmin(values)),
                "mean_final": float(np.nanmean(final_values)),
                "max": float(np.nanmax(values)),
                "positive_final_cells": int(np.count_nonzero(final_values > 0)),
            }
        )
    return pd.DataFrame.from_records(records)

Extract the Example Project

Python

project_path = RasExamples.extract_project(
    PROJECT_NAME,
    output_path=RUN_ROOT,
    suffix=PROJECT_SUFFIX,
)
ras = init_ras_project(
    project_path,
    ras_version=RAS_VERSION,
    load_results_summary=False,
)

plan_row, starting_precip_config = find_gridded_precipitation_plan(ras)
PLAN_NUMBER = str(plan_row["plan_number"]).zfill(2)
UNSTEADY_NUMBER = str(plan_row["unsteady_number"]).zfill(2)
GEOMETRY_NUMBER = str(plan_row["geometry_number"]).zfill(2)
GEOMETRY_HDF = get_plan_geometry_hdf(ras, plan_row)
BOUNDS = tuple(
    float(value)
    for value in HdfProject.get_project_bounds_latlon(GEOMETRY_HDF, buffer_percent=10.0)
)

print(f"Project folder: {ras.project_folder}")
print(f"Selected plan: p{PLAN_NUMBER}; unsteady file: u{UNSTEADY_NUMBER}; geometry: g{GEOMETRY_NUMBER}")
print(f"Starting precipitation source: {starting_precip_config.get('source')}")
print(f"Project bounds with 10% buffer: {BOUNDS}")
ras.plan_df[["plan_number", "geometry_number", "unsteady_number"]]

Model Context

Python

mesh_areas = HdfMesh.get_mesh_areas(GEOMETRY_HDF)
mesh_summary = pd.DataFrame(
    {
        "project": [PROJECT_NAME],
        "plan": [f"p{PLAN_NUMBER}"],
        "geometry": [f"g{GEOMETRY_NUMBER}"],
        "unsteady": [f"u{UNSTEADY_NUMBER}"],
        "starting_precip_source": [starting_precip_config.get("source")],
        "mesh_areas": [", ".join(mesh_areas["mesh_name"].astype(str).tolist())],
        "simulation_start": [EVENT_START],
        "simulation_end": [SIM_END],
    }
)
display(mesh_summary)

fig, ax = plt.subplots(figsize=(7, 5))
plot_mesh = mesh_areas.to_crs("EPSG:4326") if mesh_areas.crs else mesh_areas
plot_mesh.boundary.plot(ax=ax, color="#1f4e79", linewidth=1.5, label="2D mesh boundary")
west, south, east, north = BOUNDS
ax.plot(
    [west, east, east, west, west],
    [south, south, north, north, south],
    color="#b84a39",
    linewidth=1.4,
    linestyle="--",
    label="MRMS download bounds",
)
ax.set_title("Bald Eagle RasExamples model and MRMS query extent")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.grid(True, alpha=0.25)
ax.legend(loc="best")
plt.show()

Download and Inspect MRMS QPE

Python

catalog = PrecipMrms.catalog(
    bounds=BOUNDS,
    start_date=EVENT_START,
    end_date=EVENT_LAST_QPE,
    source="noaa_s3",
)
expected_files = int((EVENT_LAST_QPE - EVENT_START).total_seconds() // 3600) + 1
assert len(catalog) >= expected_files, f"Expected {expected_files} MRMS frames, found {len(catalog)}"

mrms_files = PrecipMrms.download(
    bounds=BOUNDS,
    start_time=EVENT_START,
    end_time=EVENT_LAST_QPE,
    output_dir=ras.project_folder / "P" / "mrms" / "grib2",
    source="noaa_s3",
)
assert len(mrms_files) >= expected_files

mrms_stack = PrecipMrms.load_grib2_stack(mrms_files, bounds=BOUNDS)
hyetograph = PrecipMrms.to_hyetograph(mrms_stack)
assert float(hyetograph["cumulative_depth"].iloc[-1]) > 0.0

print(f"Downloaded {len(mrms_files)} MRMS GRIB2 files")
print(f"MRMS stack shape: {mrms_stack.shape}")
print(f"Spatial-average event depth: {hyetograph['cumulative_depth'].iloc[-1]:.3f} inches")
display(catalog[["valid_time", "archive_product", "filename", "compressed"]])
display(hyetograph)

Python

fig, ax1 = plt.subplots(figsize=(9, 4))
ax1.bar(
    hyetograph["time"],
    hyetograph["incremental_depth"],
    width=0.03,
    color="#2f6f9f",
    label="Incremental depth",
)
ax1.set_ylabel("Incremental depth (in)")
ax1.set_xlabel("MRMS valid time")
ax1.grid(True, axis="y", alpha=0.25)

ax2 = ax1.twinx()
ax2.plot(
    hyetograph["time"],
    hyetograph["cumulative_depth"],
    color="#b84a39",
    marker="o",
    linewidth=2,
    label="Cumulative depth",
)
ax2.set_ylabel("Cumulative depth (in)")

handles1, labels1 = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles1 + handles2, labels1 + labels2, loc="upper left")
ax1.set_title("MRMS QPE over Bald Eagle example-project bounds")
fig.autofmt_xdate()
plt.show()

Export NetCDF and Configure HEC-RAS Gridded Precipitation

Python

mrms_netcdf = PrecipMrms.to_ras_netcdf(
    mrms_stack,
    ras.project_folder / "P" / "mrms" / "mrms_qpe.nc",
    resolution=MRMS_RESOLUTION_METERS,
)
assert mrms_netcdf.exists() and mrms_netcdf.stat().st_size > 0

RasUnsteady.set_gridded_precipitation(
    UNSTEADY_NUMBER,
    mrms_netcdf,
    interpolation="Bilinear",
    ras_object=ras,
    dataset_name="APCP_surface",
)
updated_precip_config = RasUnsteady.get_met_precipitation_config(
    UNSTEADY_NUMBER,
    ras_object=ras,
)
assert updated_precip_config["enabled"] is True
assert updated_precip_config["source"] == "GDAL Raster File(s)"
assert updated_precip_config["gdal_group"] == "APCP_surface"

print(f"NetCDF: {mrms_netcdf}")
print("Updated precipitation configuration:")
for key in ["enabled", "mode", "source", "gdal_filename", "gdal_group"]:
    print(f"  {key}: {updated_precip_config.get(key)}")
print(f"Unsteady HDF sidecar: {Path(str(ras.project_folder / f'{ras.project_name}.u{UNSTEADY_NUMBER}') + '.hdf')}")

Compute the Plan

Python

RasPlan.update_simulation_date(
    PLAN_NUMBER,
    EVENT_START,
    SIM_END,
    ras_object=ras,
)
RasPlan.update_plan_intervals(
    PLAN_NUMBER,
    output_interval="1HOUR",
    instantaneous_interval="1HOUR",
    mapping_interval="1HOUR",
    ras_object=ras,
)

compute_result = RasCmdr.compute_plan(
    PLAN_NUMBER,
    ras_object=ras,
    force_rerun=True,
    verify=True,
    use_optimal_hdf_settings=True,
    hdf_output_variables=[
        "Cell Hydraulic Depth",
        "Cell Precipitation Rate",
        "Cell Cumulative Precipitation Depth",
    ],
)
assert bool(compute_result), "HEC-RAS compute failed"

PLAN_HDF = ras.project_folder / f"{ras.project_name}.p{PLAN_NUMBER}.hdf"
assert PLAN_HDF.exists() and PLAN_HDF.stat().st_size > 0
print(f"Compute result: {compute_result}")
print(f"Plan HDF: {PLAN_HDF}")

Validate Model Output Precipitation

Python

cumulative_summary = mesh_precip_summary(
    PLAN_HDF,
    "Cell Cumulative Precipitation Depth",
)
rate_summary = mesh_precip_summary(
    PLAN_HDF,
    "Cell Precipitation Rate",
)
assert int(cumulative_summary["positive_final_cells"].sum()) > 0
assert float(cumulative_summary["max"].max()) > 0.0
assert float(rate_summary["max"].max()) > 0.0

display(cumulative_summary)
display(rate_summary)

Python

mesh_data = HdfResultsMesh.get_mesh_cells_timeseries(
    PLAN_HDF,
    var="Cell Cumulative Precipitation Depth",
)
mesh_name, dataset = next(iter(mesh_data.items()))
cell_values = dataset["Cell Cumulative Precipitation Depth"]
time_index = pd.to_datetime(dataset["time"].values)
mean_cumulative = cell_values.mean(dim="cell_id").values
max_cumulative = cell_values.max(dim="cell_id").values

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(time_index, mean_cumulative, marker="o", label="Mesh mean")
ax.plot(time_index, max_cumulative, marker="o", label="Mesh max")
ax.set_title(f"HEC-RAS output precipitation applied to {mesh_name}")
ax.set_xlabel("Simulation time")
ax.set_ylabel("Cumulative precipitation depth (in)")
ax.grid(True, alpha=0.25)
ax.legend()
fig.autofmt_xdate()
plt.show()

print(
    f"Final mesh mean: {mean_cumulative[-1]:.3f} in; "
    f"final mesh max: {max_cumulative[-1]:.3f} in"
)

Validate Hydraulic Response

Python

depth_summary = mesh_precip_summary(
    PLAN_HDF,
    "Cell Hydraulic Depth",
).rename(
    columns={
        "mean_final": "mean_final_depth",
        "max": "max_depth",
        "positive_final_cells": "wet_final_cells",
    }
)
assert float(depth_summary["max_depth"].max()) > 0.0
assert int(depth_summary["wet_final_cells"].sum()) > 0

display(depth_summary)

depth_data = HdfResultsMesh.get_mesh_cells_timeseries(
    PLAN_HDF,
    var="Cell Hydraulic Depth",
)
mesh_name, depth_dataset = next(iter(depth_data.items()))
depth_values = depth_dataset["Cell Hydraulic Depth"]
depth_time = pd.to_datetime(depth_dataset["time"].values)
mean_depth = depth_values.mean(dim="cell_id").values
max_depth = depth_values.max(dim="cell_id").values

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(depth_time, mean_depth, marker="o", label="Mesh mean depth")
ax.plot(depth_time, max_depth, marker="o", label="Mesh max depth")
ax.set_title(f"HEC-RAS hydraulic depth response in {mesh_name}")
ax.set_xlabel("Simulation time")
ax.set_ylabel("Hydraulic depth (ft)")
ax.grid(True, alpha=0.25)
ax.legend()
fig.autofmt_xdate()
plt.show()

print(
    f"Final mesh mean depth: {mean_depth[-1]:.3f} ft; "
    f"peak mesh max depth: {max_depth.max():.3f} ft"
)