"""UnsteadyUnsteadyPlan - read HEC-RAS unsteady plan HDF5 files (.p*.hdf).
Plan HDF files embed the same ``Geometry/`` group as geometry HDF files
*plus* ``Results/Unsteady/...`` time-series and summary output.
``UnsteadyUnsteadyPlan`` inherits ``Geometry`` so all geometry accessors are available.
``FlowAreaResults`` extends ``FlowArea`` with lazy time-series properties,
summary DataFrames, and computed depth / velocity methods. Raster export
methods delegate to ``rivia.geo`` via a deferred import so this module is
fully usable without rasterio or scipy installed.
"""
from __future__ import annotations
import dataclasses
import datetime as dt
import logging
import math
from abc import abstractmethod
from collections.abc import Callable, Iterable, Iterator
from pathlib import Path
from typing import TYPE_CHECKING, Any, Literal, overload
import numpy as np
import pandas as pd
from rivia.utils import log_call, parse_hec_datetime, parse_interval, timed
from ._base import _PlanHdf, _RAS_TS_FMT, _parse_hec_ts_array
from .geometry import (
_SA_ROOT,
Bridge,
CrossSection,
CrossSectionCollection,
FlowArea,
FlowAreaCollection,
Geometry,
InlineStructure,
LateralStructure,
SA2DConnection,
StorageArea,
StorageAreaCollection,
Structure,
StructureCollection,
_decode,
)
from .log import UnsteadyRuntimeLog
if TYPE_CHECKING:
import h5py
import rasterio.io
logger = logging.getLogger("rivia.hdf")
# ---------------------------------------------------------------------------
# HDF path constants
# ---------------------------------------------------------------------------
_TS_ROOT = "Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series"
_SUM_ROOT = "Results/Unsteady/Output/Output Blocks/Base Output/Summary Output"
_TS_2D = f"{_TS_ROOT}/2D Flow Areas"
_SUM_2D = f"{_SUM_ROOT}/2D Flow Areas"
_TS_SA = f"{_TS_ROOT}/Storage Areas"
_SUM_SA = f"{_SUM_ROOT}/Storage Areas"
_TIME_DS = f"{_TS_ROOT}/Time"
_TIME_STAMP_DS = f"{_TS_ROOT}/Time Date Stamp"
# Base Output structure group paths (for output="mapping")
_TS_INLINE = f"{_TS_ROOT}/Inline Structures"
_TS_LATERAL = f"{_TS_ROOT}/Lateral Structures"
_TS_BRIDGE = f"{_TS_ROOT}/Bridges"
_TS_SA_CONN = f"{_TS_ROOT}/SA 2D Area Conn"
_DSS_ROOT = (
"Results/Unsteady/Output/Output Blocks"
"/DSS Hydrograph Output/Unsteady Time Series"
)
_DSS_SA_CONN = f"{_DSS_ROOT}/SA 2D Area Conn"
_DSS_INLINE = f"{_DSS_ROOT}/Inline Structures"
_DSS_LATERAL = f"{_DSS_ROOT}/Lateral Structures"
_DSS_BRIDGE = f"{_DSS_ROOT}/Bridges"
_DSS_TIME_STAMP_DS = f"{_DSS_ROOT}/Time Date Stamp"
_TS_XS = f"{_TS_ROOT}/Cross Sections"
_DSS_XS = f"{_DSS_ROOT}/Cross Sections"
# DSS Profile Output (for output="profile")
_DSS_PROF_ROOT = (
"Results/Unsteady/Output/Output Blocks"
"/DSS Profile Output/Unsteady Time Series"
)
_DSS_PROF_XS = f"{_DSS_PROF_ROOT}/Cross Sections"
_DSS_PROF_INLINE = f"{_DSS_PROF_ROOT}/Inline Structures"
_DSS_PROF_LATERAL = f"{_DSS_PROF_ROOT}/Lateral Structures"
_DSS_PROF_BRIDGE = f"{_DSS_PROF_ROOT}/Bridges"
_DSS_PROF_SA_CONN = f"{_DSS_PROF_ROOT}/SA 2D Area Conn"
_DSS_PROF_SA = f"{_DSS_PROF_ROOT}/Storage Areas"
_DSS_PROF_TIME_STAMP_DS = f"{_DSS_PROF_ROOT}/Time Date Stamp"
_RUN_SUM = "Results/Unsteady/Summary"
_VOL_ACC = f"{_RUN_SUM}/Volume Accounting"
_VOL_1D = f"{_VOL_ACC}/Volume Accounting 1D"
_VOL_2D = f"{_VOL_ACC}/Volume Accounting 2D"
# Post Process Profiles paths (for output="post_process")
_POSTPROC_ROOT = (
"Results/Post Process/Steady/Output/Output Blocks"
"/Base Output/Post Process/Post Process Profiles"
)
_POSTPROC_PROFILE_DATES = f"{_POSTPROC_ROOT}/Profile Dates"
_POSTPROC_XS = f"{_POSTPROC_ROOT}/Cross Sections"
_POSTPROC_GEOM_ATTRS = (
"Results/Post Process/Steady/Output/Geometry Info"
"/Cross Section Attributes"
)
_POSTPROC_LATERAL = f"{_POSTPROC_ROOT}/Lateral Structures"
_POSTPROC_SA_CONN = f"{_POSTPROC_ROOT}/SA 2D Area Conn"
_POSTPROC_SA = f"{_POSTPROC_ROOT}/Storage Areas"
# ---------------------------------------------------------------------------
# _FlowAreaResultsDerived - abstract parent with all derived/computed methods
# ---------------------------------------------------------------------------
class _FlowAreaResultsDerived(FlowArea):
"""Abstract parent implementing all derived/computed methods for a 2-D flow area.
Declares the minimal set of abstract properties that
:class:`FlowAreaResults` provides as direct HDF reads, allowing this class
to be type-checked without depending on ``_ts`` / ``_sum`` group handles.
Concrete subclasses must implement the six abstract properties:
``water_surface``, ``face_velocity``, ``face_flow``, ``timestamps``,
``max_water_surface``, ``_max_water_surface``.
New methods that compute a result from HDF data (any arithmetic, pipeline
calls, or pandas wrapping beyond a bare slice) belong in this class.
New methods that directly return an HDF dataset or a thin slice/summary
of one belong in :class:`FlowAreaResults`.
"""
# ------------------------------------------------------------------
# Abstract stubs -- implemented by FlowAreaResults
# ------------------------------------------------------------------
@property
@abstractmethod
def water_surface(self) -> "h5py.Dataset":
"""Water-surface elevation time series, shape ``(n_t, n_cells + n_ghost)``."""
...
@property
@abstractmethod
def face_velocity(self) -> "h5py.Dataset":
"""Signed face-normal velocity time series, shape ``(n_t, n_faces)``."""
...
@property
@abstractmethod
def face_flow(self) -> "h5py.Dataset | None":
"""Volumetric face-flow time series, or ``None`` if not output."""
...
@property
@abstractmethod
def timestamps(self) -> pd.DatetimeIndex:
"""Result output timestamps as a ``pd.DatetimeIndex``."""
...
@property
@abstractmethod
def max_water_surface(self) -> pd.DataFrame:
"""Maximum WSE per cell, columns ``['value', 'time']``. Real cells only."""
...
@property
@abstractmethod
def _max_water_surface(self) -> pd.DataFrame:
"""Maximum WSE including ghost cell rows."""
...
# ------------------------------------------------------------------
# Unified get_* accessors
# ------------------------------------------------------------------
@overload
def get_water_surface(
self,
*,
timestep: int,
cell: None = ...,
include_ghost: bool = ...,
) -> np.ndarray: ...
@overload
def get_water_surface(
self,
*,
timestep: None = ...,
cell: int,
include_ghost: bool = ...,
) -> pd.Series: ...
@overload
def get_water_surface(
self,
*,
timestep: int,
cell: int,
include_ghost: bool = ...,
) -> float: ...
def get_water_surface(
self,
*,
timestep: int | None = None,
cell: int | None = None,
include_ghost: bool = False,
) -> np.ndarray | pd.Series | float:
"""Water-surface elevation snapshot, time series, or scalar.
Parameters
----------
timestep : int, optional
0-based index into the time dimension.
cell : int, optional
0-based real-cell index.
include_ghost : bool, optional
When ``True`` and *cell* is ``None``, include ghost-cell columns
in the returned snapshot. Raises ``ValueError`` when combined
with *cell*.
Returns
-------
ndarray, shape ``(n_cells,)``
Snapshot for all real cells when ``timestep=t, cell=None,
include_ghost=False``.
ndarray, shape ``(n_cells + n_ghost,)``
Snapshot including ghost cells when ``timestep=t, cell=None,
include_ghost=True``.
pd.Series
WSE over time indexed by :attr:`timestamps`, named
``"Water Surface"``, when ``timestep=None, cell=c``.
float
Single scalar when both *timestep* and *cell* are given.
Raises
------
ValueError
If both *timestep* and *cell* are ``None``.
If *include_ghost* is ``True`` and *cell* is not ``None``.
"""
if timestep is None and cell is None:
raise ValueError("At least one of timestep or cell must be specified.")
if include_ghost and cell is not None:
raise ValueError("include_ghost=True is incompatible with cell selection.")
if timestep is not None and cell is None:
if include_ghost:
return np.array(self.water_surface[timestep, :])
return np.array(self.water_surface[timestep, : self.n_cells])
if timestep is None and cell is not None:
values = np.array(self.water_surface[:, cell])
return pd.Series(values, index=self.timestamps, name="Water Surface")
# Both specified -> scalar
return float(self.water_surface[timestep, cell]) # type: ignore[index]
@overload
def get_depth(self, *, timestep: int, cell: None = ...) -> np.ndarray: ...
@overload
def get_depth(self, *, timestep: None = ..., cell: int) -> pd.Series: ...
@overload
def get_depth(self, *, timestep: int, cell: int) -> float: ...
def get_depth(
self,
*,
timestep: int | None = None,
cell: int | None = None,
) -> np.ndarray | pd.Series | float:
"""Water depth snapshot, time series, or scalar.
Depth is ``max(0, WSE - cell_min_elevation)``.
Parameters
----------
timestep : int, optional
0-based index into the time dimension.
cell : int, optional
0-based real-cell index.
Returns
-------
ndarray, shape ``(n_cells,)``
Depth at every real cell when ``timestep=t, cell=None``.
pd.Series
Depth over time indexed by :attr:`timestamps`, named ``"Depth"``,
when ``timestep=None, cell=c``.
float
Single scalar when both are given.
Raises
------
ValueError
If both *timestep* and *cell* are ``None``.
"""
if timestep is None and cell is None:
raise ValueError("At least one of timestep or cell must be specified.")
if timestep is not None and cell is None:
wse = np.array(self.water_surface[timestep, : self.n_cells])
return np.maximum(0.0, wse - self.cell_min_elevation)
if timestep is None and cell is not None:
wse = np.array(self.water_surface[:, cell])
depth = np.maximum(0.0, wse - self.cell_min_elevation[cell])
return pd.Series(depth, index=self.timestamps, name="Depth")
ws = float(self.water_surface[timestep, cell]) # type: ignore[index]
bed = float(self.cell_min_elevation[cell])
return max(0.0, ws - bed)
def get_max_depth(self) -> pd.DataFrame:
"""Maximum depth per cell using the time of maximum WSE.
``value = max(0, max_WSE - bed_elevation)``.
``time`` is the elapsed time (days) of maximum WSE; maximum depth
may not coincide exactly with maximum WSE.
Returns
-------
pd.DataFrame
Columns ``['value', 'time']``, index = cell index.
"""
df = self.max_water_surface.copy()
df["value"] = np.maximum(0.0, df["value"].to_numpy() - self.cell_min_elevation)
return df
def get_wet_cells(self, timestep: int, depth_min: float = 0.0) -> np.ndarray:
"""Boolean mask of wet cells for one timestep.
A cell is wet when ``WSE - cell_min_elevation > depth_min``.
Parameters
----------
timestep : int
0-based index into the time dimension.
depth_min : float, optional
Minimum depth threshold in model units. Default ``0.0``.
Returns
-------
ndarray, shape ``(n_cells,)``, dtype bool
"""
wse = np.array(self.water_surface[timestep, : self.n_cells])
return (wse - self.cell_min_elevation) > depth_min
def get_wet_faces(self, timestep: int, depth_min: float = 0.0) -> np.ndarray:
"""Boolean mask of wet faces for one timestep.
A face is wet when at least one of its adjacent cells is wet
(see :meth:`get_wet_cells`). Boundary faces are wet when their
single adjacent real cell is wet.
Parameters
----------
timestep : int
0-based index into the time dimension.
depth_min : float, optional
Minimum cell depth threshold passed to :meth:`get_wet_cells`.
Returns
-------
ndarray, shape ``(n_faces,)``, dtype bool
"""
wc = self.get_wet_cells(timestep, depth_min)
fci = self.face_cell_indexes # (n_faces, 2)
n = self.n_cells
c0_wet = np.where(fci[:, 0] >= 0, wc[np.clip(fci[:, 0], 0, n - 1)], False)
c1_wet = np.where(fci[:, 1] >= 0, wc[np.clip(fci[:, 1], 0, n - 1)], False)
return c0_wet | c1_wet
def get_cell_velocity(
self,
timestep: int,
*,
component: Literal["vector", "speed", "angle"] = "vector",
method: Literal[
"area_weighted", "length_weighted", "flow_ratio"
] = "area_weighted",
wse_interp: Literal["average", "sloped", "max"] = "average",
face_velocity_location: Literal[
"centroid", "normal_intercept"
] = "normal_intercept",
) -> np.ndarray:
"""Reconstruct cell-centre velocity for one timestep.
Uses the WLS method prescribed by the HEC-RAS Technical Reference
Manual (Section: 2D Unsteady Flow - Cell Velocity).
Parameters
----------
timestep : int
0-based index into the time dimension.
component : {"vector", "speed", "angle"}, optional
``"vector"`` (default): returns ``[Vx, Vy]`` pairs, shape
``(n_cells, 2)``.
``"speed"``: returns velocity magnitude ``sqrt(Vx^2 + Vy^2)``,
shape ``(n_cells,)``.
``"angle"``: returns flow direction in degrees clockwise from
north; ``nan`` where speed < 1e-10, shape ``(n_cells,)``.
method : {"area_weighted", "length_weighted", "flow_ratio"}, optional
``"area_weighted"`` (default, matches HEC-RAS): weights are
wetted face flow areas from the hydraulic property tables.
``"length_weighted"``: weights are face plan-view lengths.
``"flow_ratio"``: requires ``Face Flow`` output; back-calculates
flow area as ``|Q|/|V_n|``.
wse_interp : {"average", "sloped", "max"}, optional
How to estimate face WSE when ``method="area_weighted"``.
``"average"`` (default): simple mean of the two adjacent cell WSEs.
``"sloped"``: distance-weighted interpolation at the face's actual
position - see *face_velocity_location*.
``"max"``: maximum of the two adjacent cell WSEs.
face_velocity_location : {"centroid", "normal_intercept"}, optional
Position used as the face normal velocity measurement point when
``wse_interp="sloped"``.
``"normal_intercept"`` (default): the point where the
cell-centre connecting line crosses the face polyline, matching
how HEC-RAS locates its finite-difference gradient.
``"centroid"``: the geometric centroid of the face polyline.
Has no effect when ``wse_interp`` is ``"average"`` or ``"max"``.
Returns
-------
ndarray, shape ``(n_cells, 2)``
``[Vx, Vy]`` depth-averaged velocity components when
``component="vector"``.
ndarray, shape ``(n_cells,)``
Velocity magnitude when ``component="speed"`` or flow direction
(degrees clockwise from north) when ``component="angle"``.
Raises
------
KeyError
If ``method="flow_ratio"`` and ``Face Flow`` is absent.
"""
from .velocity import compute_all_cell_velocities
if method == "flow_ratio" and self.face_flow is None:
raise KeyError(
"Face Flow is not present in this HDF file. "
"Enable 'Face Flow' in HEC-RAS HDF5 Write Parameters "
"before running the simulation, or use a different method."
)
face_normal_velocity = np.array(self.face_velocity[timestep, :])
# Read all rows (real + ghost) so boundary face WSE benefits from
# ghost-cell WSE (boundary condition stage).
cell_wse = np.array(self.water_surface[timestep, :])
face_flow = (
np.array(self.face_flow[timestep, :]) if method == "flow_ratio" else None
)
cell_face_info, cell_face_values = self.cell_face_info
face_ae_info, face_ae_values = self.face_area_elevation
if wse_interp == "sloped":
# Stack real + ghost cell coordinates so the sloped face-WSE
# estimator can distance-weight boundary faces correctly.
cell_centers = np.vstack([self.cell_centers, self.ghost_cell_centers])
face_velocity_coords = (
self.face_normal_intercept
if face_velocity_location == "normal_intercept"
else self.face_centroids
)
else:
cell_centers = None
face_velocity_coords = None
vecs = compute_all_cell_velocities(
n_cells=self.n_cells,
cell_face_info=cell_face_info,
cell_face_values=cell_face_values,
face_normals=self.face_normals,
face_cell_indexes=self.face_cell_indexes,
face_ae_info=face_ae_info,
face_ae_values=face_ae_values,
face_normal_velocity=face_normal_velocity,
cell_wse=cell_wse,
method=method,
face_flow=face_flow,
wse_interp=wse_interp,
cell_centers=cell_centers,
face_velocity_coords=face_velocity_coords,
)
if component == "vector":
return vecs
vx = vecs[:, 0]
vy = vecs[:, 1]
speed = np.sqrt(vx**2 + vy**2)
if component == "speed":
return speed
# component == "angle"
angle = (90.0 - np.degrees(np.arctan2(vy, vx))) % 360.0
angle[speed < 1e-10] = np.nan
return angle
@overload
def get_face_velocity(
self, *, timestep: None = ..., face: int, component: Literal["normal"] = ...
) -> pd.Series: ...
@overload
def get_face_velocity(
self, *, timestep: int, face: None = ..., component: Literal["normal"] = ...
) -> np.ndarray: ...
@overload
def get_face_velocity(
self, *, timestep: int, face: None = ..., component: Literal["vector"]
) -> np.ndarray: ...
def get_face_velocity(
self,
*,
timestep: int | None = None,
face: int | None = None,
component: Literal["normal", "vector"] = "normal",
) -> np.ndarray | pd.Series:
"""Face-normal or reconstructed 2-D face velocity.
Parameters
----------
timestep : int, optional
0-based index into the time dimension.
face : int, optional
0-based face index.
component : {"normal", "vector"}, optional
``"normal"`` (default): signed stored face-normal velocity.
``"vector"``: reconstructed ``[Vx, Vy]`` via the RASMapper-exact
C-stencil least-squares pipeline (Steps A + 2).
Returns
-------
pd.Series
Signed face-normal velocity over time indexed by
:attr:`timestamps`, named ``"Face Velocity"``, when
``component="normal", face=f``.
ndarray, shape ``(n_faces,)``
Snapshot of stored face-normal velocities when
``component="normal", timestep=t``.
ndarray, shape ``(n_faces, 2)``
Reconstructed ``[Vx, Vy]`` at every face when
``component="vector", timestep=t``.
Raises
------
ValueError
If both *timestep* and *face* are ``None``.
If ``component="vector"`` and *face* is given (reconstructed
vector history is not supported).
"""
if timestep is None and face is None:
raise ValueError("At least one of timestep or face must be specified.")
if component == "vector":
if face is not None:
raise ValueError(
"component='vector' does not support face= selection; "
"reconstructed vector history is not available. "
"Use timestep= for a whole-mesh snapshot."
)
if timestep is None:
raise ValueError("timestep is required when component='vector'.")
from rivia.geo import _rasmapper_pipeline as _rasmap
cell_wse = np.array(self.water_surface[timestep, :])
face_normal_vel = np.array(self.face_velocity[timestep, :])
cell_face_info, cell_face_values = self.cell_face_info
_cell_face_count = cell_face_info[:, 1].astype(np.int32)
_, _, face_hconn = _rasmap.compute_face_wss(
cell_wse, self._cell_min_elevation, self.face_min_elevation,
self.face_cell_indexes, _cell_face_count,
)
face_connected = (
(face_hconn >= _rasmap.HC_BACKFILL)
& (face_hconn <= _rasmap.HC_DOWNHILL_SHALLOW)
)
face_vel_A, _ = _rasmap.reconstruct_face_velocities(
face_normal_vel, self.face_normals[:, :2],
face_connected, self.face_cell_indexes,
cell_face_info, cell_face_values,
)
return face_vel_A
# component == "normal"
if face is not None and timestep is None:
values = np.array(self.face_velocity[:, face])
return pd.Series(values, index=self.timestamps, name="Face Velocity")
# timestep given (with or without face)
return np.array(self.face_velocity[timestep, :])
@overload
def get_face_flow(
self,
*,
face: int,
timestep: None = ...,
source: Literal["stored", "derived"] = ...,
) -> pd.Series: ...
@overload
def get_face_flow(
self,
*,
face: Iterable[int],
timestep: None = ...,
source: Literal["stored", "derived"] = ...,
) -> pd.DataFrame: ...
@overload
def get_face_flow(
self,
*,
timestep: int,
face: None = ...,
source: Literal["stored", "derived"] = ...,
) -> np.ndarray: ...
def get_face_flow(
self,
*,
timestep: int | None = None,
face: int | Iterable[int] | None = None,
source: Literal["stored", "derived"] = "derived",
) -> pd.Series | pd.DataFrame | np.ndarray:
"""Volumetric flow through one or more faces, or a whole-mesh snapshot.
Parameters
----------
timestep : int, optional
0-based index into the time dimension. When given with
*face* omitted, returns a snapshot array for all faces.
face : int or iterable of int, optional
0-based face index, or an iterable of face indices. When an
iterable is supplied the return type is
:class:`~pandas.DataFrame` (columns = face indices).
source : {"derived", "stored"}, optional
``"derived"`` (default): computes ``flow = wetted_area(WS) *
face_velocity`` at every output timestep, replicating
``GetFaceFlow_PostProcessed``/``GetFaceWSELTimeSeries`` from
``RasMapperLib/RASD2FlowArea.cs``. Per-face water surface is
resolved with the hydraulic-connectivity algorithm RAS uses
(``FaceWSMethod.Ben``,
:func:`~rivia.geo._rasmapper_pipeline.compute_face_wss`).
When multiple faces are requested, :func:`compute_face_wss` is
called once per timestep (not once per timestep per face).
``"stored"``: thin slice of the stored ``Face Flow`` dataset.
Raises ``KeyError`` if that optional dataset is absent.
Returns
-------
pd.Series
Volumetric flow indexed by :attr:`timestamps`, named
``"Face Flow"``, when *face* is a single integer.
pd.DataFrame
Volumetric flow indexed by :attr:`timestamps`, columns = face
indices, when *face* is an iterable.
ndarray, shape ``(n_faces,)``
Whole-mesh flow snapshot when ``timestep=t, face=None``.
Raises
------
ValueError
If both *timestep* and *face* are ``None``.
KeyError
If ``source="stored"`` and the ``Face Flow`` dataset is absent.
Notes
-----
RAS forces ``MinWSPlotTolerance = 0`` for this computation, whereas
:func:`compute_face_wss` always uses ``0.001`` (model units). This
produces a negligible sub-millimetre discrepancy at the dry/wet
threshold.
"""
if timestep is None and face is None:
raise ValueError("At least one of timestep or face must be specified.")
# Whole-mesh snapshot
if timestep is not None and face is None:
if source == "stored":
if self.face_flow is None:
raise KeyError(
"Face Flow is not present in this HDF file. "
"Enable 'Face Flow' in HEC-RAS HDF5 Write Parameters "
"before running the simulation, or use source='derived'."
)
return np.array(self.face_flow[timestep, :])
raise NotImplementedError(
"Whole-mesh snapshot with source='derived' is not yet supported. "
"Specify face= to request derived flow for particular faces."
)
# Time-series path (face= given)
if isinstance(face, Iterable):
face_idxs: list[int] = list(face)
multi = True
else:
face_idxs = [face] # type: ignore[list-item]
multi = False
if source == "stored":
if self.face_flow is None:
raise KeyError(
"Face Flow is not present in this HDF file. "
"Enable 'Face Flow' in HEC-RAS HDF5 Write Parameters "
"before running the simulation, or use source='derived'."
)
if multi:
cols = [np.array(self.face_flow[:, f]) for f in face_idxs]
return pd.DataFrame(
np.column_stack(cols), index=self.timestamps, columns=face_idxs,
)
values = np.array(self.face_flow[:, face_idxs[0]])
return pd.Series(values, index=self.timestamps, name="Face Flow")
from rivia.geo import _rasmapper_pipeline as _rasmap
_NODATA = -9999.0 # FaceValues sentinel (matches _rasmap._NODATA)
cell_wse = np.array(self.water_surface[:, :]) # (n_t, n_cells + n_ghost)
all_face_vel = np.array(self.face_velocity[:, :]) # (n_t, n_faces)
face_normal_vel = all_face_vel[:, face_idxs] # (n_t, n_req)
cell_face_info, _ = self.cell_face_info
_cell_face_count = cell_face_info[:, 1].astype(np.int32)
cellB = self.face_cell_indexes[face_idxs, 1]
face_is_perimeter = _cell_face_count[cellB] == 1 # (n_req,)
face_lengths = self.face_lengths
n_t = cell_wse.shape[0]
n_req = len(face_idxs)
value_a = np.empty((n_t, n_req), dtype=np.float32)
value_b = np.empty((n_t, n_req), dtype=np.float32)
hconn = np.empty((n_t, n_req), dtype=np.uint8)
for t in range(n_t):
fa, fb, fh = _rasmap.compute_face_wss(
cell_wse[t], self._cell_min_elevation, self.face_min_elevation,
self.face_cell_indexes, _cell_face_count,
)
value_a[t] = fa[face_idxs]
value_b[t] = fb[face_idxs]
hconn[t] = fh[face_idxs]
ws = np.where(
value_a == _NODATA, value_b,
np.where(value_b != _NODATA, np.maximum(value_a, value_b), value_a),
)
active = (hconn != _rasmap.HC_NONE) | face_is_perimeter[np.newaxis, :]
valid = active & (ws != _NODATA)
flow = np.zeros((n_t, n_req), dtype=np.float64)
for i, f in enumerate(face_idxs):
face_valid = valid[:, i]
if not face_valid.any():
continue
elevations, areas = self.face_elevation_area(f)
face_length = float(face_lengths[f])
w = ws[face_valid, i].astype(np.float64)
area = np.empty(w.shape, dtype=np.float64)
below = w <= elevations[0]
above = w >= elevations[-1]
mid = ~below & ~above
area[below] = areas[0]
area[above] = areas[-1] + (w[above] - elevations[-1]) * face_length
area[mid] = np.interp(w[mid], elevations, areas)
flow[face_valid, i] = area * face_normal_vel[face_valid, i]
if multi:
return pd.DataFrame(flow, index=self.timestamps, columns=face_idxs)
return pd.Series(flow[:, 0], index=self.timestamps, name="Face Flow")
@overload
def get_facepoint_velocity(
self, *, timestep: None = ..., facepoint: int
) -> pd.DataFrame: ...
@overload
def get_facepoint_velocity(
self, *, timestep: int, facepoint: None = ...
) -> np.ndarray: ...
@overload
def get_facepoint_velocity(
self, *, timestep: int, facepoint: int
) -> np.ndarray: ...
def get_facepoint_velocity(
self,
*,
timestep: int | None = None,
facepoint: int | None = None,
) -> pd.DataFrame | np.ndarray:
"""Local WLS velocity at one or all facepoints.
Replicates RASMapper's lightweight ad-hoc query
``ComputeFacePointVelocity_FacePerpLeastSquares_Weighted_Local``
(``RasMapperLib/MeshFV2D.cs``) - a local weighted-least-squares fit
of the *raw stored* signed face-normal velocities (:attr:`face_velocity`,
not the reconstructed ``[Vx, Vy]`` field) at the faces meeting a
facepoint, weighted by each face's inverse plan-view length. This
is the algorithm behind RASMapper's facepoint context-menu velocity
plot - distinct from, and cheaper/less accurate than, the
RASMapper-exact pipeline used by :meth:`get_facepoint_velocity_field`.
Parameters
----------
timestep : int, optional
0-based index into the time dimension. When given without
*facepoint*, runs the local WLS at every facepoint for that
single timestep (expensive for large meshes).
facepoint : int, optional
0-based facepoint (mesh corner) index. When given without
*timestep*, returns the full time history for that facepoint.
Returns
-------
pd.DataFrame
Columns ``['vx', 'vy', 'speed']``, indexed by :attr:`timestamps`,
when ``facepoint=fp`` (and *timestep* is omitted).
ndarray, shape ``(n_facepoints, 3)``
Columns ``[vx, vy, speed]`` for every facepoint when
``timestep=t`` (and *facepoint* is omitted).
ndarray, shape ``(3,)``
``[vx, vy, speed]`` for a single facepoint at a single timestep
when both are given.
Raises
------
ValueError
If both *timestep* and *facepoint* are ``None``.
Notes
-----
The weighting matrix depends only on mesh geometry (face normals and
lengths), so for the time-history path it is built once and reused
across all timesteps.
"""
if timestep is None and facepoint is None:
raise ValueError(
"At least one of timestep or facepoint must be specified."
)
fp_face_info, fp_face_values = self.facepoint_face_orientation
def _local_wls_history(fp: int) -> pd.DataFrame:
"""Full time-history WLS for a single facepoint."""
start = int(fp_face_info[fp, 0])
count = int(fp_face_info[fp, 1])
faces = fp_face_values[start : start + count, 0].astype(np.int64)
normals = self.face_normals[faces, :2]
# Guard against degenerate zero-length faces.
inv_lengths = 1.0 / np.maximum(self.face_lengths[faces], 1e-12)
# h5py fancy-indexing requires strictly increasing, non-repeated indices.
unique_faces, inverse = np.unique(faces, return_inverse=True)
raw_velocities = np.array(self.face_velocity[:, unique_faces])
velocities = raw_velocities[:, inverse] # (n_t, k)
# FaceVelocityCoef: geometry-only normal-equations matrix + 2x2 inverse.
m11 = float(np.sum(normals[:, 0] ** 2 * inv_lengths))
m22 = float(np.sum(normals[:, 1] ** 2 * inv_lengths))
m12 = float(np.sum(normals[:, 0] * normals[:, 1] * inv_lengths))
det = m11 * m22 - m12 * m12
inv_matrix = (
np.eye(2) / count if det == 0.0
else np.array([[m22, -m12], [-m12, m11]]) / det
)
weighted_normals = normals * inv_lengths[:, None] # (k, 2)
rhs = (weighted_normals[None, :, :] * velocities[:, :, None]).sum(axis=1)
result = rhs @ inv_matrix
speed = np.linalg.norm(result, axis=1)
return pd.DataFrame(
{"vx": result[:, 0], "vy": result[:, 1], "speed": speed},
index=self.timestamps,
)
def _local_wls_snapshot(fp: int, t: int) -> np.ndarray:
"""Single-timestep WLS for one facepoint; returns [vx, vy, speed]."""
start = int(fp_face_info[fp, 0])
count = int(fp_face_info[fp, 1])
faces = fp_face_values[start : start + count, 0].astype(np.int64)
normals = self.face_normals[faces, :2]
inv_lengths = 1.0 / np.maximum(self.face_lengths[faces], 1e-12)
unique_faces, inverse = np.unique(faces, return_inverse=True)
raw_vel = np.array(self.face_velocity[t, unique_faces])
velocities = raw_vel[inverse] # (k,)
m11 = float(np.sum(normals[:, 0] ** 2 * inv_lengths))
m22 = float(np.sum(normals[:, 1] ** 2 * inv_lengths))
m12 = float(np.sum(normals[:, 0] * normals[:, 1] * inv_lengths))
det = m11 * m22 - m12 * m12
inv_matrix = (
np.eye(2) / count if det == 0.0
else np.array([[m22, -m12], [-m12, m11]]) / det
)
rhs = (normals * inv_lengths[:, None] * velocities[:, None]).sum(axis=0)
result = inv_matrix @ rhs
speed = float(np.linalg.norm(result))
return np.array([result[0], result[1], speed])
if facepoint is not None and timestep is None:
return _local_wls_history(facepoint)
if timestep is not None and facepoint is None:
n_fp = len(fp_face_info)
out = np.zeros((n_fp, 3), dtype=np.float64)
for fp in range(n_fp):
out[fp] = _local_wls_snapshot(fp, timestep)
return out
# Both given: single facepoint, single timestep
return _local_wls_snapshot(facepoint, timestep) # type: ignore[arg-type]
def get_facepoint_velocity_field(self, timestep: int) -> np.ndarray:
"""Full 2D velocity ``[Vx, Vy]`` at each mesh facepoint (corner).
Implements the RASMapper-exact pipeline (Steps A + 2 + 3):
* **Step A** -- hydraulic connectivity
(:func:`~rivia.geo._rasmap.compute_face_wss`).
* **Step 2** -- C-stencil face velocity reconstruction
(:func:`~rivia.geo._rasmap.reconstruct_face_velocities`).
* **Step 3** -- inverse-face-length weighted facepoint averaging
(:func:`~rivia.geo._rasmap.compute_facepoint_velocities`).
Each facepoint has one arc-context velocity vector per adjacent
face; these are averaged to produce a single ``[Vx, Vy]`` per
facepoint.
Replicates ``ComputeVertexVelocities`` from
``RasMapperLib/MeshFV2D.cs``.
Parameters
----------
timestep : int
0-based index into the time dimension.
Returns
-------
ndarray, shape ``(n_facepoints, 2)``
``[Vx, Vy]`` velocity at each mesh corner.
Facepoints adjacent only to dry faces receive ``[0, 0]``.
"""
from rivia.geo import _rasmapper_pipeline as _rasmap
cell_wse = np.array(self.water_surface[timestep, :])
face_normal_vel = np.array(self.face_velocity[timestep, :])
cell_face_info, cell_face_values = self.cell_face_info
_cell_face_count = cell_face_info[:, 1].astype(np.int32)
face_value_a, face_value_b, face_hconn = _rasmap.compute_face_wss(
cell_wse, self._cell_min_elevation, self.face_min_elevation,
self.face_cell_indexes, _cell_face_count,
)
face_connected = (
(face_hconn >= _rasmap.HC_BACKFILL)
& (face_hconn <= _rasmap.HC_DOWNHILL_SHALLOW)
)
face_vel_A, face_vel_B = _rasmap.reconstruct_face_velocities(
face_normal_vel, self.face_normals[:, :2],
face_connected, self.face_cell_indexes,
cell_face_info, cell_face_values,
)
fp_face_info, fp_face_values = self.facepoint_face_orientation
fp_vel_data, _ = _rasmap.compute_facepoint_velocities(
face_vel_A, face_vel_B, face_connected,
self.face_lengths,
self.face_facepoint_indexes, self.face_cell_indexes,
cell_wse, fp_face_info, fp_face_values,
face_value_a, face_value_b,
)
# fp_vel_data is a flat CSR array (total_fp_face_entries, 2).
# Use fp_face_info to slice per-facepoint entries and average them.
n_fp = len(fp_face_info)
result = np.zeros((n_fp, 2), dtype=np.float64)
for fp in range(n_fp):
start = int(fp_face_info[fp, 0])
count = int(fp_face_info[fp, 1])
if count > 0:
result[fp] = fp_vel_data[start : start + count].mean(axis=0)
return result
# ------------------------------------------------------------------
# Visualization
# ------------------------------------------------------------------
def plot_velocity(
self,
timestep: int,
cell_index: int,
*,
render_mode: Literal["horizontal", "sloping", "hybrid"] = "sloping",
buffer: int = 1,
reference_raster: str | Path | None = None,
use_depth_weights: bool = False,
shallow_to_flat: bool = False,
pixel_size: float | None = None,
n_arrows: int = 200,
ax: Any | None = None,
) -> Any:
"""Quiver plot of rasterized velocity vectors around a target cell.
Rasterizes the neighbourhood of *cell_index* using the RASMapper-exact
pipeline (``rasterize_results`` with ``variable="velocity"``), then draws
quiver arrows at the pixel centres of wet pixels. Arrows show the
final, fully interpolated ``[Vx, Vy]`` value at each pixel -- the same
values that appear in RASMapper's velocity map.
Mesh polygon outlines and cell index labels are drawn as context.
Requires ``matplotlib`` (``pip install matplotlib``).
Parameters
----------
timestep:
0-based time index.
cell_index:
Target cell. The neighbourhood is expanded by BFS from this cell.
render_mode:
``"horizontal"``, ``"sloping"`` (default), or ``"hybrid"`` --
passed to :meth:`export_raster`.
buffer:
Number of face-adjacency hops to expand from *cell_index*.
``1`` = immediate neighbours; ``2`` = two rings out.
reference_raster:
Optional path to a terrain DEM GeoTIFF. When supplied, the
pixel size and CRS are inherited from the DEM and
``use_depth_weights=True`` becomes available. When ``None``
(default), the pixel size is auto-derived from the median face
length in the neighbourhood.
use_depth_weights:
Passed to :meth:`export_raster`. ``hybrid`` mode only.
Requires *reference_raster*.
shallow_to_flat:
Passed to :meth:`export_raster`. ``hybrid`` mode only.
pixel_size:
Raster pixel size in model coordinate units. Smaller values
produce a finer grid and more potential arrow positions.
Ignored when *reference_raster* is supplied (the DEM pixel size
is used instead). Defaults to ``local_cell_size / 3``, giving
~9 pixels per face-length unit so *n_arrows* has room to work.
n_arrows:
Target number of quiver arrows. Wet pixels are subsampled with
stride ``ceil(sqrt(n_wet / n_arrows))`` so the actual count is
approximately *n_arrows*. Increase for denser plots (e.g.
``n_arrows=800``); set to a very large number to show every
wet pixel. Default ``200``.
ax:
Existing ``matplotlib.Axes`` to draw on. If ``None`` a new
figure/axes is created.
Returns
-------
matplotlib.axes.Axes
"""
try:
import matplotlib.pyplot as plt
except ImportError as exc:
raise ImportError(
"matplotlib is required for plot_velocity(). "
"Install it with: pip install matplotlib"
) from exc
from rivia.geo import raster as _raster
# -- 1. BFS neighbourhood (with ring tracking) -------------------
cell_face_info, cell_face_values = self.cell_face_info
face_cell_indexes = self.face_cell_indexes
n_cells = self.n_cells
# ring_of[c] = BFS hop distance from cell_index (0 = focus cell)
ring_of: dict[int, int] = {cell_index: 0}
neighbors: set[int] = {cell_index}
frontier: set[int] = {cell_index}
for hop in range(buffer):
next_frontier: set[int] = set()
for c in frontier:
start = int(cell_face_info[c, 0])
count = int(cell_face_info[c, 1])
for k in range(count):
fi = int(cell_face_values[start + k, 0])
for nb in (
int(face_cell_indexes[fi, 0]),
int(face_cell_indexes[fi, 1]),
):
if 0 <= nb < n_cells and nb not in neighbors:
ring_of[nb] = hop + 1
next_frontier.add(nb)
neighbors |= next_frontier
frontier = next_frontier
# -- 2. Bounding box and auto cell_size --------------------------
cell_polys = self.cell_polygons
cell_centers = self.cell_centers
all_verts = np.vstack([
cell_polys[c] for c in neighbors if len(cell_polys[c]) > 0
])
x_min = float(all_verts[:, 0].min())
x_max = float(all_verts[:, 0].max())
y_min = float(all_verts[:, 1].min())
y_max = float(all_verts[:, 1].max())
# Local cell size: median face length in the neighbourhood.
# Always computed -- used for arrow sizing and (when no reference_raster)
# as the output pixel size.
nb_face_idx: list[int] = []
for c in neighbors:
start = int(cell_face_info[c, 0])
count = int(cell_face_info[c, 1])
for k in range(count):
nb_face_idx.append(int(cell_face_values[start + k, 0]))
face_lengths = self.face_normals[:, 2]
local_cell_size = float(np.median(face_lengths[sorted(set(nb_face_idx))]))
if reference_raster is None:
# pixel_size overrides the auto value; default is local_cell_size / 3
# so each cell contains ~9 pixels and n_arrows has room to work.
_cell_size: float | None = (
pixel_size if pixel_size is not None else local_cell_size / 3.0
)
else:
_cell_size = None # reference_raster provides the pixel grid
_margin = local_cell_size
# Rectangular perimeter with a small margin
bbox_perim = np.array([
[x_min - _margin, y_min - _margin],
[x_max + _margin, y_min - _margin],
[x_max + _margin, y_max + _margin],
[x_min - _margin, y_max + _margin],
])
# -- 3. Rasterize velocity via full rasmap pipeline ---------------
# reference_raster enables pixel-level dry masking (WSE < terrain ->
# nodata) in addition to the coarser cell-level wet check.
# "velocity_vector" -> 4-band [Vx, Vy, speed, direction]; quiver
# needs all four bands.
fp_face_info, fp_face_values = self.facepoint_face_orientation
ds = _raster.rasterize_results(
variable="velocity_vector",
cell_wse=np.array(self.water_surface[timestep, :]),
cell_min_elevation=self._cell_min_elevation,
face_min_elevation=self.face_min_elevation,
face_cell_indexes=face_cell_indexes,
cell_face_info=cell_face_info,
cell_face_values=cell_face_values,
face_facepoint_indexes=self.face_facepoint_indexes,
fp_coords=self.facepoint_coordinates,
face_normals=self.face_normals,
fp_face_info=fp_face_info,
fp_face_values=fp_face_values,
cell_polygons=cell_polys,
face_normal_velocity=np.array(self.face_velocity[timestep, :]),
output_path=None,
cell_centers=cell_centers,
cell_surface_area=self.cell_surface_area,
reference_raster=reference_raster,
cell_size=_cell_size,
render_mode=render_mode,
use_depth_weights=use_depth_weights,
shallow_to_flat=shallow_to_flat,
tight_extent=False,
perimeter=bbox_perim,
)
# -- 4. Read bands and pixel coordinates -------------------------
from matplotlib.colors import Normalize
from matplotlib.path import Path as _MPath
nodata_val = float(ds.nodata)
vx_grid = ds.read(1).astype(np.float64)
vy_grid = ds.read(2).astype(np.float64)
speed_grid = ds.read(3).astype(np.float64)
transform = ds.transform
ds.close()
height, width = vx_grid.shape
rows_idx, cols_idx = np.mgrid[0:height, 0:width]
qx_grid = transform.c + (cols_idx + 0.5) * transform.a
qy_grid = transform.f + (rows_idx + 0.5) * transform.e # e < 0
wet = (speed_grid > 0) & (speed_grid != nodata_val)
# Assign each pixel to the innermost neighbourhood ring it belongs to.
# Pixels outside all neighbourhood cells are excluded from quiver.
pts = np.column_stack([qx_grid.ravel(), qy_grid.ravel()])
pixel_ring = np.full(height * width, buffer + 1, dtype=np.int32)
for c in neighbors:
poly = cell_polys[c]
if len(poly) < 3:
continue
ring = ring_of.get(c, buffer)
inside = _MPath(poly).contains_points(pts)
# Lower ring number = closer to focus cell -> higher priority
pixel_ring = np.where(inside & (ring < pixel_ring), ring, pixel_ring)
pixel_ring = pixel_ring.reshape(height, width)
wet = wet & (pixel_ring <= buffer)
# -- 5. Arrow scaling and subsampling ----------------------------
# Arrow length is mapped linearly from [sp_min, sp_max] speed to
# [0.30, 0.85] * local_cell_size so every arrow is visible.
# Outer rings are scaled down slightly to emphasise the focus area.
wet_r, wet_c = np.where(wet)
if len(wet_r) == 0:
# No wet pixels in neighbourhood -- draw polygons only
if ax is None:
_, ax = plt.subplots()
else:
speeds = speed_grid[wet_r, wet_c]
rings = pixel_ring[wet_r, wet_c]
sp_min, sp_max = float(speeds.min()), float(speeds.max())
arrow_min = local_cell_size * 0.30
arrow_max = local_cell_size * 0.85
if sp_max > sp_min + 1e-12:
t = (speeds - sp_min) / (sp_max - sp_min)
arrow_len = arrow_min + t * (arrow_max - arrow_min)
else:
arrow_len = np.full(len(speeds), arrow_max)
# Outer rings get 15 % shorter arrows per hop from the focus cell
ring_weight = np.maximum(0.55, 1.0 - rings * 0.15)
arrow_len *= ring_weight
# Unit direction * scaled length (data-coordinate arrows)
eps = 1e-12
u_norm = np.where(speeds > eps, vx_grid[wet_r, wet_c] / speeds, 0.0)
v_norm = np.where(speeds > eps, vy_grid[wet_r, wet_c] / speeds, 0.0)
u_draw = u_norm * arrow_len
v_draw = v_norm * arrow_len
# Subsample: target ~200 arrows (subsampled from neighbourhood only)
stride = max(1, int(np.ceil(np.sqrt(len(wet_r) / n_arrows))))
if stride > 1:
sel = np.arange(0, len(wet_r), stride)
wet_r, wet_c = wet_r[sel], wet_c[sel]
u_draw, v_draw, speeds = u_draw[sel], v_draw[sel], speeds[sel]
norm = Normalize(vmin=sp_min, vmax=sp_max)
# -- 6. Draw -------------------------------------------------
if ax is None:
fig, ax = plt.subplots()
else:
fig = ax.get_figure()
for c in neighbors:
poly = cell_polys[c]
if len(poly) == 0:
continue
is_target = c == cell_index
ring_x = np.append(poly[:, 0], poly[0, 0])
ring_y = np.append(poly[:, 1], poly[0, 1])
fc = "steelblue" if is_target else "lightgray"
lw = 1.8 if is_target else 0.8
ax.fill(poly[:, 0], poly[:, 1], fc=fc, alpha=0.18, ec="none")
ax.plot(ring_x, ring_y, color="dimgray", lw=lw)
cx, cy = float(cell_centers[c, 0]), float(cell_centers[c, 1])
ax.text(cx, cy, str(c), ha="center", va="center",
fontsize=7, color="black", clip_on=True)
Q = ax.quiver(
qx_grid[wet_r, wet_c], qy_grid[wet_r, wet_c],
u_draw, v_draw, speeds,
cmap="plasma", norm=norm,
scale=1.0, scale_units="xy", angles="xy",
pivot="middle",
)
fig.colorbar(Q, ax=ax, label="Speed")
ax.set_aspect("equal")
ax.set_title(
f"{self.name} t={timestep} cell={cell_index} "
f"mode={render_mode} buf={buffer}"
)
ax.set_xlabel("X")
ax.set_ylabel("Y")
return ax
# Fallback: no wet pixels -- still draw polygons
if ax is None:
_, ax = plt.subplots()
for c in neighbors:
poly = cell_polys[c]
if len(poly) == 0:
continue
is_target = c == cell_index
ring_x = np.append(poly[:, 0], poly[0, 0])
ring_y = np.append(poly[:, 1], poly[0, 1])
fc = "steelblue" if is_target else "lightgray"
lw = 1.8 if is_target else 0.8
ax.fill(poly[:, 0], poly[:, 1], fc=fc, alpha=0.18, ec="none")
ax.plot(ring_x, ring_y, color="dimgray", lw=lw)
cx, cy = float(cell_centers[c, 0]), float(cell_centers[c, 1])
ax.text(cx, cy, str(c), ha="center", va="center",
fontsize=7, color="black", clip_on=True)
ax.set_aspect("equal")
ax.set_title(
f"{self.name} t={timestep} cell={cell_index} "
f"mode={render_mode} buf={buffer} [dry]"
)
ax.set_xlabel("X")
ax.set_ylabel("Y")
return ax
# ------------------------------------------------------------------
# Raster export -- delegates to rivia.geo (deferred import)
# ------------------------------------------------------------------
@log_call(logging.INFO)
@timed()
def export_raster(
self,
variable: Literal["wse", "water_surface", "depth", "velocity", "velocity_vector"],
timestep: int | None = None,
output_path: str | Path | None = None,
*,
reference_raster: str | Path | None = None,
cell_size: float | None = None,
crs: Any | None = None,
nodata: float = -9999.0,
render_mode: Literal["horizontal", "sloping", "hybrid"] = "sloping",
use_depth_weights: bool = False,
shallow_to_flat: bool = False,
depth_threshold: float = 0.001,
tight_extent: bool = True,
) -> Path | rasterio.io.DatasetReader:
"""Rasterize a result variable using the RASMapper-exact algorithm.
Implements the pixel-perfect pipeline reverse-engineered from
``RasMapperLib/`` (decompiled C# source, HEC-RAS 6.6),
validated against RASMapper VRT exports -- median ``diff`` = 0.000000.
Parameters
----------
variable:
``"wse"`` / ``"water_surface"`` -- water-surface elevation.
``"depth"`` -- water depth (WSE minus terrain); requires *reference_raster*.
``"velocity"`` -- 1-band speed raster ``sqrt(Vx^2+Vy^2)``; requires an explicit *timestep*.
``"velocity_vector"`` -- 4-band raster ``[Vx, Vy, speed, direction_deg]``; requires an explicit *timestep*.
timestep:
0-based time index. Pass ``None`` to use the time of maximum
water-surface elevation (``"wse"``/``"water_surface"`` and
``"depth"`` only; raises ``ValueError`` for ``"velocity"``).
output_path:
Destination ``.tif`` file path. ``None`` returns an open
in-memory ``rasterio.DatasetReader``; the caller must close it.
reference_raster:
Existing GeoTIFF whose transform and CRS are inherited.
Also used as the terrain DEM for depth computation.
**Required** when ``variable="depth"``.
Mutually exclusive with *cell_size*.
cell_size:
Output pixel size in model coordinate units. Used when no
*reference_raster* is supplied; grid origin is derived from
the flow-area perimeter bounding box.
Mutually exclusive with *reference_raster*.
crs:
Output CRS. Inherited from *reference_raster* when ``None``.
nodata:
Fill value for dry / out-of-domain pixels (default ``-9999``).
render_mode:
``"sloping"`` (default) -- RASMapper "Sloping Cell Corners";
uses corner facepoints only. RasMapperLib hardcodes
``shallow_to_flat=True`` for this mode; the user-supplied value
is overridden. Matches ``store_map(render_mode="sloping")``.
``"hybrid"`` -- "Sloping Cell Corners + Face Centers";
``use_depth_weights`` and ``shallow_to_flat`` are honoured.
Matches ``store_map(render_mode="hybrid")``.
``"horizontal"`` -- flat per-cell value; facepoint interpolation
is skipped. Matches ``store_map(render_mode="horizontal")``.
use_depth_weights:
Weight face contributions by water depth. **``hybrid`` only**;
forced ``False`` for other modes. Requires *reference_raster*.
shallow_to_flat:
Render cells with no hydraulically-connected faces flat.
**``hybrid`` only** (user-configurable); forced ``True`` for
``"sloping"`` per RasMapperLib, ``False`` for ``"horizontal"``.
depth_threshold:
Minimum depth for a pixel to be considered wet (default
``0.001``). Matches ``RASResults.MinWSPlotTolerance``.
tight_extent:
When ``True`` (default), pixels outside the flow-area boundary
polygon are set to *nodata*.
Returns
-------
Path
Written GeoTIFF path when *output_path* is given.
rasterio.io.DatasetReader
Open in-memory dataset when *output_path* is ``None``.
Raises
------
ImportError
If ``rasterio`` or ``shapely`` are not installed.
ValueError
If ``variable="depth"`` and *reference_raster* is not provided.
If ``variable="velocity"`` or ``variable="velocity_vector"`` and ``timestep=None``.
If neither *reference_raster* nor *cell_size* is provided.
"""
from rivia.geo import raster as _raster
if variable in ("velocity", "velocity_vector") and timestep is None:
raise ValueError(
"timestep=None is not supported for velocity. "
"Provide an explicit timestep index."
)
if reference_raster is None and cell_size is None:
# Default: median face length (same heuristic as export_raster)
cell_size = float(np.median(self.face_normals[:, 2]))
# ---- Read HDF arrays ------------------------------------------------
if timestep is None:
cell_wse = self._max_water_surface["value"].to_numpy()
else:
cell_wse = self.get_water_surface(timestep=timestep, include_ghost=True)
face_normal_velocity: np.ndarray | None = None
if variable in ("velocity", "velocity_vector"):
face_normal_velocity = np.array(self.face_velocity[timestep, :])
cell_face_info, cell_face_values = self.cell_face_info
fp_face_info, fp_face_values = self.facepoint_face_orientation
# ---- Delegate to rasterize_results -------------------------------------
return _raster.rasterize_results(
variable=variable,
cell_wse=cell_wse,
cell_min_elevation=self._cell_min_elevation,
face_min_elevation=self.face_min_elevation,
face_cell_indexes=self.face_cell_indexes,
cell_face_info=cell_face_info,
cell_face_values=cell_face_values,
face_facepoint_indexes=self.face_facepoint_indexes,
fp_coords=self.facepoint_coordinates,
face_normals=self.face_normals,
fp_face_info=fp_face_info,
fp_face_values=fp_face_values,
cell_polygons=self.cell_polygons,
face_normal_velocity=face_normal_velocity,
output_path=output_path,
cell_centers=self.cell_centers,
cell_surface_area=self.cell_surface_area,
reference_raster=reference_raster,
cell_size=cell_size,
crs=crs,
nodata=nodata,
render_mode=render_mode,
use_depth_weights=use_depth_weights,
shallow_to_flat=shallow_to_flat,
depth_threshold=depth_threshold,
tight_extent=tight_extent,
perimeter=self.perimeter,
)
@log_call(logging.INFO)
@timed()
def export_hydraulic_rasters(
self,
timestep: int,
reference_raster: str | Path,
*,
wse_path: str | Path | None = None,
depth_path: str | Path | None = None,
velocity_path: str | Path | None = None,
nodata: float = -9999.0,
render_mode: Literal["horizontal", "sloping", "hybrid"] = "sloping",
use_depth_weights: bool = False,
shallow_to_flat: bool = False,
depth_threshold: float = 0.001,
tight_extent: bool = True,
) -> dict[str, Path | "rasterio.io.DatasetReader"]:
"""Export water-surface elevation, depth, and velocity rasters in one call.
Convenience wrapper that calls :meth:`export_raster` three times -- once
each for ``"water_surface"``, ``"depth"``, and ``"velocity"`` -- sharing
the same render settings.
Parameters
----------
timestep:
0-based time index. Required -- all three outputs need a specific
timestep (velocity has no max-value fallback).
reference_raster:
Path to the terrain DEM GeoTIFF. Required -- used to derive depth
(WSE minus DEM) and to inherit the output CRS and transform.
wse_path:
Destination ``.tif`` for the water-surface elevation raster.
``None`` returns an open in-memory ``rasterio.DatasetReader``.
depth_path:
Destination ``.tif`` for the depth raster.
``None`` returns an open in-memory ``rasterio.DatasetReader``.
velocity_path:
Destination ``.tif`` for the velocity magnitude raster.
``None`` returns an open in-memory ``rasterio.DatasetReader``.
nodata:
Fill value for dry / out-of-domain pixels (default ``-9999``).
render_mode:
Passed to :meth:`export_raster` for all three variables.
use_depth_weights:
Passed to :meth:`export_raster`. ``hybrid`` only.
shallow_to_flat:
Passed to :meth:`export_raster`. ``hybrid`` only.
depth_threshold:
Minimum depth for a pixel to be considered wet (default ``0.001``).
tight_extent:
When ``True`` (default), pixels outside the flow-area boundary are
set to *nodata*.
Returns
-------
dict with keys ``"water_surface"``, ``"depth"``, ``"velocity"``.
Each value is the written ``Path`` (when the corresponding ``*_path``
argument is given) or an open in-memory ``rasterio.DatasetReader``
(when ``None``). The caller must close any in-memory datasets.
"""
common = dict(
timestep=timestep,
reference_raster=reference_raster,
nodata=nodata,
render_mode=render_mode,
use_depth_weights=use_depth_weights,
shallow_to_flat=shallow_to_flat,
depth_threshold=depth_threshold,
tight_extent=tight_extent,
)
return {
"water_surface": self.export_raster(
"water_surface", output_path=wse_path, **common
),
"depth": self.export_raster(
"depth", output_path=depth_path, **common
),
"velocity": self.export_raster(
"velocity", output_path=velocity_path, **common
),
}
def flow_across_line(
self,
xy: np.ndarray,
*,
method: Literal["walk", "shortest_path"] = "shortest_path",
) -> pd.Series:
"""Total volumetric discharge through a user-supplied profile line.
Identifies the mesh face "fence" that best approximates *xy* via
:meth:`faces_along_line`, then sums oriented face fluxes
``face_area * face_velocity`` across that fence at every timestep,
returning a time-series of discharge.
The sign convention matches RASMapper: flow from left bank to right
bank (when facing from *xy* start to *xy* end) is **positive**.
Positive flow direction is ``rotate_ccw_90(xy[-1] - xy[0])``.
Parameters
----------
xy : ndarray, shape ``(n_pts, 2)``
Profile polyline vertices ``(x, y)`` drawn from left bank to
right bank.
method : {"shortest_path", "walk"}, optional
Face-selection method passed to :meth:`faces_along_line`.
Default is ``"shortest_path"``.
Returns
-------
pd.Series
Index: :attr:`timestamps` (datetime). Values: signed volumetric
discharge ``(m^3/s or ft^3/s depending on project units)`` across
the profile fence at each timestep. Name is ``"discharge"``.
Raises
------
ValueError
If the polyline does not intersect the mesh or no connected face
path can be found.
NotImplementedError
If ``method="walk"``.
See Also
--------
faces_along_line : identifies which faces are on the fence
get_face_flow : per-face volumetric flux time-series
Notes
-----
Each face contributes ``+face_flow`` when its stored normal agrees
with the positive-flow direction and ``-face_flow`` when it opposes
it. ``get_face_flow`` already incorporates the face area and the
correct RAS sign, so no projection onto the profile line is needed:
the face-normal velocity is the full flux component across the face.
"""
xy = np.asarray(xy, dtype=np.float64)
faces_df = self.faces_along_line(xy, method=method)
face_ids = faces_df["face"].tolist()
orientations = faces_df["orientation"].to_numpy(dtype=bool) # True -> negate
flow_df: pd.DataFrame = self.get_face_flow( # type: ignore[assignment]
face=face_ids, source="derived"
)
signs = np.where(orientations, -1.0, 1.0)
discharge = (flow_df.values * signs[np.newaxis, :]).sum(axis=1)
return pd.Series(discharge, index=flow_df.index, name="discharge")
def wse_along_line(
self,
xy: np.ndarray,
*,
timestep: int | Literal["max"] = "max",
interval: float = 1.0,
) -> pd.DataFrame:
"""Water-surface elevation sampled along a polyline.
Parameters
----------
xy : ndarray, shape ``(n_pts, 2)``
Polyline vertices ``(x, y)`` in model coordinates.
timestep : int or "max", optional
0-based time index, or ``"max"`` (default) to use the
simulation-wide maximum WSE per cell from the summary HDF.
interval : float, optional
Along-line spacing (model units) between regular sample
stations. Default ``1.0``. Cell boundary crossings are
always included regardless of this value.
Returns
-------
pd.DataFrame
Columns:
* ``station`` -- cumulative along-line distance (model units)
from the polyline start.
* ``cell`` -- 0-based mesh cell index at that station.
``-1`` for stations that fall outside the mesh.
* ``wse`` -- water-surface elevation (model units).
``NaN`` for stations outside the mesh or in dry cells
(HEC-RAS stores ``-9999`` for dry cells; converted to
``NaN`` on output).
Notes
-----
Sample stations come from two sources, merged and sorted:
1. Regular interval -- ``np.arange(0, line_length, interval)``
plus the line endpoint.
2. Cell boundary crossings -- the ``station_start`` /
``station_end`` values from :meth:`cells_along_line`
(mandatory, so every cell transition is represented
regardless of *interval*).
Points within ``1e-6`` model units of each other are collapsed
to one station.
No ground-elevation column is included; terrain-pixel crossings
and 8-stencil intra-cell interpolation are deferred to
``geo.profile`` (full RASMapper-equivalent profile).
"""
xy = np.asarray(xy, dtype=np.float64)
if xy.ndim != 2 or xy.shape[1] != 2 or len(xy) < 2:
raise ValueError("xy must be shape (n_pts, 2) with n_pts >= 2.")
cells_df = self.cells_along_line(xy)
if cells_df.empty:
return pd.DataFrame(columns=["station", "cell", "wse"])
# Line length
diffs = np.diff(xy, axis=0)
line_length = float(np.sum(np.hypot(diffs[:, 0], diffs[:, 1])))
# Station points: regular interval + mandatory cell boundary crossings
interval_stations = np.append(
np.arange(0.0, line_length, float(interval)), line_length
)
boundary_stations = np.concatenate([
cells_df["station_start"].to_numpy(),
cells_df["station_end"].to_numpy(),
])
all_sorted = np.sort(np.concatenate([interval_stations, boundary_stations]))
# Collapse near-duplicates within 1e-6 model units
if len(all_sorted) > 1:
keep = np.concatenate([[True], np.diff(all_sorted) > 1e-6])
stations = all_sorted[keep]
else:
stations = all_sorted
# WSE array (real cells only)
_NODATA = -9999.0
if timestep == "max":
cell_wse = self.max_water_surface["value"].to_numpy(dtype=np.float64)
else:
cell_wse = self.get_water_surface(
timestep=int(timestep)
).astype(np.float64)
# Vectorised binary search: for each station find its row in cells_df
station_starts = cells_df["station_start"].to_numpy()
station_ends = cells_df["station_end"].to_numpy()
cell_indices = cells_df["cell"].to_numpy(dtype=np.int64)
idxs = np.searchsorted(station_starts, stations, side="right") - 1
clamped = np.maximum(idxs, 0)
in_range = (idxs >= 0) & (stations <= station_ends[clamped] + 1e-10)
cells_out = np.where(in_range, cell_indices[clamped], np.int64(-1))
# WSE lookup; use index 0 as a safe dummy for out-of-range stations
safe_cells = np.where(in_range, cells_out, 0)
wse_raw = cell_wse[safe_cells]
wse_out = np.where(in_range & (wse_raw != _NODATA), wse_raw, np.nan)
return pd.DataFrame({"station": stations, "cell": cells_out, "wse": wse_out})
# ---------------------------------------------------------------------------
# FlowAreaResults - direct HDF reads; inherits derived methods above
# ---------------------------------------------------------------------------
[docs]
class FlowAreaResults(_FlowAreaResultsDerived):
"""Geometry *and* time-series results for one named 2-D flow area.
Inherits all geometry properties from :class:`FlowArea` and all
derived/computed methods from :class:`_FlowAreaResultsDerived`.
**Method naming conventions**
*Raw dataset properties* -- ``h5py.Dataset``; slice to control what is loaded::
area.water_surface[t] # one timestep -> ndarray (n_cells + n_ghost,)
area.water_surface[a:b] # slice -> ndarray (b-a, n_cells + n_ghost)
area.water_surface[:] # all -> ndarray (n_t, n_cells + n_ghost)
*Snapshot accessors* -- one time, all locations -> ndarray::
area.get_water_surface(timestep=t) # (n_cells,)
area.get_depth(timestep=t) # (n_cells,)
area.get_cell_velocity(t) # (n_cells, 2)
area.get_face_velocity(timestep=t, component="vector") # (n_faces, 2)
area.get_facepoint_velocity_field(t) # (n_facepoints, 2)
*Location time-series* -- one location, all times -> pandas::
area.get_water_surface(cell=c) # Series WSE over time
area.get_face_velocity(face=f) # Series velocity over time
area.get_face_flow(face=f) # Series flow over time
area.get_facepoint_velocity(facepoint=fp) # DataFrame [vx, vy, speed]
*Summary properties* -- aggregate over the full time span, all locations::
area.max_water_surface # DataFrame ['value', 'time'] per cell
area.max_face_velocity # DataFrame ['value', 'time'] per face
area.get_max_depth() # DataFrame ['value', 'time'] per cell
Parameters
----------
geom_group:
``h5py.Group`` at ``Geometry/2D Flow Areas/<name>``.
ts_group:
``h5py.Group`` at the time-series result path for this area.
sum_group:
``h5py.Group`` at the summary result path for this area.
name:
Flow area name.
n_cells:
Number of real computational cells.
"""
def __init__(
self,
geom_group: "h5py.Group",
ts_group: "h5py.Group",
sum_group: "h5py.Group",
name: str,
n_cells: int,
) -> None:
super().__init__(geom_group, name, n_cells)
self._ts = ts_group
self._sum = sum_group
def __repr__(self) -> str:
return (
f"FlowAreaResults({self.name!r},"
f" cells={self.n_cells}, faces={self.n_faces})"
)
# ------------------------------------------------------------------
# Internal HDF helper
# ------------------------------------------------------------------
def _load_summary(self, key: str, n: int | None = None) -> pd.DataFrame:
"""Load a ``(2, n_*)`` summary dataset as a tidy DataFrame.
The HDF summary datasets have shape ``(2, n_elements)`` where
``[0, :]`` = maximum/minimum values and ``[1, :]`` = elapsed-time
(in days) at which the extremum occurred. HEC-RAS stores entries
for ghost cells as well; pass *n* to clip to the first *n* entries.
Returns a DataFrame with columns ``['value', 'time']`` and
integer index corresponding to cell or face index.
"""
raw = np.array(self._sum[key]) # shape (2, n_elements)
if n is not None:
raw = raw[:, :n]
return pd.DataFrame(
{"value": raw[0], "time": raw[1]},
)
# ------------------------------------------------------------------
# Lazy time-series (h5py.Dataset - slice to control memory)
# ------------------------------------------------------------------
@property
def water_surface(self) -> "h5py.Dataset":
"""Water-surface elevation time series.
``h5py.Dataset``, shape ``(n_timesteps, n_cells + n_ghost)``.
HEC-RAS stores ghost cell WSE (boundary condition stages) in the
trailing columns. Slice with ``[:self.n_cells]`` for real cells only,
or ``[:]`` for all including ghost cells.
Slice to read: ``area.water_surface[10]``.
"""
return self._ts["Water Surface"]
@property
def face_velocity(self) -> "h5py.Dataset":
"""Signed face-normal velocity time series.
``h5py.Dataset``, shape ``(n_timesteps, n_faces)``.
"""
return self._ts["Face Velocity"]
@property
def face_flow(self) -> "h5py.Dataset | None":
"""Volumetric face-flow time series, or ``None`` if not output.
``h5py.Dataset``, shape ``(n_timesteps, n_faces)``.
"""
return self._ts.get("Face Flow")
@property
def cell_velocity(self) -> "h5py.Dataset | None":
"""HEC-RAS cell-velocity *speed* scalar, or ``None`` if not output.
``h5py.Dataset``, shape ``(n_timesteps, n_cells)``.
This is the optional output enabled in the HDF Write Parameters;
see :meth:`get_cell_velocity` for the derived vector field.
"""
return self._ts.get("Cell Velocity")
# ------------------------------------------------------------------
# Timestamps
# ------------------------------------------------------------------
@property
def timestamps(self) -> pd.DatetimeIndex:
"""Result output time stamps as a ``pd.DatetimeIndex``.
Parsed from the ``Time Date Stamp`` dataset (sibling of the
``2D Flow Areas`` group under ``Unsteady Time Series``). Used as
the index for all time-series results returned by ``get_*`` accessors.
"""
ds = self._ts.parent.parent["Time Date Stamp"]
raw = np.array(ds).astype(str)
return _parse_hec_ts_array(raw, _RAS_TS_FMT)
# ------------------------------------------------------------------
# Eager summary results (small arrays, loaded once per access)
# ------------------------------------------------------------------
@property
def max_water_surface(self) -> pd.DataFrame:
"""Maximum WSE per cell.
DataFrame with columns ``['value', 'time']``.
``value``: maximum water-surface elevation (model units).
``time``: elapsed simulation time (days) when max occurred.
Index: 0-based cell index. Real cells only (ghost rows excluded).
"""
return self._load_summary("Maximum Water Surface", n=self.n_cells)
@property
def _max_water_surface(self) -> pd.DataFrame:
"""Maximum WSE including ghost cell rows. Same layout as
:attr:`max_water_surface` but shape ``(n_cells + n_ghost,)``.
Required when indexing with raw ``face_cell_indexes`` values.
"""
return self._load_summary("Maximum Water Surface")
@property
def min_water_surface(self) -> pd.DataFrame:
"""Minimum WSE per cell. Same column layout as :attr:`max_water_surface`.
Real cells only (ghost rows excluded).
"""
return self._load_summary("Minimum Water Surface", n=self.n_cells)
@property
def _min_water_surface(self) -> pd.DataFrame:
"""Minimum WSE including ghost cell rows. Same layout as
:attr:`min_water_surface` but shape ``(n_cells + n_ghost,)``.
Required when indexing with raw ``face_cell_indexes`` values.
"""
return self._load_summary("Minimum Water Surface")
@property
def max_face_velocity(self) -> pd.DataFrame:
"""Maximum face velocity per face.
DataFrame with columns ``['value', 'time']``.
Index: 0-based face index.
"""
return self._load_summary("Maximum Face Velocity", n=self.n_faces)
# ---------------------------------------------------------------------------
# FlowAreaResultsCollection
# ---------------------------------------------------------------------------
[docs]
class FlowAreaResultsCollection(FlowAreaCollection):
"""Collection of :class:`FlowAreaResults` objects backed by a plan HDF file.
Overrides :class:`FlowAreaCollection` to return ``FlowAreaResults``
instead of plain ``FlowArea`` instances.
"""
def __getitem__(self, name: str) -> FlowAreaResults:
if name not in self._cache:
root = self._hdf.get("Geometry/2D Flow Areas")
if root is None or name not in root:
raise KeyError(
f"2D flow area {name!r} not found. Available: {self.names}"
)
n_cells = self._get_real_cell_count(name)
# Time-series group for this area
ts_path = f"{_TS_2D}/{name}"
if ts_path not in self._hdf:
raise KeyError(
f"No time-series results found for flow area {name!r} "
f"at '{ts_path}'. Has the plan been computed?"
)
# Summary group (may be absent for steady-flow plans)
sum_path = f"{_SUM_2D}/{name}"
sum_group = self._hdf.get(sum_path)
if sum_group is None:
raise KeyError(
f"No summary results found for flow area {name!r} at '{sum_path}'."
)
self._cache[name] = FlowAreaResults(
geom_group=root[name],
ts_group=self._hdf[ts_path],
sum_group=sum_group,
name=name,
n_cells=n_cells,
)
return self._cache[name] # type: ignore[return-value]
# ---------------------------------------------------------------------------
# StorageAreaResults - extends StorageArea geometry with plan results
# ---------------------------------------------------------------------------
[docs]
class StorageAreaResults(StorageArea):
"""Geometry *and* time-series results for one storage area.
Inherits all geometry properties from :class:`~rivia.hdf.StorageArea`
(:attr:`boundary`, :attr:`volume_elevation`, :meth:`volume_at_elevation`, etc.).
Time-series properties return ``pd.Series`` indexed by :attr:`timestamps`.
Parameters
----------
sa:
Parent geometry object whose fields are copied into this instance.
sa_index:
0-based column index of this SA in the flat ``(n_t, n_sa)`` datasets
(``Water Surface``, ``Flow``) stored under ``Storage Areas/``.
ts_sa_group:
``h5py.Group`` at ``-/Unsteady Time Series/Storage Areas``, or ``None``
when the plan has no SA results.
sum_sa_group:
``h5py.Group`` at ``-/Summary Output/Storage Areas``, or ``None``.
timestamps:
Output timestamps for the block this SA was loaded from; used as the
index for all ``pd.Series`` time-series properties.
"""
def __init__(
self,
sa: StorageArea,
sa_index: int,
ts_sa_group: "h5py.Group | None",
sum_sa_group: "h5py.Group | None",
timestamps: pd.DatetimeIndex,
skip_row0: bool = False,
) -> None:
super().__init__(
name=sa.name,
mode=sa.mode,
boundary=sa.boundary,
volume_elevation=sa.volume_elevation,
)
self._i = sa_index
self._ts = ts_sa_group
self._sum = sum_sa_group
self.timestamps = timestamps
self._skip_row0 = skip_row0
# per-SA subgroup: -/Storage Areas/<name>/
self._sub = ts_sa_group.get(sa.name) if ts_sa_group else None
self._cache: dict[str, np.ndarray] = {}
# ------------------------------------------------------------------
# Helpers
# ------------------------------------------------------------------
def _series(self, col: np.ndarray, name: str) -> pd.Series:
"""Wrap a 1-D time-series array as a timestamp-indexed ``pd.Series``."""
return pd.Series(col, index=self.timestamps, name=name)
def _load_flat(self, key: str) -> np.ndarray:
"""Read a plan-level dataset and return this SA's column.
HEC-RAS stores some SA time series as flat ``(n_t, n_sa)`` arrays
directly under the Storage Areas group (e.g. ``"Water Surface"``,
``"Flow"``). ``key`` is the HDF dataset name; ``self._i`` is this
SA's column index within that array. Use :meth:`_load_vars` for
datasets in the per-SA subgroup instead.
"""
if key not in self._cache:
if self._ts is None:
raise KeyError(
f"No time-series results for storage area {self.name!r}. "
"Has the plan been computed?"
)
data = np.array(self._ts[key])
col = data[1:, self._i] if self._skip_row0 else data[:, self._i]
self._cache[key] = col
return self._cache[key]
def _load_vars(self) -> np.ndarray:
"""Load and cache the ``(n_t, 6)`` Storage Area Variables array."""
if "_vars" not in self._cache:
if self._sub is None or "Storage Area Variables" not in self._sub:
raise KeyError(
f"'Storage Area Variables' not found for storage area "
f"{self.name!r}."
)
raw = np.array(self._sub["Storage Area Variables"])
self._cache["_vars"] = raw[1:] if self._skip_row0 else raw
return self._cache["_vars"]
def _load_summary(self, key: str) -> pd.DataFrame:
"""Read a summary dataset for this SA and return a one-row DataFrame.
Summary datasets have shape ``(2, n_sa)``: row 0 is the extreme value
and row 1 is the simulation time (elapsed days) when it occurred. This
method slices column ``self._i`` from both rows and packages them as
``{"value": ..., "time": ...}``. Only available for ``output="mapping"``
(Base Output summary group).
"""
if self._sum is None:
raise KeyError(
f"No summary results for storage area {self.name!r}. "
"Has the plan been computed?"
)
raw = np.array(self._sum[key]) # shape (2, n_sa)
return pd.DataFrame({"value": [float(raw[0, self._i])],
"time": [float(raw[1, self._i])]})
# ------------------------------------------------------------------
# Flat time-series (one value per timestep)
# ------------------------------------------------------------------
@property
def wse(self) -> pd.Series:
"""Water-surface elevation time series, indexed by :attr:`timestamps`."""
return self._series(self._load_flat("Water Surface"), "Water Surface")
@property
def flow(self) -> pd.Series:
"""Net inflow rate (positive = into SA), indexed by :attr:`timestamps`."""
return self._series(self._load_flat("Flow"), "Flow")
# ------------------------------------------------------------------
# Storage Area Variables columns (WSE, flows, area, volume)
# ------------------------------------------------------------------
@property
def inflow_net(self) -> pd.Series:
"""Net inflow rate, indexed by :attr:`timestamps`."""
return self._series(self._load_vars()[:, 1], "Inflow Net")
@property
def inflow(self) -> pd.Series:
"""Total inflow rate (sum of all inflow sources), indexed by :attr:`timestamps`."""
return self._series(self._load_vars()[:, 2], "Inflow")
@property
def outflow(self) -> pd.Series:
"""Total outflow rate (sum of all outflow sinks), indexed by :attr:`timestamps`."""
return self._series(self._load_vars()[:, 3], "Outflow")
@property
def surface_area(self) -> pd.Series:
"""Water-surface area time series (model area units), indexed by :attr:`timestamps`."""
return self._series(self._load_vars()[:, 4], "Surface Area")
@property
def volume(self) -> pd.Series:
"""Stored volume time series (model volume units), indexed by :attr:`timestamps`."""
return self._series(self._load_vars()[:, 5], "Volume")
# ------------------------------------------------------------------
# Connection inflows
# ------------------------------------------------------------------
@property
def connections(self) -> np.ndarray | None:
"""Inflow from each named connection.
Shape ``(n_t, n_conns)``, or ``None`` when no connection data is stored.
Column names are in :attr:`connection_names`.
"""
if "_conns" not in self._cache:
if self._sub is None:
return None
ds = self._sub.get("Connections to Storage Area")
if ds is None:
return None
raw = np.array(ds)
self._cache["_conns"] = raw[1:] if self._skip_row0 else raw
return self._cache["_conns"]
@property
def connection_names(self) -> list[str]:
"""Names of the inflow connection sources (from HDF ``Connections`` attribute).
Falls back to index-based names if the attribute is absent.
"""
if self._sub is None:
return []
ds = self._sub.get("Connections to Storage Area")
if ds is None:
return []
attr = ds.attrs.get("Connections")
if attr is None:
n = ds.shape[1] if ds.ndim > 1 else 1
return [f"connection_{i}" for i in range(n)]
return [_decode(v) for v in attr]
# ------------------------------------------------------------------
# Summary results
# ------------------------------------------------------------------
@property
def max_wse(self) -> pd.DataFrame:
"""Maximum WSE.
DataFrame with columns ``['value', 'time']``.
``value``: maximum WSE in model units.
``time``: elapsed simulation time (days) when maximum occurred.
"""
return self._load_summary("Maximum Water Surface")
@property
def min_wse(self) -> pd.DataFrame:
"""Minimum WSE. Same column layout as :attr:`max_wse`."""
return self._load_summary("Minimum Water Surface")
# ---------------------------------------------------------------------------
# StorageAreaResultsCollection
# ---------------------------------------------------------------------------
[docs]
class StorageAreaResultsCollection(StorageAreaCollection):
"""Collection of :class:`StorageAreaResults` backed by a plan HDF file.
Overrides :class:`~rivia.hdf.StorageAreaCollection` to return
``StorageAreaResults`` with both geometry *and* plan results.
The :attr:`timestamps` property exposes the output block's time axis as a
``pd.DatetimeIndex``; it is also passed into each :class:`StorageAreaResults`
item so per-SA ``pd.Series`` properties share the same index.
Parameters
----------
hdf:
Open ``h5py.File`` handle.
output : {"mapping", "output", "profile", "post_process"}, optional
Which output block to read. Defaults to ``"mapping"`` (Base Output).
"""
def __init__(
self,
hdf: h5py.File,
output: Literal["mapping", "output", "profile", "post_process"] = "mapping",
) -> None:
super().__init__(hdf)
self._output = output
self._timestamps: pd.DatetimeIndex | None = None
@property
def timestamps(self) -> pd.DatetimeIndex:
"""Output timestamps for this block as a ``pd.DatetimeIndex``.
Resolved lazily from the HDF file on first access.
Raises
------
KeyError
If the expected timestamp dataset is absent from the HDF file.
"""
if self._timestamps is None:
if self._output == "mapping":
ds = self._hdf.get(_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_TIME_STAMP_DS}'."
)
elif self._output == "output":
ds = self._hdf.get(_DSS_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_DSS_TIME_STAMP_DS}'."
)
else: # "profile" and "post_process"
ds = self._hdf.get(_DSS_PROF_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_DSS_PROF_TIME_STAMP_DS}'."
)
self._timestamps = _parse_hec_ts_array(
np.array(ds).astype(str), _RAS_TS_FMT
)
return self._timestamps
def _resolve_paths(
self,
) -> tuple[h5py.Group | None, h5py.Group | None, bool]:
"""Return ``(ts_sa_group, sum_sa_group, skip_row0)`` for ``self._output``."""
if self._output == "mapping":
return self._hdf.get(_TS_SA), self._hdf.get(_SUM_SA), False
if self._output == "output":
return self._hdf.get(f"{_DSS_ROOT}/Storage Areas"), None, False
if self._output == "profile":
return self._hdf.get(_DSS_PROF_SA), None, False
if self._output == "post_process":
return self._hdf.get(_POSTPROC_SA), None, True
raise ValueError(
f"output={self._output!r} is not valid; "
"choose 'mapping', 'output', 'profile', or 'post_process'."
)
def _load(self) -> dict[str, StorageAreaResults]: # type: ignore[override]
if self._items is not None:
return self._items # type: ignore[return-value]
if _SA_ROOT not in self._hdf:
self._items = {}
return self._items # type: ignore[return-value]
ts_sa_group, sum_sa_group, skip_row0 = self._resolve_paths()
timestamps = self.timestamps
# Re-read geometry flat arrays (same logic as StorageAreaCollection._load)
root = self._hdf[_SA_ROOT]
attrs = np.array(root["Attributes"])
poly_info = np.array(root["Polygon Info"])
poly_pts = np.array(root["Polygon Points"])
ve_info = np.array(root["Volume Elevation Info"])
ve_vals = np.array(root["Volume Elevation Values"])
items: dict[str, StorageAreaResults] = {}
for i, row in enumerate(attrs):
name = _decode(row["Name"])
mode = _decode(row["Mode"])
start_pt = int(poly_info[i, 0])
n_pts = int(poly_info[i, 1])
boundary = poly_pts[start_pt : start_pt + n_pts].astype(float)
ve_start = int(ve_info[i, 0])
ve_count = int(ve_info[i, 1])
vol_elev = (
ve_vals[ve_start : ve_start + ve_count].astype(float)
if ve_count > 0
else np.empty((0, 2), dtype=float)
)
sa = StorageArea(
name=name, mode=mode, boundary=boundary, volume_elevation=vol_elev
)
items[name] = StorageAreaResults(
sa=sa,
sa_index=i,
ts_sa_group=ts_sa_group,
sum_sa_group=sum_sa_group,
timestamps=timestamps,
skip_row0=skip_row0,
)
self._items = items # type: ignore[assignment]
return items
def __getitem__(self, name: str) -> StorageAreaResults: # type: ignore[override]
items = self._load()
if name not in items:
raise KeyError(
f"Storage area {name!r} not found. Available: {self.names}"
)
return items[name]
# ---------------------------------------------------------------------------
# _StructureResultsMixin - shared HDF result access for all structure types
# ---------------------------------------------------------------------------
class _StructureResultsMixin:
"""Mixin that adds ``Structure Variables`` HDF access to geometry dataclasses.
All four structure result classes inherit this so ``variable_names``,
``flow_total``, ``stage_hw``, ``stage_tw``, ``weir_variables``, and
``flow_gate`` are implemented once and shared.
Concrete subclasses must set ``self._g`` (the result ``h5py.Group``),
``self._cache`` (empty ``dict``), and ``self.timestamps``
(``pd.DatetimeIndex``) in their ``__init__``.
"""
if TYPE_CHECKING:
_g: h5py.Group
_cache: dict[str, np.ndarray]
timestamps: pd.DatetimeIndex
# ------------------------------------------------------------------
# Internal helpers
# ------------------------------------------------------------------
def _load(self, key: str) -> np.ndarray:
"""Read and cache ``self._g[key]`` as a NumPy array.
When ``_skip_row0`` is ``True`` (post_process output), row 0 of each
dataset is the Max WS envelope and is stripped so the remaining rows
align with ``detailed_timestamps``.
Example: ``self._load("Structure Variables")`` returns shape
``(n_t, n_vars)``; ``self._load("Bridge Variables")`` for bridges.
"""
if key not in self._cache:
raw = np.array(self._g[key])
self._cache[key] = raw[1:] if getattr(self, "_skip_row0", False) else raw
return self._cache[key]
def _series(self, col: np.ndarray, name: str) -> pd.Series:
"""Wrap a 1-D time-series array as a timestamp-indexed ``pd.Series``."""
return pd.Series(col, index=self.timestamps, name=name)
def _col_index(self, *candidates: str) -> int:
"""Index of first column whose name contains any candidate (case-insensitive).
Candidates are tried in order; the first column whose lowercased name
contains a lowercased candidate wins. Raises ``KeyError`` if nothing
matches.
Example: ``self._col_index("stage hw", "hw", "headwater")`` tries the
most specific substring first and falls back to broader ones, so it
matches ``"Stage HW"``, ``"Stage HW US"``, or any label containing
``"headwater"`` — whichever appears earliest in the candidate list.
"""
names_lower = [n.lower() for n in self.variable_names]
for cand in candidates:
cand_l = cand.lower()
for i, n in enumerate(names_lower):
if cand_l in n:
return i
raise KeyError(
f"No column matching {candidates!r} in {self.variable_names!r}"
)
# ------------------------------------------------------------------
# Structure Variables
# ------------------------------------------------------------------
@property
def variable_names(self) -> list[str]:
"""Column names from the ``Structure Variables`` ``Variable_Unit`` attribute.
Falls back to ``col_0``, ``col_1``, ... when the attribute is absent.
"""
ds = self._g["Structure Variables"]
attr = ds.attrs.get("Variable_Unit")
if attr is None:
attr = ds.attrs.get("Variables")
if attr is not None:
return [_decode(v[0]) for v in attr]
return [f"col_{i}" for i in range(ds.shape[1])]
@property
def structure_variables(self) -> "h5py.Dataset":
"""All structure variables as a lazy ``h5py.Dataset``, shape ``(n_t, n_vars)``.
Column names are in :attr:`variable_names`.
"""
return self._g["Structure Variables"]
@property
def flow_total(self) -> pd.Series:
"""Total flow through the structure, indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[
:, self._col_index("total flow", "flow")
]
return self._series(col, "Total Flow")
@property
def stage_hw(self) -> pd.Series:
"""Headwater stage (upstream side), indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[
:, self._col_index("stage hw", "hw", "headwater")
]
return self._series(col, "Stage HW")
@property
def stage_tw(self) -> pd.Series:
"""Tailwater stage (downstream side), indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[
:, self._col_index("stage tw", "tw", "tailwater")
]
return self._series(col, "Stage TW")
@property
def flow_gate_total(self) -> pd.Series | None:
"""Sum of flow through all gate groups, indexed by :attr:`timestamps`.
Returns ``None`` when no ``Total Gate Flow`` column is present in
``Structure Variables`` (i.e. the structure has no gate groups).
"""
try:
col = self._col_index("total gate flow", "gate flow")
except KeyError:
return None
return self._series(
self._load("Structure Variables")[:, col], "Total Gate Flow"
)
@property
def flow_weir(self) -> pd.Series | None:
"""Weir overflow component, indexed by :attr:`timestamps`.
Returns ``None`` when no ``Weir Flow`` column is present in
``Structure Variables`` (i.e. the structure has no weir).
"""
try:
col = self._col_index("weir flow")
except KeyError:
return None
return self._series(
self._load("Structure Variables")[:, col], "Weir Flow"
)
# ------------------------------------------------------------------
# Optional datasets (shared by Inline, Lateral, Bridge, SA2DConnection)
# ------------------------------------------------------------------
@property
def weir_variables(self) -> "h5py.Dataset | None":
"""Detailed weir hydraulics time series, or ``None`` if absent."""
return self._g.get("Weir Variables")
def flow_gate(self, gate_number: int) -> "h5py.Dataset":
"""Gate operation dataset for gate *gate_number* (1-based).
Returns a lazy ``h5py.Dataset``, shape ``(n_t, n_vars)``.
Raises
------
KeyError
If the gate number does not exist for this structure.
"""
path = f"Gate Groups/Gate #{gate_number}"
if path not in self._g:
gates_grp = self._g.get("Gate Groups")
available = list(gates_grp.keys()) if gates_grp is not None else []
raise KeyError(
f"Gate #{gate_number} not found. Available: {available}"
)
return self._g[path]
# ---------------------------------------------------------------------------
# SA2DConnectionResults - one connection between two hydraulic areas
# ---------------------------------------------------------------------------
[docs]
class SA2DConnectionResults(_StructureResultsMixin, SA2DConnection):
"""Geometry *and* time-series results for one HEC-RAS SA/2D connection.
Inherits geometry from :class:`~rivia.hdf.SA2DConnection` and shared HDF
result access (``structure_variables``, ``total_flow``, ``stage_hw``,
``stage_tw``, ``weir_variables``, ``flow_gate``) from
:class:`_StructureResultsMixin`.
Parameters
----------
geom:
Geometry object from :class:`~rivia.hdf.StructureCollection`.
group:
``h5py.Group`` at ``-/SA 2D Area Conn/<plan_name>``.
timestamps:
Output timestamps for this block; index for all ``pd.Series`` properties.
"""
def __init__(
self,
geom: SA2DConnection,
group: h5py.Group,
timestamps: pd.DatetimeIndex,
skip_row0: bool = False,
) -> None:
SA2DConnection.__init__(
self,
mode=geom.mode,
upstream_type=geom.upstream_type,
downstream_type=geom.downstream_type,
centerline=geom.centerline,
name=geom.name,
upstream_node=geom.upstream_node,
downstream_node=geom.downstream_node,
)
# Plan result group name may differ from geometry name for 2D<->2D connections
# (HEC-RAS prefixes the flow-area name, e.g. "Lower Levee" ->
# "BaldEagleCr Lower Levee").
self._plan_name: str = group.name.split("/")[-1]
self._g = group
self.timestamps = timestamps
self._skip_row0 = skip_row0
self._cache: dict[str, np.ndarray] = {}
# ------------------------------------------------------------------
# SA2D-specific result datasets
# ------------------------------------------------------------------
@property
def breaching_variables(self) -> "h5py.Dataset | None":
"""Breach geometry and flow time series, or ``None`` if not breach-capable.
Lazy ``h5py.Dataset``, shape ``(n_t, 10)``.
Columns: Stage HW, Stage TW, Bottom Width, Bottom Elevation,
Left Side Slope, Right Side Slope, Breach Flow, Breach Velocity,
Breach Flow Area, Top Elevation.
"""
return self._g.get("Breaching Variables")
@property
def headwater_cells(self) -> np.ndarray | None:
"""2-D mesh cell indices on the headwater side, or ``None`` if absent.
Shape ``(n_faces,)``.
"""
if "Headwater Cells" not in self._g:
return None
return self._load("Headwater Cells")
@property
def tailwater_cells(self) -> np.ndarray | None:
"""2-D mesh cell indices on the tailwater side, or ``None`` if absent.
For 2D<->2D connections (levees) these are stored as a flat ``int32``
dataset at the group root. For SA<->2D connections (e.g. a dam with a
storage-area headwater) they are stored as fixed-width byte strings in
``HW TW Segments/Tailwater Cells`` and are decoded here.
Shape ``(n_cells,)``.
"""
# 2D<->2D: flat int32 at group root
if "Tailwater Cells" in self._g:
return self._load("Tailwater Cells")
# SA<->2D: string-encoded cell indices in HW TW Segments subgroup
seg = self._g.get("HW TW Segments")
if seg is not None and "Tailwater Cells" in seg:
raw = seg["Tailwater Cells"][:]
return np.array(
[int(v.decode().strip()) for v in raw], dtype=np.int32
)
return None
# ---------------------------------------------------------------------------
# InlineResults - inline structure geometry + plan results
# ---------------------------------------------------------------------------
[docs]
class InlineResults(_StructureResultsMixin, InlineStructure):
"""Geometry *and* time-series results for one HEC-RAS inline structure.
Inherits geometry from :class:`~rivia.hdf.InlineStructure` and shared HDF result
access from :class:`_StructureResultsMixin`.
The HDF group is at
``Results/.../DSS Hydrograph Output/.../Inline Structures/<river reach rs>``.
Parameters
----------
geom:
Geometry object from :class:`~rivia.hdf.StructureCollection`.
group:
``h5py.Group`` at the inline structure result path.
timestamps:
Output timestamps for this block; index for all ``pd.Series`` properties.
"""
def __init__(
self, geom: InlineStructure, group: h5py.Group, timestamps: pd.DatetimeIndex
) -> None:
InlineStructure.__init__(
self,
mode=geom.mode,
upstream_type=geom.upstream_type,
downstream_type=geom.downstream_type,
centerline=geom.centerline,
location=geom.location,
upstream_node=geom.upstream_node,
downstream_node=geom.downstream_node,
weir=geom.weir,
gate_groups=geom.gate_groups,
)
self._g = group
self.timestamps = timestamps
self._cache: dict[str, np.ndarray] = {}
# ---------------------------------------------------------------------------
# LateralResults - lateral structure geometry + plan results
# ---------------------------------------------------------------------------
[docs]
class LateralResults(_StructureResultsMixin, LateralStructure):
"""Geometry *and* time-series results for one HEC-RAS lateral structure.
Inherits geometry from :class:`~rivia.hdf.LateralStructure` and shared HDF result
access from :class:`_StructureResultsMixin`.
The HDF group is at
``Results/.../DSS Hydrograph Output/.../Lateral Structures/<river reach rs>``.
``downstream_node`` is the name of the connected Storage Area or 2-D Flow
Area, or an empty string when flow exits the system.
Parameters
----------
geom:
Geometry object from :class:`~rivia.hdf.StructureCollection`.
group:
``h5py.Group`` at the lateral structure result path.
timestamps:
Output timestamps for this block; index for all ``pd.Series`` properties.
"""
def __init__(
self,
geom: LateralStructure,
group: h5py.Group,
timestamps: pd.DatetimeIndex,
skip_row0: bool = False,
) -> None:
LateralStructure.__init__(
self,
mode=geom.mode,
upstream_type=geom.upstream_type,
downstream_type=geom.downstream_type,
centerline=geom.centerline,
location=geom.location,
upstream_node=geom.upstream_node,
downstream_node=geom.downstream_node,
weir=geom.weir,
gate_groups=geom.gate_groups,
)
self._g = group
self.timestamps = timestamps
self._skip_row0 = skip_row0
self._cache: dict[str, np.ndarray] = {}
# ------------------------------------------------------------------
# Lateral-specific reach-split variables
# ------------------------------------------------------------------
@property
def flow_hw_us(self) -> pd.Series:
"""Flow at the upstream bounding cross section, indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[:, self._col_index("flow hw us")]
return self._series(col, "Flow HW US")
@property
def flow_hw_ds(self) -> pd.Series:
"""Flow at the downstream bounding cross section, indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[:, self._col_index("flow hw ds")]
return self._series(col, "Flow HW DS")
@property
def stage_hw_us(self) -> pd.Series:
"""Stage at the upstream bounding cross section, indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[:, self._col_index("stage hw us")]
return self._series(col, "Stage HW US")
@property
def stage_hw_ds(self) -> pd.Series:
"""Stage at the downstream bounding cross section, indexed by :attr:`timestamps`."""
col = self._load("Structure Variables")[:, self._col_index("stage hw ds")]
return self._series(col, "Stage HW DS")
# ---------------------------------------------------------------------------
# BridgeResults - bridge geometry + plan results
# ---------------------------------------------------------------------------
[docs]
class BridgeResults(_StructureResultsMixin, Bridge):
"""Geometry *and* time-series results for one HEC-RAS bridge structure.
Inherits geometry from :class:`~rivia.hdf.Bridge` and shared HDF result
access from :class:`_StructureResultsMixin`.
The HDF group is at
``Results/.../DSS Hydrograph Output/.../Bridge/<river reach rs>``.
Parameters
----------
geom:
Geometry object from :class:`~rivia.hdf.StructureCollection`.
group:
``h5py.Group`` at the bridge result path.
timestamps:
Output timestamps for this block; index for all ``pd.Series`` properties.
"""
def __init__(
self, geom: Bridge, group: h5py.Group, timestamps: pd.DatetimeIndex
) -> None:
Bridge.__init__(
self,
mode=geom.mode,
upstream_type=geom.upstream_type,
downstream_type=geom.downstream_type,
centerline=geom.centerline,
location=geom.location,
upstream_node=geom.upstream_node,
downstream_node=geom.downstream_node,
weir=geom.weir,
gate_groups=geom.gate_groups,
)
self._g = group
self.timestamps = timestamps
self._cache: dict[str, np.ndarray] = {}
# ------------------------------------------------------------------
# Override mixin to use "Bridge Variables" (not "Structure Variables")
# ------------------------------------------------------------------
@property
def variable_names(self) -> list[str]:
"""Column names from the ``Bridge Variables`` ``Variable_Unit`` attribute."""
ds = self._g["Bridge Variables"]
attr = ds.attrs.get("Variable_Unit")
if attr is None:
attr = ds.attrs.get("Variables")
if attr is not None:
return [_decode(v[0]) for v in attr]
return [f"col_{i}" for i in range(ds.shape[1])]
@property
def structure_variables(self) -> "h5py.Dataset":
"""All bridge variables as a lazy ``h5py.Dataset``, shape ``(n_t, 3)``."""
return self._g["Bridge Variables"]
@property
def flow_total(self) -> pd.Series:
"""Total flow through the bridge, indexed by :attr:`timestamps`."""
return self._series(
self._load("Bridge Variables")[:, self._col_index("flow")], "Flow"
)
@property
def stage_hw(self) -> pd.Series:
"""Headwater stage (upstream side), indexed by :attr:`timestamps`."""
return self._series(
self._load("Bridge Variables")[
:, self._col_index("stage hw", "hw", "headwater")
],
"Stage HW",
)
@property
def stage_tw(self) -> pd.Series:
"""Tailwater stage (downstream side), indexed by :attr:`timestamps`."""
return self._series(
self._load("Bridge Variables")[
:, self._col_index("stage tw", "tw", "tailwater")
],
"Stage TW",
)
# ---------------------------------------------------------------------------
# StructureResultsCollection - plan-enriched StructureCollection
# ---------------------------------------------------------------------------
[docs]
class StructureResultsCollection(StructureCollection):
"""Plan-enriched structure collection: all structure types with results.
Overrides :class:`~rivia.hdf.StructureCollection` so each item carries
both geometry attributes *and* time-series result access.
Parameters
----------
hdf:
Open ``h5py.File`` handle.
output : {"mapping", "output", "profile", "post_process"}, optional
Which output block to read. Defaults to ``"output"`` (DSS Hydrograph).
``"mapping"`` — Base Output (all four structure types).
``"output"`` — DSS Hydrograph Output (all four structure types).
``"profile"`` — DSS Profile Output (all four structure types).
``"post_process"`` — Post Process Profiles. Lateral and SA/2D
connections only; Inline and Bridge items fall back to plain geometry.
When no plan result group is found for a structure, the plain geometry
object is kept unchanged.
"""
def __init__(
self,
hdf: h5py.File,
output: Literal["mapping", "output", "profile", "post_process"] = "output",
) -> None:
super().__init__(hdf)
self._output = output
self._timestamps: pd.DatetimeIndex | None = None
@property
def timestamps(self) -> pd.DatetimeIndex:
"""Output timestamps for this block as a ``pd.DatetimeIndex``.
Resolved lazily from the HDF file on first access.
Raises
------
KeyError
If the expected timestamp dataset is absent from the HDF file.
"""
if self._timestamps is None:
if self._output == "mapping":
ds = self._hdf.get(_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_TIME_STAMP_DS}'."
)
elif self._output == "output":
ds = self._hdf.get(_DSS_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_DSS_TIME_STAMP_DS}'."
)
else: # "profile" and "post_process"
ds = self._hdf.get(_DSS_PROF_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_DSS_PROF_TIME_STAMP_DS}'."
)
self._timestamps = _parse_hec_ts_array(
np.array(ds).astype(str), _RAS_TS_FMT
)
return self._timestamps
def _block_paths(
self,
) -> tuple[str, str | None, str, str | None]:
"""Return ``(sa_conn_path, inline_path, lateral_path, bridge_path)``.
``None`` means the structure type is absent from this output block
(e.g. Inline and Bridge are absent from ``"post_process"``).
"""
if self._output == "mapping":
return _TS_SA_CONN, _TS_INLINE, _TS_LATERAL, _TS_BRIDGE
if self._output == "output":
return _DSS_SA_CONN, _DSS_INLINE, _DSS_LATERAL, _DSS_BRIDGE
if self._output == "profile":
return (
_DSS_PROF_SA_CONN, _DSS_PROF_INLINE,
_DSS_PROF_LATERAL, _DSS_PROF_BRIDGE,
)
if self._output == "post_process":
return _POSTPROC_SA_CONN, None, _POSTPROC_LATERAL, None
raise ValueError(
f"output={self._output!r} is not valid; "
"choose 'mapping', 'output', 'profile', or 'post_process'."
)
def _load(self) -> dict[str, Structure]: # type: ignore[override]
if self._items is not None:
return self._items
import h5py as _h5
# Build geometry items first (parent caches in self._items).
geom_items = StructureCollection._load(self)
sa_conn_path, inline_path, lateral_path, bridge_path = self._block_paths()
skip_row0 = self._output == "post_process"
# Helper: collect sub-groups from an HDF path (returns {} when absent).
def _groups(path: str | None) -> dict[str, h5py.Group]:
if path is None:
return {}
root = self._hdf.get(path)
if root is None:
return {}
return {k: v for k, v in root.items() if isinstance(v, _h5.Group)}
conn_groups = _groups(sa_conn_path)
inline_groups = _groups(inline_path)
lateral_groups = _groups(lateral_path)
bridge_groups = _groups(bridge_path)
# Resolve timestamps lazily — only read from HDF when at least one
# result group is found, so calling structures() on a plan that has
# no matching output block (e.g. no DSS Hydrograph data) does not
# raise when all structures fall back to plain geometry objects.
_ts: pd.DatetimeIndex | None = None
def _resolved_ts() -> pd.DatetimeIndex:
nonlocal _ts
if _ts is None:
_ts = self.timestamps
return _ts
items: dict[str, Structure] = {}
for key, geom in geom_items.items():
if isinstance(geom, SA2DConnection):
# Derive plan result group name from geometry fields:
# 2D<->2D (levee): "{upstream_2d_area} {connection}"
# SA<->2D / SA<->SA (one end is SA or '--'): Connection name
if (
geom.upstream_type == "2D"
and geom.downstream_type == "2D"
):
plan_key = f"{geom.upstream_node} {geom.name}"
else:
plan_key = geom.name
grp = conn_groups.get(plan_key)
items[key] = (
SA2DConnectionResults(geom, grp, _resolved_ts(), skip_row0=skip_row0)
if grp is not None else geom
)
elif isinstance(geom, InlineStructure):
plan_key = " ".join(geom.location)
grp = inline_groups.get(plan_key)
items[key] = (
InlineResults(geom, grp, _resolved_ts()) if grp is not None else geom
)
elif isinstance(geom, LateralStructure):
plan_key = " ".join(geom.location)
grp = lateral_groups.get(plan_key)
items[key] = (
LateralResults(geom, grp, _resolved_ts(), skip_row0=skip_row0)
if grp is not None else geom
)
elif isinstance(geom, Bridge):
plan_key = " ".join(geom.location)
grp = bridge_groups.get(plan_key)
items[key] = (
BridgeResults(geom, grp, _resolved_ts()) if grp is not None else geom
)
else:
items[key] = geom
self._items = items
return self._items
# ---------------------------------------------------------------------------
# CrossSectionResults / CrossSectionResultsCollection
# ---------------------------------------------------------------------------
class _CrossSectionResultsBase(CrossSection):
"""Private base for all three cross-section result variants.
Holds the HDF handle, column index, result-group root, and the
timestamps resolved by the parent collection. Provides the shared
``_load`` and ``_series`` helpers.
Concrete subclasses add the properties specific to their output block
(:class:`CrossSectionMappingResults`, :class:`CrossSectionOutputResults`,
:class:`CrossSectionPostProcessResults`).
"""
def __init__(
self,
geom: CrossSection,
hdf: "h5py.File",
index: int,
root: str,
timestamps: pd.DatetimeIndex,
) -> None:
CrossSection.__init__(
self,
river=geom.river,
reach=geom.reach,
rs=geom.rs,
name=geom.name,
left_bank=geom.left_bank,
right_bank=geom.right_bank,
len_left=geom.len_left,
len_channel=geom.len_channel,
len_right=geom.len_right,
contraction=geom.contraction,
expansion=geom.expansion,
station_elevation=geom.station_elevation,
mannings_n=geom.mannings_n,
cut_line=geom.cut_line,
centerline_polyline=geom.centerline_polyline,
)
self._hdf = hdf
self._index = index
self._root = root
self._cache: dict[str, np.ndarray] = {}
self.timestamps = timestamps
def _load(self, dataset: str) -> np.ndarray:
"""Load column ``self._index`` from ``{root}/{dataset}``, cached."""
if dataset not in self._cache:
ds = self._hdf.get(f"{self._root}/{dataset}")
if ds is None:
raise KeyError(
f"Dataset '{dataset}' not found at '{self._root}'."
)
self._cache[dataset] = np.array(ds[:, self._index])
return self._cache[dataset]
def _series(self, dataset: str, name: str) -> pd.Series:
"""Return a time-indexed ``pd.Series`` for one dataset column."""
return pd.Series(self._load(dataset), index=self.timestamps, name=name)
[docs]
class CrossSectionMappingResults(_CrossSectionResultsBase):
"""Geometry *and* results for one XS from the **Base Output** block.
Returned by ``plan.cross_sections("mapping")[key]``. All variables
written by HEC-RAS at the mapping interval are exposed as time-indexed
``pd.Series``.
Parameters
----------
geom:
Geometry object from :class:`CrossSectionCollection`.
hdf:
Open ``h5py.File`` -- kept alive by the parent ``UnsteadyPlan`` context.
index:
Column index of this XS in the ``(n_t, n_xs)`` result datasets.
root:
HDF path prefix -- ``_TS_XS``.
timestamps:
Mapping output timestamps from the parent collection.
"""
@property
def wse(self) -> pd.Series:
"""Water surface elevation time series, indexed by :attr:`timestamps`."""
return self._series("Water Surface", "Water Surface")
@property
def flow(self) -> pd.Series:
"""Flow time series, indexed by :attr:`timestamps`."""
return self._series("Flow", "Flow")
@property
def flow_cumulative(self) -> pd.Series:
"""Cumulative flow volume time series, indexed by :attr:`timestamps`."""
return self._series("Flow Volume Cumulative", "Flow Volume Cumulative")
@property
def flow_lateral(self) -> pd.Series:
"""Lateral flow time series, indexed by :attr:`timestamps`."""
return self._series("Flow Lateral", "Flow Lateral")
@property
def velocity_channel(self) -> pd.Series:
"""Channel velocity time series, indexed by :attr:`timestamps`."""
return self._series("Velocity Channel", "Velocity Channel")
@property
def velocity_total(self) -> pd.Series:
"""Total velocity time series, indexed by :attr:`timestamps`."""
return self._series("Velocity Total", "Velocity Total")
[docs]
class CrossSectionOutputResults(_CrossSectionResultsBase):
"""Geometry *and* results for one XS from the **DSS Hydrograph Output** block.
Returned by ``plan.cross_sections("output")[key]``.
Available datasets: ``wse``, ``flow``, ``flow_cumulative``.
Parameters
----------
geom:
Geometry object from :class:`CrossSectionCollection`.
hdf:
Open ``h5py.File`` -- kept alive by the parent ``UnsteadyPlan`` context.
index:
Column index of this XS in the ``(n_t, n_xs)`` result datasets.
root:
HDF path prefix -- ``_DSS_XS``.
timestamps:
DSS hydrograph output timestamps from the parent collection.
"""
@property
def wse(self) -> pd.Series:
"""Water surface elevation time series, indexed by :attr:`timestamps`."""
return self._series("Water Surface", "Water Surface")
@property
def flow(self) -> pd.Series:
"""Flow time series, indexed by :attr:`timestamps`."""
return self._series("Flow", "Flow")
@property
def flow_cumulative(self) -> pd.Series:
"""Cumulative flow volume time series, indexed by :attr:`timestamps`."""
return self._series("Flow Volume Cumulative", "Flow Volume Cumulative")
[docs]
class CrossSectionProfileResults(_CrossSectionResultsBase):
"""Geometry *and* results for one XS from the **DSS Profile Output** block.
Returned by ``plan.cross_sections("profile")[key]``.
DSS Profile Output contains raw unsteady-engine output at the Detailed
Output Interval. Available datasets: ``wse``, ``flow``. Also exposes
the ``Cross Section Attributes`` structured array unique to this block.
Parameters
----------
geom:
Geometry object from :class:`CrossSectionCollection`.
hdf:
Open ``h5py.File`` -- kept alive by the parent ``UnsteadyPlan`` context.
index:
Column index of this XS in the ``(n_t, n_xs)`` result datasets.
root:
HDF path prefix -- ``_DSS_PROF_XS``.
timestamps:
Detailed output timestamps from the parent collection.
"""
@property
def wse(self) -> pd.Series:
"""Water surface elevation time series, indexed by :attr:`timestamps`."""
return self._series("Water Surface", "Water Surface")
@property
def flow(self) -> pd.Series:
"""Flow time series, indexed by :attr:`timestamps`."""
return self._series("Flow", "Flow")
@property
def cross_section_attributes(self) -> np.void | None:
"""Structured-array row for this XS from ``Cross Section Attributes``.
Contains River, Reach, Station, and Name fields. Returns ``None``
when the dataset is absent from this HDF file.
"""
ds = self._hdf.get(f"{self._root}/Cross Section Attributes")
if ds is None:
return None
return ds[self._index]
def _resolve_inst_dataset(hdf: h5py.File, root: str, variable: str) -> str:
"""Resolve *variable* to an HDF path relative to *root*.
Tries ``{root}/{variable}`` first (top-level: ``Water Surface``, ``Flow``,
``Energy Grade``), then ``{root}/Additional Variables/{variable}``.
Parameters
----------
hdf:
Open ``h5py.File``.
root:
HDF path prefix for the Post Process Profiles XS group.
variable:
Dataset name without subgroup prefix, e.g. ``"Water Surface"``,
``"Velocity Channel"``.
Returns
-------
str
Path component to append to *root* (e.g. ``"Water Surface"`` or
``"Additional Variables/Velocity Channel"``).
Raises
------
KeyError
If *variable* is not found in either location.
"""
if hdf.get(f"{root}/{variable}") is not None:
return variable
av = f"Additional Variables/{variable}"
if hdf.get(f"{root}/{av}") is not None:
return av
raise KeyError(
f"Variable {variable!r} not found at '{root}/{variable}' "
f"or '{root}/{av}'."
)
[docs]
class CrossSectionPostProcessResults(_CrossSectionResultsBase):
"""Geometry *and* results for one XS from the **Post Process Profiles** block.
Returned by ``plan.cross_sections("post_process")[key]``.
Named properties expose every variable as a time-indexed ``pd.Series``
(timeseries only — the Max WS envelope row is excluded). The shape and
index type are identical to :class:`CrossSectionMappingResults` and
:class:`CrossSectionOutputResults`::
xs = plan.cross_sections("post_process")["Butte Cr Upper 7"]
xs = plan.cross_sections("post_process")[0]
xs.wse # pd.Series, index=pd.DatetimeIndex
xs.flow # pd.Series, index=pd.DatetimeIndex
xs.velocity_channel # pd.Series, index=pd.DatetimeIndex
To access the Max WS envelope value for a specific variable, use the
parent collection::
coll = plan.cross_sections("post_process")
coll.wse()["max_wse"].loc[(river, reach, rs)]
For multi-XS access prefer the collection's :meth:`profile_table` or named
methods (one bulk HDF read vs. one read per property call here).
Parameters
----------
geom:
Geometry object from :class:`CrossSectionCollection`.
hdf:
Open ``h5py.File`` -- kept alive by the parent ``UnsteadyPlan`` context.
index:
Column index of this XS in the ``(n_profiles, n_xs)`` result datasets.
root:
HDF path prefix -- ``_POSTPROC_XS``.
timestamps:
Instantaneous profile timestamps from the parent collection
(Max WS excluded, length ``n``).
"""
# ------------------------------------------------------------------
# Private helper
# ------------------------------------------------------------------
def _series_inst(self, variable: str, name: str) -> pd.Series:
"""Read one variable for this XS; return a time-indexed pd.Series.
Skips HDF row 0 (Max WS envelope). One direct column slice per call.
"""
dataset = _resolve_inst_dataset(self._hdf, self._root, variable)
col = np.array(self._hdf[f"{self._root}/{dataset}"][1:, self._index])
return pd.Series(col, index=self.timestamps, name=name)
# ------------------------------------------------------------------
# Named properties — top-level datasets
# ------------------------------------------------------------------
@property
def wse(self) -> pd.Series:
"""Water-surface elevation time series, indexed by :attr:`timestamps`."""
return self._series_inst("Water Surface", "Water Surface")
@property
def flow(self) -> pd.Series:
"""Flow time series, indexed by :attr:`timestamps`."""
return self._series_inst("Flow", "Flow")
@property
def energy_grade(self) -> pd.Series:
"""Energy grade line elevation time series, indexed by :attr:`timestamps`."""
return self._series_inst("Energy Grade", "Energy Grade")
# ------------------------------------------------------------------
# Named properties — Additional Variables
# ------------------------------------------------------------------
@property
def alpha(self) -> pd.Series:
"""Velocity-head correction factor alpha, indexed by :attr:`timestamps`."""
return self._series_inst("Alpha", "Alpha")
@property
def beta(self) -> pd.Series:
"""Momentum correction factor beta, indexed by :attr:`timestamps`."""
return self._series_inst("Beta", "Beta")
@property
def flow_area_channel(self) -> pd.Series:
"""Channel flow area time series (m²), indexed by :attr:`timestamps`."""
return self._series_inst("Area Flow Channel", "Area Flow Channel")
@property
def flow_area_left_ob(self) -> pd.Series:
"""Left overbank flow area time series (m²), indexed by :attr:`timestamps`."""
return self._series_inst("Area Flow Left OB", "Area Flow Left OB")
@property
def flow_area_right_ob(self) -> pd.Series:
"""Right overbank flow area time series (m²), indexed by :attr:`timestamps`."""
return self._series_inst("Area Flow Right OB", "Area Flow Right OB")
@property
def flow_area_total(self) -> pd.Series:
"""Total flow area time series (m²), indexed by :attr:`timestamps`."""
return self._series_inst("Area Flow Total", "Area Flow Total")
@property
def ineffective_area_channel(self) -> pd.Series:
"""Channel area including ineffective zones (m²), indexed by timestamps."""
return self._series_inst(
"Area including Ineffective Channel", "Area including Ineffective Channel"
)
@property
def ineffective_area_left_ob(self) -> pd.Series:
"""Left overbank area incl. ineffective zones (m²), indexed by timestamps."""
return self._series_inst(
"Area including Ineffective Left OB", "Area including Ineffective Left OB"
)
@property
def ineffective_area_right_ob(self) -> pd.Series:
"""Right overbank area incl. ineffective zones (m²), indexed by timestamps."""
return self._series_inst(
"Area including Ineffective Right OB", "Area including Ineffective Right OB"
)
@property
def ineffective_area_total(self) -> pd.Series:
"""Total area incl. ineffective zones (m²), indexed by timestamps."""
return self._series_inst(
"Area including Ineffective Total", "Area including Ineffective Total"
)
@property
def conveyance_channel(self) -> pd.Series:
"""Channel conveyance (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Conveyance Channel", "Conveyance Channel")
@property
def conveyance_left_ob(self) -> pd.Series:
"""Left overbank conveyance (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Conveyance Left OB", "Conveyance Left OB")
@property
def conveyance_right_ob(self) -> pd.Series:
"""Right overbank conveyance (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Conveyance Right OB", "Conveyance Right OB")
@property
def conveyance_total(self) -> pd.Series:
"""Total conveyance (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Conveyance Total", "Conveyance Total")
@property
def critical_energy_grade(self) -> pd.Series:
"""Critical energy grade line elevation (m), indexed by :attr:`timestamps`."""
return self._series_inst("Critical Energy Grade", "Critical Energy Grade")
@property
def critical_water_surface(self) -> pd.Series:
"""Critical water surface elevation (m), indexed by :attr:`timestamps`."""
return self._series_inst("Critical Water Surface", "Critical Water Surface")
@property
def energy_grade_slope(self) -> pd.Series:
"""Energy grade slope (m/m), indexed by :attr:`timestamps`."""
return self._series_inst("EG Slope", "EG Slope")
@property
def friction_slope(self) -> pd.Series:
"""Friction slope (m/m), indexed by :attr:`timestamps`."""
return self._series_inst("Friction Slope", "Friction Slope")
@property
def flow_channel(self) -> pd.Series:
"""Channel flow (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Flow Channel", "Flow Channel")
@property
def flow_left_ob(self) -> pd.Series:
"""Left overbank flow (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Flow Left OB", "Flow Left OB")
@property
def flow_right_ob(self) -> pd.Series:
"""Right overbank flow (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Flow Right OB", "Flow Right OB")
@property
def flow_total(self) -> pd.Series:
"""Total flow (m³/s), indexed by :attr:`timestamps`."""
return self._series_inst("Flow Total", "Flow Total")
@property
def hydraulic_depth_channel(self) -> pd.Series:
"""Channel hydraulic depth (m), indexed by :attr:`timestamps`."""
return self._series_inst("Hydraulic Depth Channel", "Hydraulic Depth Channel")
@property
def hydraulic_depth_left_ob(self) -> pd.Series:
"""Left overbank hydraulic depth (m), indexed by :attr:`timestamps`."""
return self._series_inst("Hydraulic Depth Left OB", "Hydraulic Depth Left OB")
@property
def hydraulic_depth_right_ob(self) -> pd.Series:
"""Right overbank hydraulic depth (m), indexed by :attr:`timestamps`."""
return self._series_inst(
"Hydraulic Depth Right OB", "Hydraulic Depth Right OB"
)
@property
def hydraulic_depth_total(self) -> pd.Series:
"""Total hydraulic depth (m), indexed by :attr:`timestamps`."""
return self._series_inst("Hydraulic Depth Total", "Hydraulic Depth Total")
@property
def hydraulic_radius_channel(self) -> pd.Series:
"""Channel hydraulic radius (m), indexed by :attr:`timestamps`."""
return self._series_inst("Hydraulic Radius Channel", "Hydraulic Radius Channel")
@property
def hydraulic_radius_left_ob(self) -> pd.Series:
"""Left overbank hydraulic radius (m), indexed by :attr:`timestamps`."""
return self._series_inst(
"Hydraulic Radius Left OB", "Hydraulic Radius Left OB"
)
@property
def hydraulic_radius_right_ob(self) -> pd.Series:
"""Right overbank hydraulic radius (m), indexed by :attr:`timestamps`."""
return self._series_inst(
"Hydraulic Radius Right OB", "Hydraulic Radius Right OB"
)
@property
def hydraulic_radius_total(self) -> pd.Series:
"""Total hydraulic radius (m), indexed by :attr:`timestamps`."""
return self._series_inst("Hydraulic Radius Total", "Hydraulic Radius Total")
@property
def mannings_n_channel(self) -> pd.Series:
"""Weighted channel Manning's n, indexed by :attr:`timestamps`."""
return self._series_inst("Manning n Channel", "Manning n Channel")
@property
def mannings_n_left_ob(self) -> pd.Series:
"""Left overbank Manning's n, indexed by :attr:`timestamps`."""
return self._series_inst("Manning n Left OB", "Manning n Left OB")
@property
def mannings_n_right_ob(self) -> pd.Series:
"""Right overbank Manning's n, indexed by :attr:`timestamps`."""
return self._series_inst("Manning n Right OB", "Manning n Right OB")
@property
def mannings_n_total(self) -> pd.Series:
"""Total weighted Manning's n, indexed by :attr:`timestamps`."""
return self._series_inst("Manning n Total", "Manning n Total")
@property
def max_depth_total(self) -> pd.Series:
"""Total maximum water depth (m), indexed by :attr:`timestamps`."""
return self._series_inst("Maximum Depth Total", "Maximum Depth Total")
@property
def shear(self) -> pd.Series:
"""Bed shear stress (N/m²), indexed by :attr:`timestamps`."""
return self._series_inst("Shear", "Shear")
@property
def top_width_channel(self) -> pd.Series:
"""Channel top width (m), indexed by :attr:`timestamps`."""
return self._series_inst("Top Width Channel", "Top Width Channel")
@property
def top_width_channel_with_ineffective(self) -> pd.Series:
"""Channel top width incl. ineffective areas (m), indexed by timestamps."""
return self._series_inst(
"Top Width Channel including Ineffective",
"Top Width Channel including Ineffective",
)
@property
def top_width_left_ob(self) -> pd.Series:
"""Left overbank top width (m), indexed by :attr:`timestamps`."""
return self._series_inst("Top Width Left OB", "Top Width Left OB")
@property
def top_width_left_ob_with_ineffective(self) -> pd.Series:
"""Left OB top width incl. ineffective areas (m), indexed by timestamps."""
return self._series_inst(
"Top Width Left OB including Ineffective",
"Top Width Left OB including Ineffective",
)
@property
def top_width_right_ob(self) -> pd.Series:
"""Right overbank top width (m), indexed by :attr:`timestamps`."""
return self._series_inst("Top Width Right OB", "Top Width Right OB")
@property
def top_width_right_ob_with_ineffective(self) -> pd.Series:
"""Right OB top width incl. ineffective areas (m), indexed by timestamps."""
return self._series_inst(
"Top Width Right OB including Ineffective",
"Top Width Right OB including Ineffective",
)
@property
def top_width_total(self) -> pd.Series:
"""Total top width (m), indexed by :attr:`timestamps`."""
return self._series_inst("Top Width Total", "Top Width Total")
@property
def top_width_total_with_ineffective(self) -> pd.Series:
"""Total top width incl. ineffective areas (m), indexed by timestamps."""
return self._series_inst(
"Top Width Total including Ineffective",
"Top Width Total including Ineffective",
)
@property
def velocity_channel(self) -> pd.Series:
"""Channel velocity (m/s), indexed by :attr:`timestamps`."""
return self._series_inst("Velocity Channel", "Velocity Channel")
@property
def velocity_left_ob(self) -> pd.Series:
"""Left overbank velocity (m/s), indexed by :attr:`timestamps`."""
return self._series_inst("Velocity Left OB", "Velocity Left OB")
@property
def velocity_right_ob(self) -> pd.Series:
"""Right overbank velocity (m/s), indexed by :attr:`timestamps`."""
return self._series_inst("Velocity Right OB", "Velocity Right OB")
@property
def velocity_total(self) -> pd.Series:
"""Total velocity (m/s), indexed by :attr:`timestamps`."""
return self._series_inst("Velocity Total", "Velocity Total")
@property
def wse_total(self) -> pd.Series:
"""Total WSE from Additional Variables (m), indexed by :attr:`timestamps`."""
return self._series_inst("Water Surface Total", "Water Surface Total")
@property
def wetted_perimeter_channel(self) -> pd.Series:
"""Channel wetted perimeter (m), indexed by :attr:`timestamps`."""
return self._series_inst("Wetted Perimeter Channel", "Wetted Perimeter Channel")
@property
def wetted_perimeter_left_ob(self) -> pd.Series:
"""Left overbank wetted perimeter (m), indexed by :attr:`timestamps`."""
return self._series_inst(
"Wetted Perimeter Left OB", "Wetted Perimeter Left OB"
)
@property
def wetted_perimeter_right_ob(self) -> pd.Series:
"""Right overbank wetted perimeter (m), indexed by :attr:`timestamps`."""
return self._series_inst(
"Wetted Perimeter Right OB", "Wetted Perimeter Right OB"
)
@property
def wetted_perimeter_total(self) -> pd.Series:
"""Total wetted perimeter (m), indexed by :attr:`timestamps`."""
return self._series_inst("Wetted Perimeter Total", "Wetted Perimeter Total")
[docs]
class CrossSectionResultsCollection(CrossSectionCollection):
"""Plan-enriched cross section collection with time-series results.
Parameterised over the concrete result class and the HDF path used to
map cross sections to column indices, so one implementation serves all
three output blocks.
Parameters
----------
hdf:
Open ``h5py.File`` handle.
root:
HDF path to the cross section result group:
``_TS_XS``, ``_DSS_XS``, or ``_POSTPROC_XS``.
result_cls:
Concrete result class to instantiate per cross section --
:class:`CrossSectionMappingResults`,
:class:`CrossSectionOutputResults`, or
:class:`CrossSectionPostProcessResults`.
attrs_path:
HDF path to the ``Cross Section Attributes`` structured array used
to map ``(river, reach, station)`` -> column index. Defaults to
``f"{root}/Cross Section Attributes"``, which is correct for Base
Output and DSS blocks. Pass ``_POSTPROC_GEOM_ATTRS`` for the Post
Process block, where attributes live outside the XS result group.
timestamps_fn:
Zero-argument callable that returns the ``pd.DatetimeIndex`` for
this block. Resolved lazily on first access; the result is cached
on :attr:`timestamps`.
"""
def __init__(
self,
hdf: "h5py.File",
root: str,
result_cls: type[_CrossSectionResultsBase] = CrossSectionMappingResults,
attrs_path: str | None = None,
timestamps_fn: Callable[[], pd.DatetimeIndex] | None = None,
) -> None:
super().__init__(hdf)
self._root = root
self._result_cls = result_cls
self._attrs_path = attrs_path or f"{root}/Cross Section Attributes"
self._result_items: dict[str, _CrossSectionResultsBase] | None = None
self._timestamps_fn = timestamps_fn
self._timestamps: pd.DatetimeIndex | None = None
@property
def timestamps(self) -> pd.DatetimeIndex:
"""Timestamps for this result block as a ``pd.DatetimeIndex``.
Resolved lazily from the callable passed at construction time.
Raises
------
AttributeError
If no ``timestamps_fn`` was supplied.
"""
if self._timestamps is None:
if self._timestamps_fn is None:
raise AttributeError(
"No timestamps_fn was supplied for this collection."
)
self._timestamps = self._timestamps_fn()
return self._timestamps
def _load_results(self) -> dict[str, _CrossSectionResultsBase]:
"""Join geometry XS objects to their column index in the result datasets.
HEC-RAS result blocks store XS data as ``(n_t, n_xs)`` arrays where
the column order is given by a ``Cross Section Attributes`` structured
array (River/Reach/Station fields). This method builds a
``(river, reach, station) → column_index`` lookup from that array, then
matches it against the geometry items from
:meth:`~CrossSectionCollection._load`. XS present in geometry but
absent from the result block (e.g. structures written as XS in the
geometry file) are silently excluded from the returned dict.
"""
if self._result_items is not None:
return self._result_items
geom_items = CrossSectionCollection._load(self)
attrs_ds = self._hdf.get(self._attrs_path)
if attrs_ds is None:
self._result_items = {}
return self._result_items
result_attrs = np.array(attrs_ds)
fn = attrs_ds.dtype.names
result_index: dict[tuple[str, str, str], int] = {}
for i, row in enumerate(result_attrs):
r = _decode(row["River"]) if "River" in fn else ""
rc = _decode(row["Reach"]) if "Reach" in fn else ""
st = _decode(row["Station"]) if "Station" in fn else ""
result_index[(r, rc, st)] = i
ts = self.timestamps
items: dict[str, _CrossSectionResultsBase] = {}
for key, geom in geom_items.items():
idx = result_index.get((geom.river, geom.reach, geom.rs))
if idx is not None:
items[key] = self._result_cls(
geom, self._hdf, idx, self._root, ts
)
self._result_items = items
return self._result_items
@overload
def __getitem__(self, key: int) -> _CrossSectionResultsBase: ...
@overload
def __getitem__(self, key: str) -> _CrossSectionResultsBase: ...
@overload
def __getitem__(self, key: tuple[str, str, str]) -> _CrossSectionResultsBase: ...
def __getitem__(
self, key: int | str | tuple[str, str, str]
) -> _CrossSectionResultsBase:
items = self._load_results()
if isinstance(key, int):
keys = list(items)
try:
return items[keys[key]]
except IndexError:
raise IndexError(
f"Index {key} out of range (n={len(items)})"
) from None
if isinstance(key, tuple):
str_key = self._loc_index.get(key)
if str_key is None:
raise KeyError(f"Cross section {key!r} not found.")
if str_key not in items:
raise KeyError(
f"Cross section {key!r} has no results in this plan."
)
return items[str_key]
if key not in items:
raise KeyError(
f"Cross section {key!r} not found. Available: {self.names}"
)
return items[key]
def __len__(self) -> int:
return len(self._load_results())
def __iter__(self) -> Iterator[_CrossSectionResultsBase]:
return iter(self._load_results().values())
@property
def names(self) -> list[str]:
"""Cross section keys available in this result block.
Returns
-------
list[str]
Keys of the form ``"<river> <reach> <station>"``.
"""
return list(self._load_results().keys())
[docs]
class CrossSectionPostProcessResultsCollection(CrossSectionResultsCollection):
"""Cross-section results for the **Post Process Profiles** block.
Returned by ``plan.cross_sections("post_process")``.
All result data is read through this collection rather than through
per-XS objects. The generic engine is :meth:`profile_table`; named
convenience methods (``wse``, ``flow``, ``velocity_channel``, …)
delegate to it.
The profile axis of every returned ``pd.DataFrame`` is labeled
``["max_wse", 0, 1, …, n-1]`` where ``"max_wse"`` is the Max WS
envelope (HDF index 0) and integers ``0 … n-1`` are the instantaneous
profiles aligned with :attr:`~CrossSectionResultsCollection.timestamps`.
Parameters
----------
hdf:
Open ``h5py.File`` handle.
root:
HDF path prefix -- ``_POSTPROC_XS``.
result_cls:
:class:`CrossSectionPostProcessResults` (geometry carrier).
attrs_path:
``_POSTPROC_GEOM_ATTRS`` (attributes live outside the XS group).
timestamps_fn:
Callable returning ``instantaneous_timestamps``.
"""
# ------------------------------------------------------------------
# Engine
# ------------------------------------------------------------------
def _resolve_dataset(self, variable: str) -> str:
"""Resolve *variable* to an HDF path relative to ``self._root``.
Delegates to :func:`_resolve_inst_dataset`.
Parameters
----------
variable:
Dataset name, e.g. ``"Water Surface"``, ``"Velocity Channel"``.
Returns
-------
str
Path component to append to ``self._root``.
Raises
------
KeyError
If *variable* is not found in either location.
"""
return _resolve_inst_dataset(self._hdf, self._root, variable)
[docs]
def profile_table(self, variable: str) -> pd.DataFrame:
"""Return a location × profiles ``pd.DataFrame`` for one variable.
Reads the full ``(n_profiles, n_xs)`` HDF dataset in a single call
and selects only the columns belonging to cross sections present in
this collection.
Parameters
----------
variable:
Dataset name, e.g. ``"Water Surface"``, ``"Velocity Channel"``,
``"Flow Total"``. Resolved with friendly fallback: the literal
name is tried first, then ``Additional Variables/<variable>``.
Returns
-------
pd.DataFrame
* **Index** — ``pd.MultiIndex`` with levels ``(River, Reach, RS)``,
one row per cross section in collection order.
* **Columns** — ``["max_wse", 0, 1, …, n-1]``. ``"max_wse"``
is the Max WS envelope (HDF row 0); integers map 1-to-1 into
:attr:`~CrossSectionResultsCollection.timestamps`.
* **Values** — ``float``.
Raises
------
KeyError
If *variable* is not found in the HDF.
ValueError
If the profile count in the HDF does not match
``len(timestamps) + 1``.
"""
items = self._load_results()
if not items:
return pd.DataFrame()
dataset = self._resolve_dataset(variable)
if not hasattr(self, "_array_cache"):
self._array_cache: dict[str, np.ndarray] = {}
if dataset not in self._array_cache:
self._array_cache[dataset] = np.asarray(
self._hdf[f"{self._root}/{dataset}"]
)
data = self._array_cache[dataset] # (n_profiles, n_xs)
results = list(items.values())
values = data[:, [r._index for r in results]].T # (n_xs, n_profiles)
index = pd.MultiIndex.from_tuples(
[(r.river, r.reach, r.rs) for r in results],
names=["River", "Reach", "RS"],
)
n_ts = len(self.timestamps)
if data.shape[0] == n_ts + 1:
columns: list | pd.Index = ["max_wse", *range(n_ts)]
elif data.shape[0] == n_ts:
columns = self.timestamps
else:
raise ValueError(
f"{dataset!r} has {data.shape[0]} profiles; expected "
f"{n_ts + 1} (instantaneous + Max WS) from timestamps."
)
return pd.DataFrame(values, index=index, columns=columns)
# ------------------------------------------------------------------
# Named accessors -- top-level datasets
# ------------------------------------------------------------------
@property
def wse(self) -> pd.DataFrame:
"""Water-surface elevation, location × profiles (m).
Returns
-------
pd.DataFrame
Columns ``["max_wse", 0, 1, …]``.
"""
return self.profile_table("Water Surface")
@property
def flow(self) -> pd.DataFrame:
"""Flow, location × profiles (m³/s).
Returns
-------
pd.DataFrame
Columns ``["max_wse", 0, 1, …]``.
"""
return self.profile_table("Flow")
@property
def energy_grade(self) -> pd.DataFrame:
"""Energy grade line elevation, location × profiles (m).
Returns
-------
pd.DataFrame
Columns ``["max_wse", 0, 1, …]``.
"""
return self.profile_table("Energy Grade")
# ------------------------------------------------------------------
# Named accessors -- Additional Variables
# ------------------------------------------------------------------
@property
def alpha(self) -> pd.DataFrame:
"""Velocity-head correction factor alpha, location × profiles."""
return self.profile_table("Alpha")
@property
def beta(self) -> pd.DataFrame:
"""Momentum correction factor beta, location × profiles."""
return self.profile_table("Beta")
@property
def flow_area_channel(self) -> pd.DataFrame:
"""Channel flow area, location × profiles (m²)."""
return self.profile_table("Area Flow Channel")
@property
def flow_area_left_ob(self) -> pd.DataFrame:
"""Left overbank flow area, location × profiles (m²)."""
return self.profile_table("Area Flow Left OB")
@property
def flow_area_right_ob(self) -> pd.DataFrame:
"""Right overbank flow area, location × profiles (m²)."""
return self.profile_table("Area Flow Right OB")
@property
def flow_area_total(self) -> pd.DataFrame:
"""Total flow area, location × profiles (m²)."""
return self.profile_table("Area Flow Total")
@property
def ineffective_area_channel(self) -> pd.DataFrame:
"""Channel area including ineffective zones, location × profiles (m²)."""
return self.profile_table("Area including Ineffective Channel")
@property
def ineffective_area_left_ob(self) -> pd.DataFrame:
"""Left overbank area including ineffective zones, location × profiles (m²)."""
return self.profile_table("Area including Ineffective Left OB")
@property
def ineffective_area_right_ob(self) -> pd.DataFrame:
"""Right overbank area including ineffective zones, location × profiles (m²)."""
return self.profile_table("Area including Ineffective Right OB")
@property
def ineffective_area_total(self) -> pd.DataFrame:
"""Total area including ineffective zones, location × profiles (m²)."""
return self.profile_table("Area including Ineffective Total")
@property
def conveyance_channel(self) -> pd.DataFrame:
"""Channel conveyance, location × profiles (m³/s)."""
return self.profile_table("Conveyance Channel")
@property
def conveyance_left_ob(self) -> pd.DataFrame:
"""Left overbank conveyance, location × profiles (m³/s)."""
return self.profile_table("Conveyance Left OB")
@property
def conveyance_right_ob(self) -> pd.DataFrame:
"""Right overbank conveyance, location × profiles (m³/s)."""
return self.profile_table("Conveyance Right OB")
@property
def conveyance_total(self) -> pd.DataFrame:
"""Total conveyance, location × profiles (m³/s)."""
return self.profile_table("Conveyance Total")
@property
def critical_energy_grade(self) -> pd.DataFrame:
"""Critical energy grade line elevation, location × profiles (m)."""
return self.profile_table("Critical Energy Grade")
@property
def critical_water_surface(self) -> pd.DataFrame:
"""Critical water surface elevation, location × profiles (m)."""
return self.profile_table("Critical Water Surface")
@property
def energy_grade_slope(self) -> pd.DataFrame:
"""Energy grade slope, location × profiles (m/m)."""
return self.profile_table("EG Slope")
@property
def friction_slope(self) -> pd.DataFrame:
"""Friction slope, location × profiles (m/m)."""
return self.profile_table("Friction Slope")
@property
def flow_channel(self) -> pd.DataFrame:
"""Channel flow, location × profiles (m³/s)."""
return self.profile_table("Flow Channel")
@property
def flow_left_ob(self) -> pd.DataFrame:
"""Left overbank flow, location × profiles (m³/s)."""
return self.profile_table("Flow Left OB")
@property
def flow_right_ob(self) -> pd.DataFrame:
"""Right overbank flow, location × profiles (m³/s)."""
return self.profile_table("Flow Right OB")
@property
def flow_total(self) -> pd.DataFrame:
"""Total flow, location × profiles (m³/s)."""
return self.profile_table("Flow Total")
@property
def hydraulic_depth_channel(self) -> pd.DataFrame:
"""Channel hydraulic depth, location × profiles (m)."""
return self.profile_table("Hydraulic Depth Channel")
@property
def hydraulic_depth_left_ob(self) -> pd.DataFrame:
"""Left overbank hydraulic depth, location × profiles (m)."""
return self.profile_table("Hydraulic Depth Left OB")
@property
def hydraulic_depth_right_ob(self) -> pd.DataFrame:
"""Right overbank hydraulic depth, location × profiles (m)."""
return self.profile_table("Hydraulic Depth Right OB")
@property
def hydraulic_depth_total(self) -> pd.DataFrame:
"""Total hydraulic depth, location × profiles (m)."""
return self.profile_table("Hydraulic Depth Total")
@property
def hydraulic_radius_channel(self) -> pd.DataFrame:
"""Channel hydraulic radius, location × profiles (m)."""
return self.profile_table("Hydraulic Radius Channel")
@property
def hydraulic_radius_left_ob(self) -> pd.DataFrame:
"""Left overbank hydraulic radius, location × profiles (m)."""
return self.profile_table("Hydraulic Radius Left OB")
@property
def hydraulic_radius_right_ob(self) -> pd.DataFrame:
"""Right overbank hydraulic radius, location × profiles (m)."""
return self.profile_table("Hydraulic Radius Right OB")
@property
def hydraulic_radius_total(self) -> pd.DataFrame:
"""Total hydraulic radius, location × profiles (m)."""
return self.profile_table("Hydraulic Radius Total")
@property
def mannings_n_channel(self) -> pd.DataFrame:
"""Weighted channel Manning's n, location × profiles."""
return self.profile_table("Manning n Channel")
@property
def mannings_n_left_ob(self) -> pd.DataFrame:
"""Left overbank Manning's n, location × profiles."""
return self.profile_table("Manning n Left OB")
@property
def mannings_n_right_ob(self) -> pd.DataFrame:
"""Right overbank Manning's n, location × profiles."""
return self.profile_table("Manning n Right OB")
@property
def mannings_n_total(self) -> pd.DataFrame:
"""Total weighted Manning's n, location × profiles."""
return self.profile_table("Manning n Total")
@property
def max_depth_total(self) -> pd.DataFrame:
"""Total maximum water depth, location × profiles (m)."""
return self.profile_table("Maximum Depth Total")
@property
def shear(self) -> pd.DataFrame:
"""Bed shear stress, location × profiles (N/m²)."""
return self.profile_table("Shear")
@property
def top_width_channel(self) -> pd.DataFrame:
"""Channel top width, location × profiles (m)."""
return self.profile_table("Top Width Channel")
@property
def top_width_channel_with_ineffective(self) -> pd.DataFrame:
"""Channel top width including ineffective areas, location × profiles (m)."""
return self.profile_table("Top Width Channel including Ineffective")
@property
def top_width_left_ob(self) -> pd.DataFrame:
"""Left overbank top width, location × profiles (m)."""
return self.profile_table("Top Width Left OB")
@property
def top_width_left_ob_with_ineffective(self) -> pd.DataFrame:
"""Left overbank top width incl. ineffective areas, location × profiles (m)."""
return self.profile_table("Top Width Left OB including Ineffective")
@property
def top_width_right_ob(self) -> pd.DataFrame:
"""Right overbank top width, location × profiles (m)."""
return self.profile_table("Top Width Right OB")
@property
def top_width_right_ob_with_ineffective(self) -> pd.DataFrame:
"""Right overbank top width incl. ineffective areas, location × profiles (m)."""
return self.profile_table("Top Width Right OB including Ineffective")
@property
def top_width_total(self) -> pd.DataFrame:
"""Total top width, location × profiles (m)."""
return self.profile_table("Top Width Total")
@property
def top_width_total_with_ineffective(self) -> pd.DataFrame:
"""Total top width including ineffective areas, location × profiles (m)."""
return self.profile_table("Top Width Total including Ineffective")
@property
def velocity_channel(self) -> pd.DataFrame:
"""Channel velocity, location × profiles (m/s)."""
return self.profile_table("Velocity Channel")
@property
def velocity_left_ob(self) -> pd.DataFrame:
"""Left overbank velocity, location × profiles (m/s)."""
return self.profile_table("Velocity Left OB")
@property
def velocity_right_ob(self) -> pd.DataFrame:
"""Right overbank velocity, location × profiles (m/s)."""
return self.profile_table("Velocity Right OB")
@property
def velocity_total(self) -> pd.DataFrame:
"""Total velocity, location × profiles (m/s)."""
return self.profile_table("Velocity Total")
@property
def wse_total(self) -> pd.DataFrame:
"""Total WSE from Additional Variables, location × profiles (m)."""
return self.profile_table("Water Surface Total")
@property
def wetted_perimeter_channel(self) -> pd.DataFrame:
"""Channel wetted perimeter, location × profiles (m)."""
return self.profile_table("Wetted Perimeter Channel")
@property
def wetted_perimeter_left_ob(self) -> pd.DataFrame:
"""Left overbank wetted perimeter, location × profiles (m)."""
return self.profile_table("Wetted Perimeter Left OB")
@property
def wetted_perimeter_right_ob(self) -> pd.DataFrame:
"""Right overbank wetted perimeter, location × profiles (m)."""
return self.profile_table("Wetted Perimeter Right OB")
@property
def wetted_perimeter_total(self) -> pd.DataFrame:
"""Total wetted perimeter, location × profiles (m)."""
return self.profile_table("Wetted Perimeter Total")
# ---------------------------------------------------------------------------
# UnsteadyPlan - public entry point
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# Simulation summary dataclasses
# ---------------------------------------------------------------------------
[docs]
@dataclasses.dataclass
class RunStatus:
"""Overall run metadata from ``Results/Unsteady/Summary``.
Attributes
----------
solution:
HEC-RAS solution status string, e.g.
``"Unsteady Finished Successfully"`` or ``"Unsteady Went Unstable"``.
run_window:
Wall-clock window during which the simulation ran, as a raw string,
e.g. ``"26NOV2025 15:58:44 to 26NOV2025 16:03:41"``.
compute_time_total:
Total wall-clock computation time in ``"HH:MM:SS"`` format.
compute_time_dss:
Time spent writing DSS output in ``"HH:MM:SS"`` format.
max_cores:
Maximum number of CPU cores used.
max_wsel_error:
Maximum water-surface elevation error (model units), or ``None``
when the simulation went unstable before convergence.
time_unstable:
Elapsed simulation time (days) when the solution went unstable,
or ``None`` if the run finished successfully.
timestamp_unstable:
HEC-RAS date/time string for the instability, or ``None`` if the
run finished successfully (HEC-RAS writes ``"Not Applicable"``).
"""
solution: str
run_window: str
compute_time_total: str
compute_time_dss: str
max_cores: int
max_wsel_error: float | None
time_unstable: float | None
timestamp_unstable: str | None
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a dict with short, meaningful keys."""
return {
"solution": self.solution,
"run_window": self.run_window,
"compute_time_total": self.compute_time_total,
"compute_time_dss": self.compute_time_dss,
"max_cores": self.max_cores,
"max_wsel_error": self.max_wsel_error,
"time_unstable": self.time_unstable,
"timestamp_unstable": self.timestamp_unstable,
}
[docs]
def parse_run_window(
self,
) -> tuple[dt.datetime, dt.datetime] | None:
"""Parse :attr:`run_window` into ``(start, end)`` datetimes.
Splits the raw string on ``" to "`` and parses each half with
the HEC-RAS timestamp format ``"%d%b%Y %H:%M:%S"``.
Returns
-------
tuple[datetime, datetime] or None
``(start, end)`` as timezone-naive :class:`datetime.datetime`
objects, or ``None`` when the string is missing, malformed,
or cannot be parsed.
Examples
--------
::
s = hdf.compute_summary()
window = s.run.parse_run_window()
if window:
start, end = window
print(end - start) # wall-clock duration
"""
try:
left, right = self.run_window.split(" to ", maxsplit=1)
return (
parse_hec_datetime(left.strip(), fmt=_RAS_TS_FMT),
parse_hec_datetime(right.strip(), fmt=_RAS_TS_FMT),
)
except (ValueError, AttributeError):
return None
[docs]
@dataclasses.dataclass
class VolumeAccounting:
"""Overall volume accounting from ``Results/Unsteady/Summary/Volume Accounting``.
Covers the full model (1D + 2D combined).
Attributes
----------
units:
Volume units string written by HEC-RAS, e.g. ``"Acre Feet"`` or
``"1000 m^3"``.
vol_start:
Total storage volume at the start of the simulation.
vol_end:
Total storage volume at the end of the simulation.
inflow:
Total boundary flux of water into the model.
outflow:
Total boundary flux of water out of the model.
error:
Volume balance error (``vol_start + inflow - outflow - vol_end``).
error_pct:
Volume balance error as a percentage of total inflow.
"""
units: str
vol_start: float
vol_end: float
inflow: float
outflow: float
error: float
error_pct: float
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a dict with short, meaningful keys."""
return {
"units": self.units,
"vol_start": self.vol_start,
"vol_end": self.vol_end,
"inflow": self.inflow,
"outflow": self.outflow,
"error": self.error,
"error_pct": self.error_pct,
}
[docs]
@dataclasses.dataclass
class VolumeAccounting1D:
"""1-D component volume accounting from ``Volume Accounting 1D``.
Attributes
----------
units:
Volume units string, e.g. ``"Acre Feet"`` or ``"1000 m^3"``.
reach_vol_start:
Total 1-D reach storage at simulation start.
reach_vol_end:
Total 1-D reach storage at simulation end.
sa_vol_start:
Storage-area volume at simulation start.
sa_vol_end:
Storage-area volume at simulation end.
flow_us_in:
Cumulative upstream inflow across all reaches.
flow_ds_out:
Cumulative downstream outflow across all reaches.
hydro_lat:
Cumulative lateral hydrograph exchange (positive = into model).
hydro_sa:
Cumulative storage-area hydrograph exchange.
diversions:
Cumulative diversions (negative = water removed).
groundwater:
Cumulative groundwater exchange.
precip_excess:
Cumulative precipitation excess (rainfall-runoff applied to reaches).
"""
units: str
reach_vol_start: float
reach_vol_end: float
sa_vol_start: float
sa_vol_end: float
flow_us_in: float
flow_ds_out: float
hydro_lat: float
hydro_sa: float
diversions: float
groundwater: float
precip_excess: float
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a dict with short, meaningful keys."""
return {
"units": self.units,
"reach_vol_start": self.reach_vol_start,
"reach_vol_end": self.reach_vol_end,
"sa_vol_start": self.sa_vol_start,
"sa_vol_end": self.sa_vol_end,
"flow_us_in": self.flow_us_in,
"flow_ds_out": self.flow_ds_out,
"hydro_lat": self.hydro_lat,
"hydro_sa": self.hydro_sa,
"diversions": self.diversions,
"groundwater": self.groundwater,
"precip_excess": self.precip_excess,
}
[docs]
@dataclasses.dataclass
class VolumeAccounting2DArea:
"""Volume accounting for one named 2-D flow area.
One instance per child group under
``Results/Unsteady/Summary/Volume Accounting/Volume Accounting 2D``.
Attributes
----------
units:
Volume units string, e.g. ``"Acre Feet"`` or ``"1000 m^3"``.
vol_start:
2-D flow area storage at simulation start.
vol_end:
2-D flow area storage at simulation end.
cum_inflow:
Cumulative inflow into this 2-D area.
cum_outflow:
Cumulative outflow out of this 2-D area.
error:
Volume balance error for this area.
error_pct:
Volume balance error as a percentage of cumulative inflow.
"""
units: str
vol_start: float
vol_end: float
cum_inflow: float
cum_outflow: float
error: float
error_pct: float
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a dict with short, meaningful keys."""
return {
"units": self.units,
"vol_start": self.vol_start,
"vol_end": self.vol_end,
"cum_inflow": self.cum_inflow,
"cum_outflow": self.cum_outflow,
"error": self.error,
"error_pct": self.error_pct,
}
[docs]
@dataclasses.dataclass
class ComputeSummary:
"""Full unsteady simulation summary for a plan HDF file.
Aggregates all summary groups under ``Results/Unsteady/Summary``.
Attributes
----------
run:
Overall run metadata (solution status, timing, stability).
volume:
Model-wide volume accounting (1D + 2D combined).
volume_1d:
Volume accounting for 1-D components only.
volume_2d:
Per-area volume accounting for each 2-D flow area. Empty dict
when the model has no 2-D flow areas.
"""
run: RunStatus
volume: VolumeAccounting
volume_1d: VolumeAccounting1D
volume_2d: dict[str, VolumeAccounting2DArea]
[docs]
def to_dict(self) -> dict[str, Any]:
"""Return a nested dict with short, meaningful keys.
Top-level keys: ``"run"``, ``"volume"``, ``"volume_1d"``,
``"volume_2d"``. The ``"volume_2d"`` value is a dict keyed by
2-D flow area name.
"""
return {
"run": self.run.to_dict(),
"volume": self.volume.to_dict(),
"volume_1d": self.volume_1d.to_dict(),
"volume_2d": {
name: area.to_dict() for name, area in self.volume_2d.items()
},
}
[docs]
def errors(self) -> dict[str, float]:
"""Return key convergence and volume-balance error metrics.
Derives four scalar error measures from this summary that are useful
for quick quality-control checks after a run.
Returns
-------
dict[str, float]
A flat dictionary with the following keys:
``"wse"``
Maximum water-surface elevation error across all cells
(model units -- feet or metres). ``None`` when the simulation
went unstable before the run converged.
``"volume_error_pct"``
Overall volume balance error as a percentage of total inflow,
read directly from the HEC-RAS ``Results/Unsteady/Summary``
volume accounting block.
``"volume_1d_error_pct"``
1-D volume balance error (%) computed from the individual
1-D flux terms::
error = -(flow_us_in + net_other_fluxes + storage_change)
supply = reach_vol_start + sa_vol_start + flow_us_in
+ hydro_lat + hydro_sa + diversions
+ groundwater + precip_excess
error_pct = error / supply * 100
where *net_other_fluxes* = ``-flow_ds_out + hydro_lat +
hydro_sa + diversions + groundwater + precip_excess``.
``supply`` is the total 1-D water budget (initial storage plus
all non-DS-outflow fluxes); for a 1D-only model this matches
the HEC-RAS overall ``"Error Percent"`` attribute.
Returns ``nan`` when total 1-D supply is zero.
``"volume_2d_error_pct"``
Largest-magnitude 2-D volume balance error (%) across all
2-D flow areas, preserving the sign of the worst-case area.
Zero when the model has no 2-D flow areas.
See Also
--------
ok : Quality-control pass/fail check with threshold logging.
"""
v1d = self.volume_1d
outflow = (
v1d.flow_ds_out - v1d.hydro_lat - v1d.hydro_sa
- v1d.diversions - v1d.groundwater - v1d.precip_excess
)
storage_change = (
(v1d.reach_vol_start + v1d.sa_vol_start)
- (v1d.reach_vol_end + v1d.sa_vol_end)
)
error_1d = outflow - v1d.flow_us_in - storage_change
supply_1d = (
v1d.reach_vol_start + v1d.sa_vol_start
+ v1d.flow_us_in
+ v1d.hydro_lat
+ v1d.hydro_sa
+ v1d.diversions
+ v1d.groundwater
+ v1d.precip_excess
)
volume_1d_error_pct = (
error_1d / supply_1d * 100 if supply_1d != 0 else float("nan")
)
if (self.volume_2d
and not np.isclose(error_1d, 0.0)
and np.isclose(v1d.flow_ds_out, 0.0)):
logger.warning(
"Model has 2-D flow areas but zero 1-D downstream outflow; "
"1-D error percentages may be unreliable. Setting 1-D error to zero."
)
volume_1d_error_pct = 0.0
max_2d_error_pct = 0.0
for area in self.volume_2d.values():
if abs(area.error_pct) > abs(max_2d_error_pct):
max_2d_error_pct = area.error_pct
return {
"wse": self.run.max_wsel_error,
"volume_error_pct": self.volume.error_pct,
"volume_1d_error_pct": volume_1d_error_pct,
"volume_2d_error_pct": max_2d_error_pct,
}
[docs]
def ok(
self,
wse_threshold: float = 0.5,
volume_error_pct_threshold: float = 1.0,
volume_1d_error_pct_threshold: float = 1.0,
volume_2d_error_pct_threshold: float = 1.0,
) -> bool:
"""Return ``True`` if the simulation completed stably and all error
metrics are within their thresholds.
Checks stability first, then each error metric against its threshold.
All failures are collected and logged before returning so the caller
sees the complete picture in one pass.
If the simulation went unstable every failure is logged at
``CRITICAL``; otherwise threshold violations are logged at
``WARNING``.
Parameters
----------
wse_threshold:
Maximum allowable water-surface elevation error in model units
(feet or metres). Default: ``0.5``.
volume_error_pct_threshold:
Maximum allowable overall volume balance error (%).
Default: ``1.0``.
volume_1d_error_pct_threshold:
Maximum allowable 1-D volume balance error (%).
Default: ``1.0``.
volume_2d_error_pct_threshold:
Maximum allowable 2-D volume balance error (%) per 2-D flow area.
Default: ``1.0``.
Returns
-------
bool
``True`` if the simulation is stable and all metrics are within
bounds; ``False`` otherwise.
See Also
--------
errors : Raw error metric values.
"""
failures: list[str] = []
unstable = self.run.time_unstable is not None
if unstable:
failures.append(
f"simulation went unstable at t={self.run.time_unstable}"
)
errs = self.errors()
wse = errs["wse"]
if wse is None or math.isnan(wse) or abs(wse) > wse_threshold:
failures.append(
f"WSE error {wse!r} exceeds threshold {wse_threshold}"
)
vol_err = errs["volume_error_pct"]
if math.isnan(vol_err) or abs(vol_err) > volume_error_pct_threshold:
failures.append(
f"overall volume error {vol_err:.4f}% exceeds threshold "
f"{volume_error_pct_threshold}%"
)
vol_1d_err = errs["volume_1d_error_pct"]
if math.isnan(vol_1d_err) or abs(vol_1d_err) > volume_1d_error_pct_threshold:
vol_1d_str = "nan" if math.isnan(vol_1d_err) else f"{vol_1d_err:.4f}%"
failures.append(
f"1D volume error {vol_1d_str} exceeds threshold "
f"{volume_1d_error_pct_threshold}%"
)
for area_name, area in self.volume_2d.items():
exceeds = (
math.isnan(area.error_pct)
or abs(area.error_pct) > volume_2d_error_pct_threshold
)
if exceeds:
failures.append(
f"2D area '{area_name}' volume error {area.error_pct:.4f}% "
f"exceeds threshold {volume_2d_error_pct_threshold}%"
)
if not failures:
return True
log = logger.critical if unstable else logger.warning
for msg in failures:
log("compute_ok: %s", msg)
return False
[docs]
class UnsteadyPlan(_PlanHdf, Geometry):
"""Read HEC-RAS plan HDF5 output files (``*.p*.hdf``).
A plan HDF file contains the same ``Geometry/`` data as a geometry HDF
file, *plus* ``Results/Unsteady/...`` time-series and summary output.
Parameters
----------
filename:
Path to the plan HDF file. The ``.hdf`` suffix is appended
automatically if absent.
Examples
--------
::
with UnsteadyPlan("MyModel.p01") as hdf:
ts = hdf.mapping_timestamps
area = hdf.flow_areas["spillway"]
wse = area.water_surface[10] # one timestep
depth = area.get_depth(timestep=10)
speed = area.get_cell_velocity(10, component="speed")
max_d = area.get_max_depth()
# requires rasterio + scipy:
area.export_raster("depth", "depth.tif", timestep=None,
cell_size=5.0, crs="EPSG:26910")
# or get an in-memory dataset:
ds = area.export_raster("depth", timestep=None, cell_size=5.0)
arr = ds.read(1)
ds.close()
"""
def __init__(self, filename: str | Path, program_directory: str | Path | None = None) -> None:
super().__init__(filename)
self._program_directory = Path(program_directory) if program_directory else None
self._geom_view: Geometry | None = None
self._plan_flow_areas: FlowAreaResultsCollection | None = None
self._plan_storage_areas_cache: dict[str, StorageAreaResultsCollection] = {}
self._plan_structures_cache: dict[str, StructureResultsCollection] = {}
self._plan_cross_sections_cache: dict[str, CrossSectionResultsCollection] = {}
# ------------------------------------------------------------------
# Runtime log
# ------------------------------------------------------------------
[docs]
def runtime_log(self) -> UnsteadyRuntimeLog:
"""Read the runtime compute log from ``Results/Summary/``.
Returns
-------
UnsteadyRuntimeLog
Log container with the full text/RTF compute messages, the
compute-process table, and unsteady-specific parsing methods
such as :meth:`~UnsteadyRuntimeLog.max_iterations` and
:meth:`~UnsteadyRuntimeLog.adaptive_timesteps`.
Raises
------
KeyError
If ``Results/Summary`` is absent from the HDF file.
"""
return UnsteadyRuntimeLog(*self._runtime_log_raw())
# ------------------------------------------------------------------
# File metadata
# ------------------------------------------------------------------
[docs]
def compute_summary(self) -> ComputeSummary:
"""Return the full unsteady simulation summary.
Reads all groups under ``Results/Unsteady/Summary`` and returns a
:class:`ComputeSummary` dataclass aggregating run metadata,
overall volume accounting, 1-D volume accounting, and per-area
2-D volume accounting.
Call :meth:`ComputeSummary.to_dict` on the result for a
nested dict with short, meaningful keys.
Returns
-------
ComputeSummary
Raises
------
KeyError
If ``Results/Unsteady/Summary`` is absent -- e.g. the file is
a steady-flow plan or the simulation has not been run yet.
Examples
--------
::
with UnsteadyPlan("MyModel.p01") as hdf:
s = hdf.compute_summary()
print(s.run.solution)
print(s.volume.error_pct)
d = s.to_dict()
"""
grp = self._hdf.get(_RUN_SUM)
if grp is None:
raise KeyError(
f"'{_RUN_SUM}' not found. "
"Ensure this is an unsteady-flow plan HDF file that has been run."
)
def _str(attrs: Any, key: str) -> str:
v = attrs[key]
return v.decode() if isinstance(v, (bytes, np.bytes_)) else str(v)
def _float_or_none(attrs: Any, key: str) -> float | None:
v = float(attrs[key])
return None if np.isnan(v) else v
# -- RunStatus ---------------------------------------------------
a = grp.attrs
unstable_ts = _str(a, "Time Stamp Solution Went Unstable")
run = RunStatus(
solution=_str(a, "Solution"),
run_window=_str(a, "Run Time Window"),
compute_time_total=_str(a, "Computation Time Total"),
compute_time_dss=_str(a, "Computation Time DSS"),
max_cores=int(a["Maximum number of cores"]),
max_wsel_error=_float_or_none(a, "Maximum WSEL Error"),
time_unstable=_float_or_none(a, "Time Solution Went Unstable"),
timestamp_unstable=(
None if unstable_ts == "Not Applicable" else unstable_ts
),
)
# -- VolumeAccounting (overall) -----------------------------------
va_grp = grp["Volume Accounting"]
a = va_grp.attrs
volume = VolumeAccounting(
units=_str(a, "Vol Accounting in"),
vol_start=float(a["Volume Starting"]),
vol_end=float(a["Volume Ending"]),
inflow=float(a["Total Boundary Flux of Water In"]),
outflow=float(a["Total Boundary Flux of Water Out"]),
error=float(a["Error"]),
error_pct=float(a["Error Percent"]),
)
# -- VolumeAccounting1D ------------------------------------------
a = va_grp["Volume Accounting 1D"].attrs
precip_key = next(k for k in a if k.startswith("Precip Excess"))
volume_1d = VolumeAccounting1D(
units=_str(a, "Vol Accounting in"),
reach_vol_start=float(a["Reach Start 1D"]),
reach_vol_end=float(a["Reach Final 1D"]),
sa_vol_start=float(a["SA Starting"]),
sa_vol_end=float(a["SA Final"]),
flow_us_in=float(a["Flow US In"]),
flow_ds_out=float(a["Flow DS Out"]),
hydro_lat=float(a["Hydro Lat"]),
hydro_sa=float(a["Hydro SA"]),
diversions=float(a["Diversions"]),
groundwater=float(a["Groundwater"]),
precip_excess=float(a[precip_key]),
)
# -- VolumeAccounting2D (optional) --------------------------------
volume_2d: dict[str, VolumeAccounting2DArea] = {}
vd_grp = va_grp.get("Volume Accounting 2D")
if vd_grp is not None:
for area_name, area_grp in vd_grp.items():
a = area_grp.attrs
volume_2d[area_name] = VolumeAccounting2DArea(
units=_str(a, "Vol Accounting in"),
vol_start=float(a["Vol Starting"]),
vol_end=float(a["Vol Ending"]),
cum_inflow=float(a["Cum Inflow"]),
cum_outflow=float(a["Cum Outflow"]),
error=float(a["Error"]),
error_pct=float(a["Error Percent"]),
)
return ComputeSummary(
run=run,
volume=volume,
volume_1d=volume_1d,
volume_2d=volume_2d,
)
[docs]
def compute_errors(self) -> dict[str, float]:
"""Return key convergence and volume-balance error metrics.
Convenience wrapper: calls :meth:`compute_summary` and delegates to
:meth:`ComputeSummary.errors`. See that method for full documentation.
See Also
--------
ComputeSummary.errors : Full documentation and return-value details.
compute_ok : Quality-control pass/fail check with threshold logging.
"""
return self.compute_summary().errors()
[docs]
def compute_ok(
self,
wse_threshold: float = 0.5,
volume_error_pct_threshold: float = 1.0,
volume_1d_error_pct_threshold: float = 1.0,
volume_2d_error_pct_threshold: float = 1.0,
) -> bool:
"""Return ``True`` if the simulation completed stably and all error
metrics are within their thresholds.
Convenience wrapper: calls :meth:`compute_summary` and delegates to
:meth:`ComputeSummary.ok`. See that method for full documentation,
including parameter descriptions and logging behaviour.
Examples
--------
Quick pass/fail check after a run::
if not hdf.compute_ok():
raise RuntimeError("Simulation quality check failed.")
Custom thresholds::
if not hdf.compute_ok(wse_threshold=0.1, volume_error_pct_threshold=0.5):
logger.warning("Tight-tolerance check failed.")
See Also
--------
ComputeSummary.ok : Full documentation and parameter descriptions.
compute_errors : Raw error metric values.
"""
return self.compute_summary().ok(
wse_threshold=wse_threshold,
volume_error_pct_threshold=volume_error_pct_threshold,
volume_1d_error_pct_threshold=volume_1d_error_pct_threshold,
volume_2d_error_pct_threshold=volume_2d_error_pct_threshold,
)
@property
def ras_version(self) -> str:
"""HEC-RAS version string from the plan HDF root attribute.
Returns the ``File Version`` root attribute, e.g.
``'HEC-RAS 6.6 September 2024'``.
"""
raw = self._hdf.attrs["File Version"]
return raw.decode() if isinstance(raw, (bytes, np.bytes_)) else str(raw)
@property
def computation_interval(self) -> dt.timedelta | None:
"""Base computation time step, or ``None`` if absent.
Read from ``Plan Data/Plan Information`` attribute
``Computation Time Step Base``. HEC-RAS stores the value as a
concatenated number-unit string, e.g. ``'20SEC'``, ``'5MIN'``,
``'1HR'``, ``'1HOUR'``, ``'1DAY'``.
"""
raw = self._plan_info_attr("Computation Time Step Base")
return None if raw is None else parse_interval(raw)
@property
def mapping_interval(self) -> dt.timedelta | None:
"""Mapping output interval, or ``None`` if absent.
Read from ``Plan Data/Plan Information`` attribute
``Base Output Interval``. HEC-RAS stores the value as a
concatenated number-unit string, e.g. ``'5MIN'``.
"""
raw = self._plan_info_attr("Base Output Interval")
return None if raw is None else parse_interval(raw)
# ------------------------------------------------------------------
# Time stamps
# ------------------------------------------------------------------
@property
def mapping_timestamps(self) -> pd.DatetimeIndex:
"""Simulation output time stamps as a ``pd.DatetimeIndex``.
Parsed from the ``Time Date Stamp`` dataset written by HEC-RAS.
Format: ``DD Mon YYYY HH:MM:SS`` (e.g. ``03Jan2000 00:00:00``).
"""
ds = self._hdf.get(_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp dataset not found at '{_TIME_STAMP_DS}'. "
"Ensure this is an unsteady-flow plan HDF file."
)
raw = np.array(ds).astype(str)
return _parse_hec_ts_array(raw, _RAS_TS_FMT)
@property
def output_timestamps(self) -> pd.DatetimeIndex:
"""DSS hydrograph output time stamps as a ``pd.DatetimeIndex``.
Parsed from ``Results/.../DSS Hydrograph Output/Unsteady Time
Series/Time Date Stamp``. Format: ``DD Mon YYYY HH:MM:SS``.
"""
ds = self._hdf.get(_DSS_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp dataset not found at '{_DSS_TIME_STAMP_DS}'. "
"Ensure DSS hydrograph output was written for this plan."
)
raw = np.array(ds).astype(str)
return _parse_hec_ts_array(raw, _RAS_TS_FMT)
@property
def output_interval(self) -> dt.timedelta | None:
"""DSS hydrograph output interval, or ``None`` if absent.
Derived from the difference between the first two hydrograph
time stamps.
"""
ds = self._hdf.get(_DSS_TIME_STAMP_DS)
if ds is None or len(ds) < 2:
return None
raw = np.array(ds[:2]).astype(str)
ts = _parse_hec_ts_array(raw, _RAS_TS_FMT)
return ts[1] - ts[0]
@property
def detailed_interval(self) -> dt.timedelta | None:
"""Detailed output interval, or ``None`` if absent.
Derived from the difference between the first two timestamps in the
DSS Profile Output ``Time Date Stamp`` dataset. Used by both
``"profile"`` and ``"post_process"`` output blocks (same HEC-RAS
Detailed Output Interval setting).
"""
ds = self._hdf.get(_DSS_PROF_TIME_STAMP_DS)
if ds is None or len(ds) < 2:
return None
raw = np.array(ds[:2]).astype(str)
ts = _parse_hec_ts_array(raw, _RAS_TS_FMT)
return ts[1] - ts[0]
@property
def detailed_timestamps(self) -> pd.DatetimeIndex:
"""Detailed output interval timestamps as a ``pd.DatetimeIndex``.
Reads from the DSS Profile Output ``Time Date Stamp`` dataset.
Both ``cross_sections("profile")`` / ``storage_areas("profile")``
and ``cross_sections("post_process")`` / ``storage_areas("post_process")``
share this interval (same HEC-RAS Detailed Output Interval setting).
Raises
------
KeyError
If DSS Profile Output was not written for this plan.
"""
ds = self._hdf.get(_DSS_PROF_TIME_STAMP_DS)
if ds is None:
raise KeyError(
f"Time Date Stamp not found at '{_DSS_PROF_TIME_STAMP_DS}'. "
"Ensure DSS Profile Output was written for this plan."
)
raw = np.array(ds).astype(str)
return _parse_hec_ts_array(raw, _RAS_TS_FMT)
@property
def n_mapping_timestamps(self) -> int | None:
"""Number of mapping output time steps, or ``None`` for steady-flow plans.
Reads only the dataset shape -- no timestamp data is loaded.
"""
ds = self._hdf.get(_TIME_STAMP_DS)
return None if ds is None else len(ds)
@property
def n_output_timestamps(self) -> int | None:
"""Number of DSS hydrograph output time steps, or ``None`` if absent.
Reads only the dataset shape -- no timestamp data is loaded.
"""
ds = self._hdf.get(_DSS_TIME_STAMP_DS)
return None if ds is None else len(ds)
@property
def n_detailed_timestamps(self) -> int | None:
"""Number of detailed output time steps, or ``None`` if absent.
Reads only the shape of the DSS Profile ``Time Date Stamp`` dataset.
"""
ds = self._hdf.get(_DSS_PROF_TIME_STAMP_DS)
return None if ds is None else len(ds)
# ------------------------------------------------------------------
# Collections (override Geometry equivalents with results-aware types)
# ------------------------------------------------------------------
@property
def flow_areas(self) -> FlowAreaResultsCollection:
"""Access 2-D flow areas with both geometry and results data."""
if self._plan_flow_areas is None:
self._plan_flow_areas = FlowAreaResultsCollection(self._hdf)
return self._plan_flow_areas
[docs]
def storage_areas(
self,
output: Literal["mapping", "output", "profile", "post_process"] = "mapping",
) -> StorageAreaResultsCollection:
"""Access storage areas with geometry and plan results.
Parameters
----------
output : {"mapping", "output", "profile", "post_process"}, optional
Which output block to read. Defaults to ``"mapping"``.
``"mapping"`` — Base Output (includes summary max/min WSE).
``"output"`` — DSS Hydrograph Output.
``"profile"`` — DSS Profile Output; timestamps = ``detailed_timestamps``.
``"post_process"`` — Post Process Profiles; timestamps =
``detailed_timestamps``; row 0 (Max WS envelope) is excluded from
all time-series arrays. Summary ``max_wse`` / ``min_wse`` raise
``KeyError`` (summary group absent from this block).
Returns
-------
StorageAreaResultsCollection
Exposes a ``timestamps`` property with the block's ``pd.DatetimeIndex``,
shared as the index of every per-SA ``pd.Series`` property.
"""
if output not in self._plan_storage_areas_cache:
self._plan_storage_areas_cache[output] = StorageAreaResultsCollection(
self._hdf, output=output
)
return self._plan_storage_areas_cache[output]
[docs]
def structures(
self,
output: Literal["mapping", "output", "profile", "post_process"] = "output",
) -> StructureResultsCollection:
"""Access all structures with geometry *and* plan results.
Returns a :class:`StructureResultsCollection` where each item is
upgraded to the matching results class when plan output is present.
Parameters
----------
output : {"mapping", "output", "profile", "post_process"}, optional
Which output block to read. Defaults to ``"output"``.
``"mapping"`` — Base Output (all four structure types).
``"output"`` — DSS Hydrograph Output (all four structure types).
``"profile"`` — DSS Profile Output (all four structure types).
``"post_process"`` — Post Process Profiles; Lateral and SA/2D
connections only. Inline and Bridge items fall back to plain
geometry objects.
Returns
-------
StructureResultsCollection
Exposes a ``timestamps`` property with the block's ``pd.DatetimeIndex``,
shared as the index of every per-structure ``pd.Series`` property.
Use :attr:`~rivia.hdf.StructureCollection.connections`,
:attr:`~rivia.hdf.StructureCollection.inlines`,
:attr:`~rivia.hdf.StructureCollection.laterals`, and
:attr:`~rivia.hdf.StructureCollection.bridges` for filtered access.
"""
if output not in self._plan_structures_cache:
self._plan_structures_cache[output] = StructureResultsCollection(
self._hdf, output=output
)
return self._plan_structures_cache[output]
[docs]
def cross_sections(
self,
output: Literal["mapping", "output", "profile", "post_process"] = "mapping",
) -> CrossSectionResultsCollection:
"""1-D cross sections with geometry and time-series results.
Parameters
----------
output : {"mapping", "output", "profile", "post_process"}, optional
Which output block to read.
``"mapping"`` (default) — Base Output at the mapping interval.
Returns :class:`CrossSectionResultsCollection` whose items are
:class:`CrossSectionMappingResults`. Per-XS ``pd.Series``
properties: ``wse``, ``flow``, ``flow_lateral``,
``velocity_channel``, ``velocity_total``, ``flow_cumulative``.
``"output"`` — DSS Hydrograph Output at the hydrograph interval.
Returns :class:`CrossSectionResultsCollection` whose items are
:class:`CrossSectionOutputResults`. Per-XS ``pd.Series``
properties: ``wse``, ``flow``, ``flow_cumulative``.
``"profile"`` — DSS Profile Output at the detailed interval.
Returns :class:`CrossSectionResultsCollection` whose items are
:class:`CrossSectionProfileResults`. Per-XS ``pd.Series``
properties: ``wse``, ``flow``, ``cross_section_attributes``.
``"post_process"`` — Post Process Profiles. Returns a
:class:`CrossSectionPostProcessResultsCollection` with collection-level
``profile_table(variable)`` method and named properties (``wse``,
``flow``, ``velocity_channel``, …) each returning a
location × profiles ``pd.DataFrame``. ``[key]`` returns a
:class:`CrossSectionPostProcessResults` whose named properties
(``wse``, ``flow``, ``velocity_channel``, …) each return a
``pd.Series`` indexed by ``pd.DatetimeIndex`` (timeseries only,
Max WS envelope excluded).
Returns
-------
CrossSectionResultsCollection
Collection supporting ``[key]``, integer index, and ``names``.
Timestamps are available as ``coll.timestamps``.
Raises
------
ValueError
If *output* is not one of the four recognised values.
"""
if output not in self._plan_cross_sections_cache:
if output == "mapping":
coll: CrossSectionResultsCollection = CrossSectionResultsCollection(
self._hdf, _TS_XS,
result_cls=CrossSectionMappingResults,
timestamps_fn=lambda: self.mapping_timestamps,
)
elif output == "output":
coll = CrossSectionResultsCollection(
self._hdf, _DSS_XS,
result_cls=CrossSectionOutputResults,
timestamps_fn=lambda: self.output_timestamps,
)
elif output == "profile":
coll = CrossSectionResultsCollection(
self._hdf, _DSS_PROF_XS,
result_cls=CrossSectionProfileResults,
timestamps_fn=lambda: self.detailed_timestamps,
)
elif output == "post_process":
coll = CrossSectionPostProcessResultsCollection(
self._hdf, _POSTPROC_XS,
result_cls=CrossSectionPostProcessResults,
attrs_path=_POSTPROC_GEOM_ATTRS,
timestamps_fn=lambda: self.detailed_timestamps,
)
else:
raise ValueError(
f"output={output!r} is not valid; "
"choose 'mapping', 'output', 'profile', or 'post_process'."
)
self._plan_cross_sections_cache[output] = coll
return self._plan_cross_sections_cache[output]
[docs]
def sa2d_connections(
self,
output: Literal["mapping", "output", "profile", "post_process"] = "output",
) -> dict[str, SA2DConnection]:
"""SA/2D hydraulic connections keyed by geometry name.
Convenience delegate to ``structures(output).connections``.
Parameters
----------
output : {"mapping", "output", "profile", "post_process"}, optional
Passed through to :meth:`structures`. Defaults to ``"output"``.
Returns
-------
dict[str, SA2DConnection]
Items are :class:`SA2DConnectionResults` when plan output is
present, plain :class:`SA2DConnection` otherwise.
"""
return self.structures(output).connections