Merge pull request #55 from fema-ffrd/Feature/ras_1d

Feature/ras 1d

docs/source/RasGeomHdf.rst

@@ -15,4 +15,8 @@ RasGeomHdf
              structures,
              get_geom_attrs,
              get_geom_structures_attrs,
-             get_geom_2d_flow_area_attrs
+             get_geom_2d_flow_area_attrs,
+             cross_sections_elevations,
+             cross_sections
+             river_reaches
+

docs/source/RasPlanHdf.rst

@@ -18,4 +18,13 @@ RasPlanHdf
              get_meteorology_precip_attrs,
              get_results_unsteady_attrs,
              get_results_unsteady_summary_attrs,
-             get_results_volume_accounting_attrs
+             get_results_volume_accounting_attrs,
+             cross_sections_additional_velocity_total,
+             cross_sections_additional_area_total,
+             cross_sections_additional_enc_station_right,
+             cross_sections_additional_enc_station_left,
+             cross_sections_energy_grade,
+             cross_sections_flow,
+             cross_sections_wsel,
+             steady_flow_names,
+             steady_profile_xs_output

pyproject.toml

@@ -12,8 +12,8 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
 ]
-version = "0.3.0"
-dependencies = ["h5py", "geopandas", "pyarrow", "xarray"]
+version = "0.4.0"
+dependencies = ["h5py", "geopandas>=0.14,<0.15", "pyarrow", "xarray"]
 
 [project.optional-dependencies]
 dev = ["pre-commit", "ruff", "pytest", "pytest-cov"]

src/rashdf/geom.py

@@ -1,27 +1,27 @@
 """HEC-RAS Geometry HDF class."""
 
-from .base import RasHdf
-from .utils import (
-    convert_ras_hdf_string,
-    get_first_hdf_group,
-    hdf5_attrs_to_dict,
-    convert_ras_hdf_value,
-)
+from typing import Dict, List, Optional
 
 import numpy as np
 import pandas as pd
 from geopandas import GeoDataFrame
 from pyproj import CRS
 from shapely import (
-    Polygon,
-    Point,
     LineString,
     MultiLineString,
     MultiPolygon,
+    Point,
+    Polygon,
     polygonize_full,
 )
 
-from typing import Dict, List, Optional
+from .base import RasHdf
+from .utils import (
+    convert_ras_hdf_string,
+    convert_ras_hdf_value,
+    get_first_hdf_group,
+    hdf5_attrs_to_dict,
+)
 
 
 class RasGeomHdf(RasHdf):
@@ -465,11 +465,97 @@ def classification_polygons(self) -> GeoDataFrame:  # noqa D102
     def terrain_modifications(self) -> GeoDataFrame:  # noqa D102
         raise NotImplementedError
 
-    def cross_sections(self) -> GeoDataFrame:  # noqa D102
-        raise NotImplementedError
+    def cross_sections(self, datetime_to_str: bool = False) -> GeoDataFrame:
+        """Return the model 1D cross sections.
 
-    def river_reaches(self) -> GeoDataFrame:  # noqa D102
-        raise NotImplementedError
+        Returns
+        -------
+        GeoDataFrame
+            A GeoDataFrame containing the model 1D cross sections if they exist.
+        """
+        if "/Geometry/Cross Sections" not in self:
+            return GeoDataFrame()
+
+        xs_data = self["/Geometry/Cross Sections"]
+        v_conv_val = np.vectorize(convert_ras_hdf_value)
+        xs_attrs = xs_data["Attributes"][()]
+        xs_dict = {"xs_id": range(xs_attrs.shape[0])}
+        xs_dict.update(
+            {name: v_conv_val(xs_attrs[name]) for name in xs_attrs.dtype.names}
+        )
+        geoms = list()
+        for pnt_start, pnt_cnt, part_start, part_cnt in xs_data["Polyline Info"][()]:
+            points = xs_data["Polyline Points"][()][pnt_start : pnt_start + pnt_cnt]
+            if part_cnt == 1:
+                geoms.append(LineString(points))
+            else:
+                parts = xs_data["Polyline Parts"][()][
+                    part_start : part_start + part_cnt
+                ]
+                geoms.append(
+                    MultiLineString(
+                        list(
+                            points[part_pnt_start : part_pnt_start + part_pnt_cnt]
+                            for part_pnt_start, part_pnt_cnt in parts
+                        )
+                    )
+                )
+        xs_gdf = GeoDataFrame(
+            xs_dict,
+            geometry=geoms,
+            crs=self.projection(),
+        )
+        if datetime_to_str:
+            xs_gdf["Last Edited"] = xs_gdf["Last Edited"].apply(
+                lambda x: pd.Timestamp.isoformat(x)
+            )
+        return xs_gdf
+
+    def river_reaches(self, datetime_to_str: bool = False) -> GeoDataFrame:
+        """Return the model 1D river reach lines.
+
+        Returns
+        -------
+        GeoDataFrame
+            A GeoDataFrame containing the model 1D river reach lines if they exist.
+        """
+        if "/Geometry/River Centerlines" not in self:
+            return GeoDataFrame()
+
+        river_data = self["/Geometry/River Centerlines"]
+        v_conv_val = np.vectorize(convert_ras_hdf_value)
+        river_attrs = river_data["Attributes"][()]
+        river_dict = {"river_id": range(river_attrs.shape[0])}
+        river_dict.update(
+            {name: v_conv_val(river_attrs[name]) for name in river_attrs.dtype.names}
+        )
+        geoms = list()
+        for pnt_start, pnt_cnt, part_start, part_cnt in river_data["Polyline Info"][()]:
+            points = river_data["Polyline Points"][()][pnt_start : pnt_start + pnt_cnt]
+            if part_cnt == 1:
+                geoms.append(LineString(points))
+            else:
+                parts = river_data["Polyline Parts"][()][
+                    part_start : part_start + part_cnt
+                ]
+                geoms.append(
+                    MultiLineString(
+                        list(
+                            points[part_pnt_start : part_pnt_start + part_pnt_cnt]
+                            for part_pnt_start, part_pnt_cnt in parts
+                        )
+                    )
+                )
+        river_gdf = GeoDataFrame(
+            river_dict,
+            geometry=geoms,
+            crs=self.projection(),
+        )
+        if datetime_to_str:
+            river_gdf["Last Edited"] = river_gdf["Last Edited"].apply(
+                lambda x: pd.Timestamp.isoformat(x)
+            )
+        return river_gdf
 
     def flowpaths(self) -> GeoDataFrame:  # noqa D102
         raise NotImplementedError
@@ -488,3 +574,33 @@ def ineffective_points(self) -> GeoDataFrame:  # noqa D102
 
     def blocked_obstructions(self) -> GeoDataFrame:  # noqa D102
         raise NotImplementedError
+
+    def cross_sections_elevations(self, round_to: int = 2) -> pd.DataFrame:
+        """Return the model cross section elevation information.
+
+        Returns
+        -------
+        DataFrame
+            A DataFrame containing the model cross section elevation information if they exist.
+        """
+        path = "/Geometry/Cross Sections"
+        if path not in self:
+            return pd.DataFrame()
+
+        xselev_data = self[path]
+        xs_df = self.cross_sections()
+        elevations = list()
+        for part_start, part_cnt in xselev_data["Station Elevation Info"][()]:
+            xzdata = xselev_data["Station Elevation Values"][()][
+                part_start : part_start + part_cnt
+            ]
+            elevations.append(xzdata)
+
+        xs_elev_df = xs_df[
+            ["xs_id", "River", "Reach", "RS", "Left Bank", "Right Bank"]
+        ].copy()
+        xs_elev_df["Left Bank"] = xs_elev_df["Left Bank"].round(round_to).astype(str)
+        xs_elev_df["Right Bank"] = xs_elev_df["Right Bank"].round(round_to).astype(str)
+        xs_elev_df["elevation info"] = elevations
+
+        return xs_elev_df

src/rashdf/plan.py

@@ -25,6 +25,26 @@ class RasPlanHdfError(Exception):
     pass
 
 
+class XsSteadyOutputVar(Enum):
+    """Summary of steady cross section output variables."""
+
+    ENERGY_GRADE = "Energy Grade"
+    FLOW = "Flow"
+    WATER_SURFACE = "Water Surface"
+    ENCROACHMENT_STATION_LEFT = "Encroachment Station Left"
+    ENCROACHMENT_STATION_RIGHT = "Encroachment Station Right"
+    AREA_INEFFECTIVE_TOTAL = "Area including Ineffective Total"
+    VELOCITY_TOTAL = "Velocity Total"
+
+
+XS_STEADY_OUTPUT_ADDITIONAL = [
+    XsSteadyOutputVar.ENCROACHMENT_STATION_LEFT,
+    XsSteadyOutputVar.ENCROACHMENT_STATION_RIGHT,
+    XsSteadyOutputVar.AREA_INEFFECTIVE_TOTAL,
+    XsSteadyOutputVar.VELOCITY_TOTAL,
+]
+
+
 class SummaryOutputVar(Enum):
     """Summary output variables."""
 
@@ -144,6 +164,12 @@ class RasPlanHdf(RasGeomHdf):
     )
     UNSTEADY_TIME_SERIES_PATH = f"{BASE_OUTPUT_PATH}/Unsteady Time Series"
 
+    RESULTS_STEADY_PATH = "Results/Steady"
+    BASE_STEADY_PATH = f"{RESULTS_STEADY_PATH}/Output/Output Blocks/Base Output"
+    STEADY_PROFILES_PATH = f"{BASE_STEADY_PATH}/Steady Profiles"
+    STEADY_XS_PATH = f"{STEADY_PROFILES_PATH}/Cross Sections"
+    STEADY_XS_ADDITIONAL_PATH = f"{STEADY_XS_PATH}/Additional Variables"
+
     def __init__(self, name: str, **kwargs):
         """Open a HEC-RAS Plan HDF file.
 
@@ -957,3 +983,127 @@ def get_results_volume_accounting_attrs(self) -> Dict:
 
     def enroachment_points(self) -> GeoDataFrame:  # noqa: D102
         raise NotImplementedError
+
+    def steady_flow_names(self) -> list:
+        """Return the profile information for each steady flow event.
+
+        Returns
+        -------
+        DataFrame
+            A Dataframe containing the profile names for each event
+        """
+        if self.STEADY_PROFILES_PATH not in self:
+            return pd.DataFrame()
+
+        profile_data = self[self.STEADY_PROFILES_PATH]
+        profile_attrs = profile_data["Profile Names"][()]
+
+        return [x.decode("utf-8") for x in profile_attrs]
+
+    def steady_profile_xs_output(
+        self, var: XsSteadyOutputVar, round_to: int = 2
+    ) -> DataFrame:
+        """Create a Dataframe from steady cross section results based on path.
+
+        Parameters
+        ----------
+        var : XsSteadyOutputVar
+            The name of the table in the steady results that information is to be retireved from.
+
+        round_to : int, optional
+            Number of decimal places to round output data to.
+
+        Returns
+        -------
+            Dataframe with desired hdf data.
+        """
+        if var in XS_STEADY_OUTPUT_ADDITIONAL:
+            path = f"{self.STEADY_XS_ADDITIONAL_PATH}/{var.value}"
+        else:
+            path = f"{self.STEADY_XS_PATH}/{var.value}"
+        if path not in self:
+            return DataFrame()
+
+        profiles = self.steady_flow_names()
+
+        steady_data = self[path]
+        df = DataFrame(steady_data, index=profiles)
+        df_t = df.T.copy()
+        for p in profiles:
+            df_t[p] = df_t[p].apply(lambda x: round(x, round_to))
+
+        return df_t
+
+    def cross_sections_wsel(self) -> DataFrame:
+        """Return the water surface elevation information for each 1D Cross Section.
+
+        Returns
+        -------
+        DataFrame
+            A Dataframe containing the water surface elevations for each cross section and event
+        """
+        return self.steady_profile_xs_output(XsSteadyOutputVar.WATER_SURFACE)
+
+    def cross_sections_flow(self) -> DataFrame:
+        """Return the Flow information for each 1D Cross Section.
+
+        Returns
+        -------
+        DataFrame
+            A Dataframe containing the flow for each cross section and event
+        """
+        return self.steady_profile_xs_output(XsSteadyOutputVar.FLOW)
+
+    def cross_sections_energy_grade(self) -> DataFrame:
+        """Return the energy grade information for each 1D Cross Section.
+
+        Returns
+        -------
+        DataFrame
+            A Dataframe containing the water surface elevations for each cross section and event
+        """
+        return self.steady_profile_xs_output(XsSteadyOutputVar.ENERGY_GRADE)
+
+    def cross_sections_additional_enc_station_left(self) -> DataFrame:
+        """Return the left side encroachment information for a floodway plan hdf.
+
+        Returns
+        -------
+        DataFrame
+            A DataFrame containing the cross sections left side encroachment stations
+        """
+        return self.steady_profile_xs_output(
+            XsSteadyOutputVar.ENCROACHMENT_STATION_LEFT
+        )
+
+    def cross_sections_additional_enc_station_right(self) -> DataFrame:
+        """Return the right side encroachment information for a floodway plan hdf.
+
+        Returns
+        -------
+        DataFrame
+            A DataFrame containing the cross sections right side encroachment stations
+        """
+        return self.steady_profile_xs_output(
+            XsSteadyOutputVar.ENCROACHMENT_STATION_RIGHT
+        )
+
+    def cross_sections_additional_area_total(self) -> DataFrame:
+        """Return the 1D cross section area for each profile.
+
+        Returns
+        -------
+        DataFrame
+            A DataFrame containing the wet area inside the cross sections
+        """
+        return self.steady_profile_xs_output(XsSteadyOutputVar.AREA_INEFFECTIVE_TOTAL)
+
+    def cross_sections_additional_velocity_total(self) -> DataFrame:
+        """Return the 1D cross section velocity for each profile.
+
+        Returns
+        -------
+        DataFrame
+            A DataFrame containing the velocity inside the cross sections
+        """
+        return self.steady_profile_xs_output(XsSteadyOutputVar.VELOCITY_TOTAL)

tests/__init__.py

@@ -1,6 +1,6 @@
 import h5py
 from geopandas import GeoDataFrame
-
+import json
 import hashlib
 from pathlib import Path
 
@@ -20,3 +20,8 @@ def _gdf_matches_json(gdf: GeoDataFrame, json_file: Path) -> bool:
 def get_sha1_hash(path: Path) -> str:
     with open(path, "rb") as f:
         return hashlib.sha1(f.read()).hexdigest()
+
+
+def _gdf_matches_json_alt(gdf: GeoDataFrame, json_file: Path) -> bool:
+    with open(json_file) as j:
+        return json.loads(gdf.to_json()) == json.load(j)