# Installed-package mode. Requires the Atlas14HyetographGrid API.
from pathlib import Path
import ras_commander
print(ras_commander.__file__)
Development mode: install this checkout in your environment from a terminal with
uv pip install -e ".[notebooks]". Restart the kernel afterward. Do not activate a
second import-path override. The shared-pattern API must be present in your installed
version; this notebook accompanies its introduction.
Atlas 14 Shared Temporal Pattern for an Existing Model¶
Apply one NOAA temporal curve scaled by each cell's storm-total depth to an existing 2D HEC-RAS model. This produces spatial variation in rainfall depth while holding the normalized temporal pattern constant. It differs from notebooks 726 and 727, which construct alternating blocks from multiple duration depths. Neither method is a universal substitute for the other; choose the method required by the study.
Workflow: select a plan and read its footprint → crop an original NOAA ASCII depth raster → generate interval-depth NetCDF → check totals, masks and time bounds → serialize intervals as Float32 → explicitly bilinear-project with HEC-Vortex to SHG DSS → inspect the resulting grids → optionally attach them to a new plan in a staged project copy → optionally qualify native preprocessing and its materialized forcing. This notebook never runs the final hydraulic solver.
The native NetCDF is the preferred retained design-storm product. DSS is a receiving model export. The historical-reproduction detail is Float32 per interval before spatial reprojection; projecting the total first can change the final floating-point values. Matching another delivery also requires its exact grid extent, CRS, interval bounds, NoData convention, and numerical runtime. This notebook does not promise bit-identical results for an arbitrary target model.
Inputs and dependencies¶
Use an existing model with a geometry HDF and a known CRS. Read the NOAA raster metadata to establish scale, frequency basis and CRS; do not infer them from a name. The example settings below describe Texas 100-year/24-hour PDS data, whose original ASCII values are thousandths of an inch. Select the temporal region from the NOAA regional map for your model; a state name alone does not select the region. For a domain crossing temporal regions, one shared curve is an explicit modeling assumption, not an automatic regional mosaic.
Install ras-commander[notebooks] and Rasterio in the same environment. DSS export
also needs HEC-Vortex/Jython; inspection uses ras-commander's built-in HEC-DSS Java
Monolith backend. See the precipitation guide for DSS setup and acknowledgments.
The built-in backend does not use pydsstools.
NOAA sources: depth raster GIS downloads, Texas temporal distributions. HEC documents frequency grids with a temporal pattern and alternating-block frequency storms.
from datetime import datetime, timedelta
from dataclasses import asdict, fields
import json
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from pyproj import CRS
from ras_commander import init_ras_project, RasPlan, RasUnsteady, RasMap, RasPreprocess
from ras_commander.RasProject import stage_project
from ras_commander.hdf import HdfMesh, HdfLandCover, HdfPlan
from ras_commander.precip import Atlas14HyetographGrid, VortexCli
from ras_commander.dss import RasDss
# Edit these paths and choices for your existing model. Inputs stay read-only.
SOURCE_PROJECT = Path("C:/Models/MyModel")
SOURCE_PLAN = "01"
RAS_VERSION = "6.6" # Receiving runtime; document any upgrade from the source version.
ASC_PATH = Path("C:/Rainfall/tx100yr24ha.asc")
SOURCE_CRS = "EPSG:4269"
SCALE_TO_INCHES = 0.001
DEPTH_BASIS = "PDS"
ARI_YEARS, DURATION_HOURS, INTERVAL_MINUTES = 100, 24, 30
STATE, REGION, QUARTILE, PROBABILITY = "tx", 3, "All Cases", "50%"
# 50% is a temporal-table column, NOT a peak-placement parameter or event AEP.
SOURCE_URL = "https://hdsc.nws.noaa.gov/pub/hdsc/data/tx/tx100yr24ha.zip"
repo_root = Path.cwd().parent if Path.cwd().name == "examples" else Path.cwd()
RUN_ROOT = repo_root / "working" / "example_731_shared_temporal"
RUN_ROOT.mkdir(parents=True, exist_ok=True)
RUN_DSS_EXPORT = False
ATTACH_TO_COPIED_MODEL = False
RUN_NATIVE_PRECOMPUTE = False # Explicit Windows precompute qualification; this never runs the final solver.
PRECOMPUTE_MAX_WAIT = 600
MATERIALIZED_HDF_NODATA_POLICY = None # Set to 'zero' or 'nan' only after observing the receiving runtime.
REPLACE_LEGACY_PRECIPITATION = False # Explicitly replace one reviewed area-rainfall boundary on the copy.
SANITIZE_LANDCOVER_NAMES = False # Optional native 6.6 labels-only compatibility migration on the copy.
PAD_DRY_TAIL = False # Explicitly opt in if the plan continues after the storm.
VORTEX_VERSION = "" # Record the installed release before enabling export.
VORTEX_PATH = None # Or Path("C:/Tools/vortex")
1. Select the model and its forcing footprint¶
Use the plan DataFrame to resolve the geometry. Crop in the ASCII source CRS; keep a spatial buffer for bilinear interpolation near the perimeter. The plan's simulation start anchors this synthetic design event, not an observed storm date.
Optional public eBFE example¶
The public East Galveston Bay model can be obtained through RasExamples:
from ras_commander import RasExamples
organized = RasExamples.organize_ebfe_model(
"east-galveston-bay", output_path=Path("C:/Models/eBFE"),
extract_ras_nested=True, include_results=False,
)
SOURCE_PROJECT = organized / "RAS Model" / "EastGalvBay.prj"
SOURCE_PLAN = "03"
RAS_VERSION = "6.6"
PAD_DRY_TAIL = True
REPLACE_LEGACY_PRECIPITATION = True
SANITIZE_LANDCOVER_NAMES = True
MATERIALIZED_HDF_NODATA_POLICY = "zero" # Observed on this qualified 6.6 receiving run.
The source download is about 10.7 GB and requires a Deflate64-capable 7-Zip
installation. Delivered result archives stay unexpanded by default. Plan 03 is
1PAC, a 48-hour HEC-RAS 5.0.7 model with an HMS precipitation-excess boundary.
Using 6.6 is an explicit receiving-runtime upgrade. The replacement below supplies
NOAA gross rainfall, not the original excess rainfall; review infiltration/losses
before using resulting flood levels for engineering decisions. This public eBFE
source is distinct from the RBFS East Galveston deliverable. Download and
organization do not qualify a model for computation.
The delivered land-cover names contain slashes rejected by the 6.6 geometry preprocessor. The explicit compatibility option uses HEC-RAS's native sanitizer on the copied sidecar, records its V1-to-V2 serialization migration, and verifies unchanged class codes, numerical parameters and companion raster bytes. It does not rebuild or resample the land-cover raster. This optional API is qualified for the 6.6 runtime only.
ras = init_ras_project(SOURCE_PROJECT, ras_version=RAS_VERSION, load_results_summary=False)
rows = ras.plan_df.loc[ras.plan_df.plan_number.astype(str).str.zfill(2).eq(SOURCE_PLAN)]
assert len(rows) == 1, "Select exactly one existing plan"
plan = rows.iloc[0]
assert pd.notna(plan["unsteady_number"]), "An unsteady-flow plan is required"
geometry_hdf = Path(str(plan["Geom Path"]) + ".hdf")
assert geometry_hdf.is_file(), "Existing geometry HDF required; no preprocessing is run"
footprint = HdfMesh.get_mesh_areas(geometry_hdf)
assert not footprint.empty and footprint.crs is not None
source_bounds = footprint.to_crs(SOURCE_CRS).total_bounds.tolist()
parts = [part.strip() for part in str(plan["Simulation Date"]).split(",")]
assert len(parts) == 4, "Expected start date/time and end date/time"
plan_times = []
for date_text, clock_text in [(parts[0], parts[1]), (parts[2], parts[3])]:
clock_text = clock_text.replace(":", "")
assert len(clock_text) == 4 and clock_text.isdigit(), "Expected HHMM plan clock"
hour, minute = int(clock_text[:2]), int(clock_text[2:])
assert (0 <= hour < 24 and 0 <= minute < 60) or (hour == 24 and minute == 0)
plan_times.append(datetime.strptime(date_text, "%d%b%Y") + timedelta(hours=hour, minutes=minute))
simulation_start, simulation_end = plan_times
storm_end = simulation_start + timedelta(hours=DURATION_HOURS)
assert simulation_end >= storm_end, "Plan must span the entire design storm"
print(footprint[["mesh_name"]])
print("Source-CRS bounds:", source_bounds)
print("Storm:", simulation_start, "through", storm_end)
2. Generate the shared-pattern design storm¶
For native cell depth $D(x,y)$ and cumulative NOAA fractions $F(t)$, each interval is $P_i(x,y)=D(x,y)[F(t_{i+1})-F(t_i)]$. No alternating-block rearrangement, areal reduction, or PDS/AMS conversion is applied. Intensity is interval depth divided by interval hours. A 100-year rainfall frequency does not establish a 1% flood AEP.
The ASCII route keeps source identity, CRS, scale and NoData explicit. The alternative
Atlas14HyetographGrid.generate(bounds=..., ...) uses the remote grid service;
its WGS84 bounds and source product must not be substituted silently when reproducing
an original ASCII workflow.
raw_nc = Atlas14HyetographGrid.generate_from_asc_file(
ASC_PATH, scale_factor=SCALE_TO_INCHES, source_crs=SOURCE_CRS,
ari_years=ARI_YEARS, storm_duration_hours=DURATION_HOURS,
timestep_minutes=INTERVAL_MINUTES, state=STATE, region=REGION,
quartile=QUARTILE, probability_column=PROBABILITY,
bounds=source_bounds, buffer_percent=10.0, depth_frequency_basis=DEPTH_BASIS,
output_netcdf=RUN_ROOT / "shared_pattern_native.nc",
)
with xr.open_dataset(raw_nc, engine="scipy", decode_times=False) as opened:
storm = opened.load()
print(storm)
3. Check conservation, missing cells and synchronization¶
Missing cells remain missing. Do not use a NaN-skipping sum that turns an entirely missing cell into a dry cell. Native-grid storm totals and normalized curves should agree before any reprojection. Reprojection can change pointwise totals at new cell locations; it must be evaluated on its own grid.
increments = storm.precip_incremental.values
native_total = storm.storm_total.values
fractions = storm.temporal_fraction.values
bounds_hours = storm.time_bounds.values
valid = np.isfinite(native_total)
assert valid.any()
np.testing.assert_array_equal(np.isfinite(increments), np.broadcast_to(valid, increments.shape))
assert (increments[:, valid] >= 0).all()
np.testing.assert_allclose(increments.sum(axis=0), native_total, rtol=1e-12, atol=1e-12, equal_nan=True)
np.testing.assert_allclose(increments[:, valid], fractions[:, None] * native_total[valid], rtol=1e-12, atol=1e-12)
np.testing.assert_allclose(fractions.sum(), 1.0, rtol=0, atol=1e-12)
np.testing.assert_allclose(np.diff(bounds_hours, axis=1), INTERVAL_MINUTES / 60)
np.testing.assert_allclose(bounds_hours[1:, 0], bounds_hours[:-1, 1])
np.testing.assert_allclose(bounds_hours[[0, -1], [0, 1]], [0, DURATION_HOURS])
np.testing.assert_allclose(storm.time.values, bounds_hours[:, 1])
assert CRS.from_wkt(storm.spatial_ref.attrs["crs_wkt"]).equals(CRS.from_user_input(SOURCE_CRS))
print("Native QC passed:", increments.shape[0], "intervals;", valid.sum(), "valid cells")
# Plot actual cell coordinates and interval boundaries, not array-index timestamps.
row, col = np.argwhere(valid)[len(np.argwhere(valid)) // 2]
fig, axes = plt.subplots(1, 3, figsize=(16, 4))
mesh = axes[0].pcolormesh(storm.x, storm.y, native_total, shading="auto", cmap="Blues")
footprint.to_crs(SOURCE_CRS).boundary.plot(ax=axes[0], color="black", linewidth=1)
fig.colorbar(mesh, ax=axes[0], label="Storm total (in)")
axes[0].set_title("NOAA cell depths and model footprint")
axes[0].set_xlabel("Source CRS x"); axes[0].set_ylabel("Source CRS y")
dt = INTERVAL_MINUTES / 60
axes[1].bar(bounds_hours[:, 0], increments[:, row, col] / dt, width=dt, align="edge")
axes[1].set(xlabel="Hours since storm start", ylabel="Intensity (in/hr)", title="Sample-cell hyetograph")
axes[2].plot(np.r_[0, bounds_hours[:, 1]], np.r_[0, np.cumsum(increments[:, row, col])])
axes[2].set(xlabel="Hours since storm start", ylabel="Cumulative depth (in)", title="Same sample cell and clock")
plt.tight_layout()
4. Prepare Float32 interval NetCDF with the model clock¶
Serialize each native-grid interval depth, then let Vortex reproject each frame. Do not reproject the storm-total grid first. Keep the numeric time bounds and their units together. No leading missing interval is invented: normal design rainfall starts at the model start. Reproducing a delivery with a leading missing record or offset clock is a separate compatibility choice that must be documented explicitly.
If the source plan outlasts the storm, PAD_DRY_TAIL=True appends explicit dry
intervals at valid cells while retaining spatial NoData. This assumes no additional
rain after the design event. Otherwise export stops rather than shortening the plan.
Vortex classification: set the interval variable's long_name="precipitation"
and standard_name="lwe_thickness_of_precipitation_amount" explicitly. Vortex
0.13.3 classifies the variable before applying the requested DSS C-part; setting
C=PRECIPITATION alone can still produce a PER-AVER fallback for an unrecognized
variable. Retain CF cell_methods="time: sum" and verify IN/PER-CUM after export.
See Vortex DSS writer.
forcing = storm[["precip_incremental", "time_bounds", "spatial_ref"]].copy(deep=True)
end_hours = (simulation_end - simulation_start).total_seconds() / 3600
if end_hours > DURATION_HOURS:
if not PAD_DRY_TAIL:
raise ValueError("Plan continues after storm: explicitly enable PAD_DRY_TAIL or review the plan window")
count = int(round(end_hours / dt))
assert np.isclose(count * dt, end_hours, rtol=0, atol=1e-9), "Plan end must align with the rainfall interval"
all_bounds = np.column_stack((np.arange(count) * dt, (np.arange(count) + 1) * dt))
padded = np.full((count, *native_total.shape), np.nan, dtype=np.float64)
padded[:, valid] = 0.0
padded[:len(increments)] = increments
forcing = forcing.drop_dims("time").assign_coords(time=all_bounds[:, 1])
forcing["precip_incremental"] = (("time", "y", "x"), padded, dict(storm.precip_incremental.attrs))
forcing["time_bounds"] = (("time", "bounds"), all_bounds)
forcing["precip_incremental"] = forcing.precip_incremental.astype("float32")
forcing.precip_incremental.attrs.update(units="in", cell_methods="time: sum", long_name="precipitation", standard_name="lwe_thickness_of_precipitation_amount", grid_mapping="spatial_ref")
clock_units = "hours since " + simulation_start.strftime("%Y-%m-%d %H:%M:%S")
forcing.time.attrs.update(units=clock_units, calendar="standard", bounds="time_bounds", axis="T", standard_name="time")
forcing.time_bounds.attrs.update(units=clock_units, calendar="standard")
forcing.attrs.update(source_url=SOURCE_URL, spatial_resampling="Bilinear", target_grid="SHG, 2000 m", precision_order="native interval -> float32 -> spatial reprojection")
float32_nc = RUN_ROOT / "shared_pattern_float32_model_clock.nc"
assert not float32_nc.exists(), "Use a new run folder rather than overwriting evidence"
forcing.to_netcdf(float32_nc, engine="scipy")
with xr.open_dataset(float32_nc, engine="scipy", decode_times=False) as check:
np.testing.assert_array_equal(check.precip_incremental.values, forcing.precip_incremental.values)
np.testing.assert_array_equal(check.time_bounds.values, forcing.time_bounds.values)
assert check.precip_incremental.dtype == np.dtype("float32")
provenance = dict(storm.attrs)
provenance.update(source_url=SOURCE_URL, plan=SOURCE_PLAN, source_project=str(SOURCE_PROJECT),
simulation_start=str(simulation_start), simulation_end=str(simulation_end),
resampling="Bilinear", target_cell_size_m=2000, float32_before_warp=True,
areal_reduction="none", dry_tail=PAD_DRY_TAIL)
(RUN_ROOT / "provenance.json").write_text(json.dumps(provenance, indent=2, default=str), encoding="utf-8")
5. Optional DSS export and read-back¶
Set RUN_DSS_EXPORT=True only after the native QC is satisfactory. Bilinear is
explicit: underlying programmatic Vortex callers can otherwise select a different
fallback. Vortex determines the output extent from the input domain here; an exact
delivery comparison additionally needs the delivered grid origin and extent.
Record the installed Vortex version with the run evidence.
Read every resulting record with RasDss.read_grid. Check chronology, units, mask,
geometry and totals before attaching it. These are rainfall-data checks, not a
hydraulic validation. The export uses interval amounts (IN, PER-CUM), not
intensity or cumulative-since-storm-start values.
dss_path = RUN_ROOT / "shared_pattern_shg.dss"
if RUN_DSS_EXPORT:
assert VORTEX_VERSION.strip(), "Record the installed Vortex release in VORTEX_VERSION"
provenance["vortex_version"] = VORTEX_VERSION
provenance["vortex_path"] = str(VortexCli.find_vortex(VORTEX_PATH))
(RUN_ROOT / "provenance.json").write_text(json.dumps(provenance, indent=2, default=str), encoding="utf-8")
assert not dss_path.exists(), "Use a new DSS file"
VortexCli.import_gridded(
input_files=[float32_nc], output_dss=dss_path, variables=["precip_incremental"],
target_wkt="SHG", target_cell_size=2000, resampling_method="Bilinear",
dss_parts={"A": "SHG", "B": "DESIGN", "C": "PRECIPITATION", "F": "ATLAS14_SHARED"},
vortex_path=VORTEX_PATH,
)
catalog = RasDss.get_catalog(dss_path)
frames = [RasDss.read_grid(dss_path, p) for p in catalog.pathname]
assert frames, "Vortex produced no grid records"
ordered = sorted(zip(catalog.pathname, frames), key=lambda pair: pair[1]["start_time"])
dss_paths, frames = zip(*ordered)
assert len(frames) == len(forcing.time)
for i, frame in enumerate(frames):
expected = [simulation_start + timedelta(hours=float(h)) for h in forcing.time_bounds.values[i]]
assert pd.Timestamp(frame["start_time"]) == expected[0]
assert pd.Timestamp(frame["end_time"]) == expected[1]
origin_keys = ("lower_left_cell", "origin", "projection")
for frame in frames:
assert all(frame["metadata"][key] == frames[0]["metadata"][key] for key in origin_keys), "Grid origin changed between intervals"
assert frame["units"].upper() == "IN" and frame["data_type"].upper() == "PER-CUM"
assert frame["cell_size"] == 2000 and frame["shape"] == frames[0]["shape"]
assert CRS.from_user_input(frame["crs"]).equals(CRS.from_user_input(frames[0]["crs"]))
dss_values = np.stack([frame["data"] for frame in frames])
assert np.isfinite(dss_values).any() and np.nanmin(dss_values) >= 0
np.testing.assert_array_equal(np.isfinite(dss_values), np.broadcast_to(np.isfinite(dss_values[0]), dss_values.shape))
# Public grid arrays are south-first; raster masks are north-first.
from rasterio.features import geometry_mask
from rasterio.transform import from_origin
meta = frames[0]["metadata"]
cell = float(frames[0]["cell_size"])
left, bottom = meta["origin"]
height, width = frames[0]["shape"]
right, top = left + width * cell, bottom + height * cell
projected_footprint = footprint.to_crs(frames[0]["crs"])
west, south, east, north = projected_footprint.total_bounds
assert left <= west and bottom <= south and right >= east and top >= north, "DSS must cover the entire model footprint"
model_cells = geometry_mask(projected_footprint.geometry, out_shape=(height, width),
transform=from_origin(left, top, cell, cell), invert=True, all_touched=True)[::-1]
assert model_cells.any() and np.isfinite(dss_values[:, model_cells]).all(), "No missing rainfall over model cells"
# Float32 serialization/reprojection allows small rounding differences.
# Check the shared curve independently of the new spatial cell depths.
wet = dss_values[:len(increments), model_cells].astype(np.float64)
projected_total = wet.sum(axis=0)
np.testing.assert_allclose(wet, fractions[:, None] * projected_total[None, :],
rtol=5e-6, atol=1e-7, err_msg="Export changed the shared temporal curve")
if PAD_DRY_TAIL and len(frames) > len(increments):
assert np.nanmax(np.abs(dss_values[len(increments):])) == 0
print("DSS intervals and masks passed; inspect metadata and spatial coverage:")
print({k: v for k, v in frames[0].items() if k != "data"})
print("Projected cell totals (in):", np.nanmin(dss_values.sum(axis=0)), np.nanmax(dss_values.sum(axis=0)))
else:
print("Retained NetCDF ready; DSS export not requested")
6. Optional attachment to a copied model¶
This step stages a separate project, clones both unsteady flow and plan, and uses
configure_gridded_dss_precipitation. It does not edit the original project or
change its simulation window. The staging API checks referenced assets; resolve
missing or external dependencies before treating the copy as portable. Review
inherited boundary flows, infiltration and other meteorology for consistency with
the design event. Conflicting precipitation boundaries must be reviewed explicitly. With
REPLACE_LEGACY_PRECIPITATION=True, this example removes exactly one superseded
area-rainfall block from the staged copy before cloning the new flow and plan.
The default refuses removal. Multiple legacy rainfall blocks require a separately
reviewed preparation step; they are not silently deleted. The source tree is
unchanged, while its original flow copy in the staged tree is also superseded.
After attachment, inspect rainfall coverage and totals in the receiving model. The optional precompute qualification below checks the receiving temporary HDF; follow the broader qualification checks in notebooks 727/728 before computing. A successful DSS write is not proof that the hydraulic model consumes it correctly. In particular, inspect interval alignment in the preprocessed temporary HDF when using HEC-RAS 6.6; period-average grid ingestion can introduce a one-interval timing shift. Do not compensate by silently shifting the source timestamps or accept a shifted rainfall series as a passing check.
if ATTACH_TO_COPIED_MODEL:
assert RUN_DSS_EXPORT and dss_path.is_file(), "Export and inspect the DSS first"
staged = stage_project(SOURCE_PROJECT, RUN_ROOT / "model_copy")
removal = None
# Inspect before any plan/flow edits: deletion uses the staging evidence.
boundaries = RasUnsteady.inspect_boundary_blocks(staged, unsteady_number=plan["unsteady_number"])
legacy = boundaries.loc[boundaries["bc_type"].eq("Precipitation Hydrograph")]
if not legacy.empty:
print(legacy[["boundary_id", "area_2d", "bc_type", "boundary_location_raw"]].to_string(index=False))
if not REPLACE_LEGACY_PRECIPITATION or len(legacy) != 1:
raise ValueError("Review legacy rainfall and explicitly opt in to replacing one block; the source is unchanged.")
boundary = legacy.iloc[0]
removal = RasUnsteady.delete_boundary(
staged, unsteady_number=plan["unsteady_number"],
boundary_id=boundary["boundary_id"],
expected_source_sha256=boundary["owner_sha256"],
expected_block_sha256=boundary["block_sha256"],
expected_bc_type=boundary["bc_type"],
expected_location_raw=boundary["boundary_location_raw"], dry_run=False,
)
assert removal.state == "applied"
print("Superseded legacy rainfall on copy:", removal.boundary_location_raw)
# Staging initially resolves the delivered version. Use the chosen receiver.
copied = init_ras_project(staged.ras_object.project_folder, RAS_VERSION,
ras_object=True, load_results_summary=False)
landcover_migrations = []
if SANITIZE_LANDCOVER_NAMES:
layers = RasMap.list_landcover_layers(copied.project_folder)
if layers.empty:
raise ValueError("No registered land-cover layer to sanitize; review the compatibility option.")
for layer_path in layers["resolved_path"]:
registered = Path(layer_path)
if registered.suffix.lower() == ".hdf":
sidecar, exact_raster = registered, None # Native HDF association resolves the TIFF.
elif registered.suffix.lower() in {".tif", ".tiff"}:
sidecar, exact_raster = registered.with_suffix(".hdf"), registered
else:
raise ValueError(f"Unsupported registered land-cover path: {registered}")
if not sidecar.is_file():
raise FileNotFoundError(f"Review the native land-cover sidecar: {sidecar}")
migrated = HdfLandCover.sanitize_classification_names(
sidecar, raster_path=exact_raster, hecras_version=RAS_VERSION,
)
landcover_migrations.append(asdict(migrated))
new_unsteady = RasPlan.clone_unsteady(plan["unsteady_number"], new_title="Atlas14 shared pattern", ras_object=copied)
new_plan = RasPlan.clone_plan(SOURCE_PLAN, new_plan_shortid="Atlas14 shared", new_title="Atlas14 shared pattern",
unsteady_flow=new_unsteady, ras_object=copied)
RasUnsteady.configure_gridded_dss_precipitation(
new_unsteady, dss_filename=str(dss_path.resolve()), dss_pathname=dss_paths[0],
ratio=1.0, interpolation="Bilinear", ras_object=copied,
)
config = RasUnsteady.get_met_precipitation_config(new_unsteady, ras_object=copied)
assert config["mode"] == "Gridded" and config["source"] == "DSS"
configured_dss = Path(config["dss_filename"])
if not configured_dss.is_absolute():
configured_dss = copied.project_folder / configured_dss
assert configured_dss.resolve() == dss_path.resolve()
assert config["dss_pathname"] == dss_paths[0]
attachment_receipt = {
"source_project": str(SOURCE_PROJECT), "source_plan": SOURCE_PLAN,
"receiving_version": RAS_VERSION, "copied_plan": new_plan, "copied_flow": new_unsteady,
"staging": {field.name: getattr(staged, field.name) for field in fields(staged)
if field.name not in {"assets", "ras_object"}},
"superseded_boundary": asdict(removal) if removal is not None else None,
"precipitation_config": config,
"landcover_name_migrations": landcover_migrations,
}
(RUN_ROOT / "attachment-receipt.json").write_text(
json.dumps(attachment_receipt, indent=2, default=str), encoding="utf-8",
)
print("Prepared copied plan:", new_plan, "in", copied.project_folder)
print("Retain the sibling DSS with this workspace; its reference is absolute.")
else:
print("Original model unchanged. Attachment is an explicit optional step.")
7. Optional native precompute qualification¶
Set RUN_NATIVE_PRECOMPUTE=True only after attachment and DSS read-back pass. The
public RasPreprocess.preprocess_plan() API starts Windows HEC-RAS, waits for its
preprocess handoff, stops the owned process, and validates the completed temporary
HDF. This section opens that HDF only after the API returns; never inspect an
actively written HDF from another process.
The read-back binds the materialized interval depths and interval-end timestamps to the already-reviewed DSS frames. It also checks the receiving configuration, grid attributes, spatial row order, and an explicitly observed NoData policy. It creates a precompute receipt for a qualified solver input, not a solver result or an engineering acceptance decision.
if RUN_NATIVE_PRECOMPUTE:
assert ATTACH_TO_COPIED_MODEL and RUN_DSS_EXPORT, (
"Native precompute requires the staged copy and reviewed DSS frames."
)
assert MATERIALIZED_HDF_NODATA_POLICY in {"zero", "nan"}, (
"Set MATERIALIZED_HDF_NODATA_POLICY from an observed receiving-HDF result."
)
prepared = RasPreprocess.preprocess_plan(
new_plan, ras_object=copied, max_wait=PRECOMPUTE_MAX_WAIT,
)
assert prepared, prepared.error
# RasPreprocess has completed its owned Windows lifecycle before this HDF read.
materialized = HdfPlan.get_plan_met_precip_values(prepared.tmp_hdf_path)
hdf_values = np.asarray(materialized["values"], dtype=np.float64)
raw_hdf_times = np.asarray(materialized["timestamps"])
hdf_times = pd.DatetimeIndex([
pd.Timestamp(datetime.strptime(
bytes(value).decode("ascii").strip() if isinstance(value, (bytes, np.bytes_)) else str(value).strip(),
"%d%b%Y %H:%M:%S.%f",
))
for value in raw_hdf_times
])
expected_ends = pd.DatetimeIndex([pd.Timestamp(frame["end_time"]) for frame in frames])
assert hdf_values.shape[0] == len(hdf_times) == len(frames)
assert hdf_values.shape[1] == int(np.prod(frames[0]["shape"]))
assert hdf_times.equals(expected_ends), "Materialized interval-end timestamps differ from DSS"
received = materialized["attributes"]
assert received["Mode"] == "Gridded" and received["Source"] == "DSS"
assert received["Data Type"].casefold() == "per-cum"
assert received["Units"].casefold() in {"in", "inch", "inches"}
assert float(received["Ratio"]) == 1.0
assert received["Interpolation Method"].casefold() == "bilinear"
assert received["DSS Pathname"] == config["dss_pathname"] == dss_paths[0]
received_dss = Path(received["DSS Filename"])
if not received_dss.is_absolute():
received_dss = prepared.tmp_hdf_path.parent / received_dss
assert received_dss.resolve() == dss_path.resolve()
rows, cols = (int(value) for value in frames[0]["shape"])
left, bottom = (float(value) for value in frames[0]["metadata"]["origin"])
cell = float(frames[0]["cell_size"])
assert int(received["Raster Rows"]) == rows and int(received["Raster Cols"]) == cols
assert float(received["Raster Cellsize"]) == cell
assert float(received["Raster Left"]) == left and float(received["Raster Top"]) == bottom + rows * cell
assert CRS.from_user_input(received["Projection"]).equals(CRS.from_user_input(frames[0]["crs"]))
orientations = {
"south_first_c": lambda grid: grid.reshape(-1, order="C"),
"north_first_c": lambda grid: np.flipud(grid).reshape(-1, order="C"),
"south_first_f": lambda grid: grid.reshape(-1, order="F"),
"north_first_f": lambda grid: np.flipud(grid).reshape(-1, order="F"),
}
matches = []
for name, transform in orientations.items():
expected = np.stack([transform(grid) for grid in dss_values]).astype(np.float64)
finite = np.isfinite(expected)
nodata_ok = (
np.allclose(hdf_values[~finite], 0.0, rtol=0, atol=5e-6)
if MATERIALIZED_HDF_NODATA_POLICY == "zero"
else np.isnan(hdf_values[~finite]).all()
)
values_ok = np.isfinite(hdf_values[finite]).all() and np.allclose(
hdf_values[finite], expected[finite], rtol=0, atol=5e-6,
)
if nodata_ok and values_ok:
matches.append((name, expected, finite))
assert len(matches) == 1, "Materialized grid orientation is absent or ambiguous"
orientation, expected_values, finite_values = matches[0]
finite_cells = finite_values.all(axis=0)
assert finite_cells.any(), "No finite source-grid cells remain for the cumulative comparison"
interval_error = np.abs(hdf_values[:, finite_cells] - expected_values[:, finite_cells])
cumulative_hdf = np.cumsum(hdf_values[:, finite_cells], axis=0)
cumulative_expected = np.cumsum(expected_values[:, finite_cells], axis=0)
cumulative_error = np.abs(cumulative_hdf - cumulative_expected)
np.testing.assert_allclose(cumulative_hdf, cumulative_expected, rtol=0, atol=5e-6)
tail = hdf_times > storm_end
if PAD_DRY_TAIL:
assert tail.any() and np.allclose(hdf_values[tail][:, finite_cells], 0.0, rtol=0, atol=5e-6)
else:
assert not tail.any()
precompute_receipt = {
"prepared": asdict(prepared),
"tmp_hdf": str(prepared.tmp_hdf_path),
"values_shape": list(hdf_values.shape),
"timestamp_count": int(len(hdf_times)),
"dss_pathname": dss_paths[0],
"orientation": orientation,
"nodata_policy": MATERIALIZED_HDF_NODATA_POLICY,
"interval_max_abs_error": float(np.max(interval_error)),
"cumulative_max_abs_error": float(np.max(cumulative_error)),
"dry_tail_rows": int(tail.sum()),
}
(RUN_ROOT / "precompute-receipt.json").write_text(
json.dumps(precompute_receipt, indent=2, default=str), encoding="utf-8",
)
print("Native precompute and materialized forcing passed:", prepared.tmp_hdf_path)
else:
print("Native precompute not run. Enable only after attachment and DSS read-back review.")