Several fixes to processing of cross-sections (#153)

* changes to setup_crosssections from Santa Cruz project * improvements to xyz and point x-sections from Santa Cruz * improvements to yz x-sections from Santa Cruz * 3 points are sufficient to define x-section * changes to xyz cross-sections workflow * finalizing xyz test * formatting using black * include maxdist * Update changelog.rst * New author in pyproject.toml * removed dummy tests * removed 'special' x-section type --------- Co-authored-by: veenstrajelmer <[email protected]>
Deltares · Oct 4, 2024 · 2d66900 · 2d66900
1 parent 8862432
commit 2d66900
Show file tree

Hide file tree

Showing 5 changed files with 96 additions and 21 deletions.
diff --git a/docs/changelog.rst b/docs/changelog.rst
@@ -14,17 +14,20 @@ This makes the code more robust and future-proof.
 Added
 -----
 - Setup 1D laterals at branches from points and polygons. (PR #81)
-- Adding tests for Mesh workflows. (PR #133)
+- Added tests for Mesh workflows. (PR #133)
+- Added tests for xyz cross-section workflows (PR #153).
 
 Changed
 -------
 - Change default spacing in setup_channels from ``None`` to ``np.inf``. (PR #133)
+- Added ``maxdist`` variable to setup_rivers and setup_channels. (PR #153)
 
 Fixed
 -----
 - Bugfixing of reading of frictions (global), crosssections and boundaries when update. (PR #81)
 - Fixing bug related to changes to pandas TimeDelta formatting, see also https://pandas.pydata.org/docs/whatsnew/v2.2.0.html#other-deprecations. (PR #129)
 - Fixing setup_links1d2d for 2d to 1d direction. (PR #133)
+- Several bugfixes related to processing of cross-sections (PR #153)
 - Support for geopandas v1 (PR #158)
 
 v0.2.0 (20 November 2023)

diff --git a/hydromt_delft3dfm/dflowfm.py b/hydromt_delft3dfm/dflowfm.py
@@ -217,6 +217,7 @@ def setup_channels(
         crosssections_type: str = None,
         spacing: float = np.inf,
         snap_offset: float = 0.0,
+        maxdist: float = 1.0,
         allow_intersection_snapping: bool = True,
     ):
         """Prepare the 1D channels and adds to branches 1D network.
@@ -274,6 +275,10 @@ def setup_channels(
         snap_offset: float, optional
             Snapping tolerance to automatically connecting branches.
             By default 0.0, no snapping is applied.
+        maxdist: float, optional
+            Maximum distance allowed for crosssections to be applied on branches.
+            Only used for `crosssections_type` = point.
+            By default 1.0.
         allow_intersection_snapping: bool, optional
             Switch to choose whether snapping of multiple branch ends are allowed when
             ``snap_offset`` is used.
@@ -344,6 +349,7 @@ def setup_channels(
             region=region,
             crosssections_fn=crosssections_fn,
             crosssections_type=crosssections_type,
+            maxdist=maxdist,
         )
 
         # add crosssections to exisiting ones and update geoms
@@ -668,6 +674,7 @@ def setup_rivers(
         crosssections_fn: Union[int, list] = None,
         crosssections_type: Union[int, list] = None,
         snap_offset: float = 0.0,
+        maxdist: float = 1.0,
         allow_intersection_snapping: bool = True,
     ):
         """Prepare the 1D rivers and adds to 1D branches.
@@ -743,6 +750,10 @@ def setup_rivers(
         snap_offset: float, optional
             Snapping tolerance to automatically connecting branches.
             By default 0.0, no snapping is applied.
+        maxdist: float, optional
+            Maximum distance allowed for crosssections to be applied on branches.
+            Only used for `crosssections_type` = point.
+            By default 1.0.
         allow_intersection_snapping: bool, optional
             Switch to choose whether snapping of multiple branch ends are allowed when
             ``snap_offset`` is used.
@@ -811,6 +822,7 @@ def setup_rivers(
                 region=region,
                 crosssections_fn=crs_fn,
                 crosssections_type=crs_type,
+                maxdist=maxdist,
             )
             crosssections = workflows.add_crosssections(
                 self.geoms.get("crosssections"), crosssections
@@ -1139,6 +1151,7 @@ def _setup_crosssections(
         crosssections_fn: str = None,
         crosssections_type: str = "branch",
         midpoint=True,
+        maxdist=1.0,
     ) -> gpd.GeoDataFrame:
         """Prepare 1D crosssections from branches, points and xyz.
 
@@ -1184,20 +1197,24 @@ def _setup_crosssections(
             * Optional variables:
             If ``crosssections_type`` = "point"
 
-            * Required variables: crsid, shape, shift
+            * Required variables: crsid, shape, shift, closed
             * Optional variables:
-                if shape = 'rectangle': 'width', 'height', 'closed'
-                if shape = 'trapezoid': 'width', 't_width', 'height', 'closed'
-                if shape = 'yz': 'yzcount','ycoordinates','zcoordinates','closed'
+                if shape = 'rectangle': 'width', 'height'
+                if shape = 'trapezoid': 'width', 't_width', 'height'
+                if shape = 'yz': 'yzcount','ycoordinates','zcoordinates'
                 if shape = 'zw': 'numlevels', 'levels', 'flowwidths','totalwidths',
-                'closed'.
+                    'fricitonid', 'frictiontype', 'frictionvalue'
                 if shape = 'zwRiver': Not Supported
                 Note that list input must be strings seperated by a whitespace ''.
             By default None, crosssections will be set from branches
         crosssections_type : {'branch', 'xyz', 'point'}
             Type of crosssections read from crosssections_fn. One of
             ['branch', 'xyz', 'point'].
             By default `branch`.
+        maxdist: float, optional
+            Maximum distance allowed for crosssections to be applied on branches.
+            Only used for `crosssections_type` = point.
+            By default 1.0.
 
         Returns
         -------
@@ -1290,9 +1307,8 @@ def _setup_crosssections(
             # set crsloc and crsdef attributes to crosssections
             self.logger.info(f"Preparing 1D point crossections from {crosssections_fn}")
             gdf_cs = workflows.set_point_crosssections(
-                branches, gdf_cs, maxdist=self._network_snap_offset
+                branches, gdf_cs, maxdist=maxdist
             )
-
         else:
             raise NotImplementedError(
                 f"Method {crosssections_type} is not implemented."

diff --git a/hydromt_delft3dfm/workflows/crosssections.py b/hydromt_delft3dfm/workflows/crosssections.py
@@ -421,8 +421,19 @@ def set_xyz_crosssections(
     crosssections.loc[:, "z"] = crosssections.z
     crosssections.loc[:, "order"] = crosssections.loc[:, "order"].astype("int")
 
+    # count number of points per cross-section id
+    crosssections = crosssections.reset_index(drop=True)
+    xyzcounts = crosssections.copy().groupby("crsid")["crsid"].transform("size")
+    _invalid_ids = crosssections[xyzcounts < 3].index
+    if not _invalid_ids.empty:
+        crosssections = crosssections.drop(_invalid_ids)
+        logger.warning(
+            f"Crosssection with id: {list(_invalid_ids)}"
+            "is dropped: invalid crosssections with less than 3 survey points. "
+        )
+
     # convert xyz crosssection into yz profile
-    crosssections = crosssections.groupby(level=0).apply(xyzp2xyzl, (["order"]))
+    crosssections = crosssections.groupby("crsid").apply(xyzp2xyzl, (["order"]))
     crosssections.crs = branches.crs
 
     # snap to branch
@@ -590,13 +601,19 @@ def set_point_crosssections(
         logger.error("No crossections are set up.")
         return pd.DataFrame()
 
-    # get branch friction
-    crosssections = crosssections.merge(
-        branches[["frictionid", "frictiontype", "frictionvalue"]],
+    # get branch friction (regard crosssections')
+    _friction_cols = ["frictionid", "frictiontype", "frictionvalue"]
+    crosssections = crosssections.drop(
+        columns=[c for c in _friction_cols if c in crosssections.columns],
+    ).merge(
+        branches[_friction_cols],
         left_on="branch_id",
         right_index=True,
     )
 
+    # get "closed" in the correct format
+    crosssections["closed"].replace({1: "yes", 0: "no"}, inplace=True)
+
     # NOTE: below is removed because in case of multiple structures
     # at the same location,
     # there can be multiple crossections
@@ -738,7 +755,8 @@ def set_point_crosssections(
     crosssections_["crsdef_thalweg"] = 0.0
 
     # support both string and boolean for closed column
-    crosssections_["crsdef_closed"].replace({"yes": 1, "no": 0}, inplace=True)
+    if "crsdef_closed" in crosssections_.columns:
+        crosssections_["crsdef_closed"].replace({"yes": 1, "no": 0}, inplace=True)
 
     crosssections_ = gpd.GeoDataFrame(crosssections_, crs=branches.crs)
 
@@ -844,14 +862,21 @@ def _set_trapezoid_crs(crosssections: gpd.GeoDataFrame):
     crsdefs = []
     crslocs = []
     for c in crosssections.itertuples():
-        levels = f"0 {c.height:.6f}"
-        flowwidths = f"{c.width:.6f} {c.t_width:.6f}"
+        # add closed crosssection definition
+        if c.closed == "yes":
+            levels = f"0 {c.height:.6f} {c.height+0.01:.6f}"
+            flowwidths = f"{c.width:.6f} {c.t_width:.6f} 0"
+            numlevels = 3
+        else:
+            levels = f"0 {c.height:.6f}"
+            flowwidths = f"{c.width:.6f} {c.t_width:.6f}"
+            numlevels = 3
         crsdefs.append(
             {
                 "crsdef_id": c.definitionid,
                 "crsdef_type": "zw",
                 "crsdef_branchid": c.branch_id,
-                "crsdef_numlevels": 2,
+                "crsdef_numlevels": numlevels,
                 "crsdef_levels": levels,
                 "crsdef_flowwidths": flowwidths,
                 "crsdef_totalwidths": flowwidths,
@@ -939,7 +964,8 @@ def _set_yz_crs(crosssections: gpd.GeoDataFrame):
                 "crsdef_yzcount": c.yzcount,
                 "crsdef_ycoordinates": c.ycoordinates,
                 "crsdef_zcoordinates": c.zcoordinates,
-                "crsdef_frictionid": c.frictionid,
+                "crsdef_frictionpositions": f"0 { c.ycoordinates.split(' ')[-1]}",
+                "crsdef_frictionids": c.frictionid,
                 "frictiontype": c.frictiontype,
                 "frictionvalue": c.frictionvalue,
             }
@@ -968,20 +994,20 @@ def _set_yz_crs(crosssections: gpd.GeoDataFrame):
     return crosssections_
 
 
-def xyzp2xyzl(xyz: pd.DataFrame, sort_by: list = ["x", "y"]):
+def xyzp2xyzl(xyz: pd.DataFrame, sort_by: list = ["x", "y"]) -> gpd.GeoDataFrame:
     """Convert xyz points to xyz lines.
 
     Parameters
     ----------
     xyz: pd.DataFrame
-        The xyz points.
+        A DataFrame with the xyz-points stored in columns ['x', 'y', 'z'].
     sort_by: list, optional
         List of attributes to sort by. Defaults to ["x", "y"].
 
     Returns
     -------
-    gpd.GeoSeries
-        The xyz lines.
+    gpd.GeoDataframe
+        A GeoDataframe with the xy-profile as a LineString and a column for z.
     """
     sort_by = [s.lower() for s in sort_by]
 

diff --git a/pyproject.toml b/pyproject.toml
@@ -9,6 +9,7 @@ authors = [
 	{ name = "Xiaohan Li", email = "[email protected]" },
 	{ name = "Hélène Boisgontier", email = "[email protected]" },
 	{ name = "Jelmer Veenstra", email = "[email protected]" },
+ 	{ name = "Sebastian Hartgring", email = "[email protected]" },
 ]
 dependencies = [
 	"hydromt>=0.10.0, <1", # TODO: move to hydromt>=1: https://github.com/Deltares/hydromt_delft3dfm/issues/137

diff --git a/tests/test_workflows_crosssections.py b/tests/test_workflows_crosssections.py
@@ -0,0 +1,29 @@
+from os.path import abspath, dirname, join
+import geopandas as gpd
+import pandas as pd
+import hydromt_delft3dfm.workflows.crosssections as xsec # needed to access private functions
+from shapely import Point, LineString
+
+EXAMPLEDIR = join(dirname(abspath(__file__)), "..", "examples")
+
+
+def test_set_xyz_crosssections():
+    # test on a single branch with 3 cross-sections (id = 1 .. 3)
+    branches = gpd.GeoDataFrame(data={"frictionid": ["Manning_0.023"], "frictiontype": ["Manning"], "frictionvalue": [0.023]}, geometry=[LineString([[0, 0], [1000, 1000]])], crs=28992)
+
+    xyz_crosssections_1 = gpd.GeoDataFrame(
+        geometry=[Point(-10, 1), Point(0, 1), Point(10, 2), Point(20, 2)],
+        data={"crsid": [1, 1, 1, 1], "order": [1, 2, 3, 4], "z": [10, 0, 1, 11]}, crs=28992)
+    xyz_crosssections_2 = gpd.GeoDataFrame(
+        geometry=[Point(890, 900), Point(900, 895), Point(910, 890), Point(920, 890)],
+        data={"crsid": [2, 2, 2, 2], "order": [1, 2, 3, 4], "z": [9, -1, 0, 10]}, crs=28992)
+    # cross-section 3 deliberately has too few points
+    xyz_crosssections_3 = gpd.GeoDataFrame(
+        geometry=[Point(995, 995), Point(1005, 1005)],
+        data={"crsid": [3, 3], "order": [1, 2], "z": [8, 8]}, crs=28992)
+
+    xyz_crosssections = pd.concat([xyz_crosssections_1, xyz_crosssections_2, xyz_crosssections_3]).reset_index(drop=True)
+    crosssections = xsec.set_xyz_crosssections(branches=branches, crosssections=xyz_crosssections)
+
+    assert len(crosssections) == 2 # check if the cross-section containing only 2 points is dropped
+    #TODO more checks can be added