Merge pull request #346 from nens/golnesa-variable-friction

Golnesa variable friction

.github/workflows/test.yml

@@ -95,4 +95,4 @@ jobs:
       - name: Run integration tests
         shell: bash
         run: |
-          pytest -c /dev/null integration_tests
+          pytest integration_tests

CHANGES.rst

@@ -4,7 +4,13 @@ Changelog of threedigrid-builder
 1.12.3 (unreleased)
 -------------------
 
-- Nothing changed yet.
+- Add single vegetation and variable friction/vegetation for 1D elements. 
+
+- Add a table for open YZ profile to support 1D variable friction and vegetation.
+
+- Update unit tests and the integration test (bergemeer test).
+
+- Remove the spatialie4 version of the integration test (bergemeer test).
 
 
 1.12.2 (2023-12-07)

integration_tests/test_bergermeer.py

@@ -20,9 +20,7 @@ def count_unique(arr):
     return dict(zip(*np.unique(arr, return_counts=True)))
 
 
-@pytest.mark.parametrize(
-    "filename", ["v2_bergermeer.sqlite", "v2_bergermeer_spatialite4.sqlite"]
-)
+@pytest.mark.parametrize("filename", ["v2_bergermeer.sqlite"])
 def test_integration(tmp_path, filename):
     shutil.copyfile(unittests_data_path / filename, tmp_path / filename)
     progress_callback = Mock()
@@ -108,6 +106,8 @@ def test_integration(tmp_path, filename):
             -9999: 1,
             CrossSectionShape.CIRCLE: 8,
             CrossSectionShape.RECTANGLE: 3,
+            CrossSectionShape.TABULATED_TRAPEZIUM: 1,
+            CrossSectionShape.TABULATED_YZ: 1,
         }
 
         ## EMBEDDED NODES

setup.py

@@ -52,7 +52,7 @@ def get_version():
 
 install_requires = [
     "numpy>=1.15,<1.25.0",
-    "threedi-schema>=0.217.0",
+    "threedi-schema==0.219.*",
     "shapely>=2",
     "pyproj>=3",
     "condenser[geo]>=0.1.1",

threedigrid_builder/base/lines.py

@@ -38,6 +38,8 @@ class Line:
     frict_type2: int
     frict_value1: float  # For 1D-2Dgw: resistance factor-out (hydr. cond / thickness)
     frict_value2: float  # For 1D-2Dgw: resistance factor-in (hydr. cond / thickness)
+    veg_coef1: float  # the product of vegetation properties
+    veg_coef2: float  # the product of vegetation properties
     cross_weight: float  # For 1D-2Dgw: the exchange length (tables: cross_width)
     discharge_coefficient_positive: float
     discharge_coefficient_negative: float

threedigrid_builder/grid/cross_section_definitions.py

@@ -17,6 +17,11 @@ class CrossSectionDefinition:
     shape: CrossSectionShape
     height: str  # space-separated list of floats
     width: str  # space-separated list of floats
+    friction_values: str  # space-separated list of floats
+    vegetation_stem_densities: str  # space-separated list of floats
+    vegetation_stem_diameters: str  # space-separated list of floats
+    vegetation_heights: str  # space-separated list of floats
+    vegetation_drag_coefficients: str  # space-separated list of floats
 
 
 class CrossSectionDefinitions(Array[CrossSectionDefinition]):
@@ -43,36 +48,54 @@ def convert(self, ids):
             content_pk=ids,
             code=self.code[idx],
             count=0,
+            count_yz=0,
         )
         if len(result) == 0:
             return result
 
-        offset = 0
         tables = []
+        tables_yz = []
 
-        # Numpy array views on width/height based on idx
-        width_idx = self.width[idx]
-        height_idx = self.height[idx]
-
-        for i, shape in enumerate(self.shape[idx]):
+        for i, self_i in enumerate(idx):
+            shape = self.shape[self_i]
             tabulator = tabulators[shape]
-            result.shape[i], result.width_1d[i], result.height_1d[i], table = tabulator(
-                shape, width_idx[i], height_idx[i]
-            )
+            (
+                result.shape[i],
+                result.width_1d[i],
+                result.height_1d[i],
+                table,
+                yz,
+            ) = tabulator(shape, self.width[self_i], self.height[self_i])
             if table is not None:
                 result.count[i] = len(table)
-                result.offset[i] = offset
-                offset += len(table)
                 tables.append(table)
+            if yz is not None:
+                result.count_yz[i] = len(yz)
+                yz = set_friction_vegetation_values(
+                    yz,
+                    self.friction_values[self_i],
+                    self.vegetation_stem_densities[self_i],
+                    self.vegetation_stem_diameters[self_i],
+                    self.vegetation_heights[self_i],
+                    self.vegetation_drag_coefficients[self_i],
+                )
+                tables_yz.append(yz)
 
         result.offset[:] = np.roll(np.cumsum(result.count), 1)
         result.offset[0] = 0
+        result.offset_yz[:] = np.roll(np.cumsum(result.count_yz), 1)
+        result.offset_yz[0] = 0
 
         if len(tables) > 0:
             result.tables = np.concatenate(tables, axis=0)
         else:
             result.tables = np.empty((0, 2))
 
+        if len(tables_yz) > 0:
+            result.tables_yz = np.concatenate(tables_yz, axis=0)
+        else:
+            result.tables_yz = np.empty((0, 4))
+
         return result
 
 
@@ -85,11 +108,14 @@ class CrossSection:
     height_1d: float
     offset: int
     count: int
+    offset_yz: int
+    count_yz: int
     # tables: Tuple[float, float] has different length so is specified on CrossSections
 
 
 class CrossSections(Array[CrossSection]):
     tables = None
+    tables_yz = None
 
 
 def tabulate_builtin(shape, width, height):
@@ -103,7 +129,7 @@ def tabulate_builtin(shape, width, height):
         height (str): ignored
 
     Returns:
-        tuple:  shape, width_1d (float), None, None
+        tuple:  shape, width_1d (float), None, None, None
     """
     try:
         width = float(width)
@@ -112,7 +138,7 @@ def tabulate_builtin(shape, width, height):
             f"Unable to parse cross section definition width (got: '{width}')."
         )
 
-    return shape, width, None, None
+    return shape, width, None, None, None
 
 
 def tabulate_egg(shape, width, height):
@@ -125,7 +151,7 @@ def tabulate_egg(shape, width, height):
 
     Returns:
         tuple:  TABULATED_TRAPEZIUM, width_1d (float),
-                height_1d (float), table (ndarray of shape (M, 2))
+                height_1d (float), table (ndarray of shape (M, 2)), None
     """
     NUM_INCREMENTS = 16
 
@@ -151,17 +177,17 @@ def tabulate_egg(shape, width, height):
     widths = np.sqrt(p / q) * 2
 
     table = np.array([heights, widths]).T
-    return CrossSectionShape.TABULATED_TRAPEZIUM, width, height, table
+    return CrossSectionShape.TABULATED_TRAPEZIUM, width, height, table, None
 
 
 def tabulate_inverted_egg(shape, width, height):
     """Tabulate the egg shape, upside down.
 
     See tabulate_egg.
     """
-    type_, width, height, table = tabulate_egg(shape, width, height)
+    type_, width, height, table, _ = tabulate_egg(shape, width, height)
     table[:, 1] = table[::-1, 1]
-    return type_, width, height, table
+    return type_, width, height, table, None
 
 
 def tabulate_closed_rectangle(shape, width, height):
@@ -174,7 +200,7 @@ def tabulate_closed_rectangle(shape, width, height):
 
     Returns:
         tuple:  TABULATED_RECTANGLE, width_1d (float),
-                height (float), table (ndarray of shape (M, 2))
+                height (float), table (ndarray of shape (M, 2)), None
     """
     try:
         width = float(width)
@@ -185,7 +211,7 @@ def tabulate_closed_rectangle(shape, width, height):
             f"(got: '{width}', '{height}')."
         )
     table = np.array([[0.0, width], [height, 0.0]], order="F")
-    return CrossSectionShape.TABULATED_RECTANGLE, width, height, table
+    return CrossSectionShape.TABULATED_RECTANGLE, width, height, table, None
 
 
 def _parse_tabulated(width, height):
@@ -220,7 +246,7 @@ def tabulate_tabulated(shape, width, height):
 
     Returns:
         tuple:  shape, width_1d (float),
-                height_1d (float), table (ndarray of shape (M, 2))
+                height_1d (float), table (ndarray of shape (M, 2)), None
     """
     widths, heights = _parse_tabulated(width, height)
     if len(heights) > 1 and np.any(np.diff(heights) < 0.0):
@@ -229,7 +255,7 @@ def tabulate_tabulated(shape, width, height):
             f"(got: {height})."
         )
 
-    return shape, np.max(widths), np.max(heights), np.array([heights, widths]).T
+    return shape, np.max(widths), np.max(heights), np.array([heights, widths]).T, None
 
 
 def tabulate_yz(shape, width, height):
@@ -242,7 +268,7 @@ def tabulate_yz(shape, width, height):
 
     Returns:
         tuple:  shape, width_1d (float),
-                height_1d (float), table (ndarray of shape (M, 2))
+                height_1d (float), table (ndarray of shape (M, 2)), yz (ndarray of shape (M, 4))
     """
     ys, zs = _parse_tabulated(width, height)
     is_closed = ys[0] == ys[-1] and zs[0] == zs[-1]
@@ -275,6 +301,13 @@ def tabulate_yz(shape, width, height):
     if is_closed:
         ys = ys[:-1]
         zs = zs[:-1]
+        yz = None
+        shape_return = CrossSectionShape.TABULATED_TRAPEZIUM
+    else:
+        yz = np.zeros((len(ys), 4), dtype=float)
+        yz[:, 0] = ys
+        yz[:, 1] = zs
+        shape_return = CrossSectionShape.TABULATED_YZ
 
     # Adapt non-unique height coordinates. Why?
     # Because if a segment of the profile is exactly horizontal, we need 2 widths
@@ -318,7 +351,7 @@ def tabulate_yz(shape, width, height):
         table = table[1:]
 
     height_1d, width_1d = table.max(axis=0).tolist()
-    return CrossSectionShape.TABULATED_TRAPEZIUM, width_1d, height_1d, table
+    return shape_return, width_1d, height_1d, table, yz
 
 
 tabulators = {
@@ -331,3 +364,27 @@ def tabulate_yz(shape, width, height):
     CrossSectionShape.TABULATED_YZ: tabulate_yz,
     CrossSectionShape.INVERTED_EGG: tabulate_inverted_egg,
 }
+
+
+def set_friction_vegetation_values(
+    yz,
+    friction_values,
+    vegetation_stem_densities,
+    vegetation_stem_diameters,
+    vegetation_heights,
+    vegetation_drag_coefficients,
+):
+    """Convert friction and vegetation properties from list into arrays, if available,
+    and add to yz"""
+    if friction_values is not None:
+        fric = np.array([float(x) for x in friction_values.split(" ")])
+        yz[:-1, 2] = fric
+
+    if vegetation_drag_coefficients is not None:
+        veg_stemden = np.array([float(x) for x in vegetation_stem_densities.split(" ")])
+        veg_stemdia = np.array([float(x) for x in vegetation_stem_diameters.split(" ")])
+        veg_hght = np.array([float(x) for x in vegetation_heights.split(" ")])
+        veg_drag = np.array([float(x) for x in vegetation_drag_coefficients.split(" ")])
+        yz[:-1, 3] = veg_stemden * veg_stemdia * veg_hght * veg_drag
+
+    return yz

threedigrid_builder/grid/cross_section_locations.py

@@ -17,6 +17,10 @@ class CrossSectionLocation:
     bank_level: float
     friction_type: FrictionType
     friction_value: float
+    vegetation_stem_density: float
+    vegetation_stem_diameter: float
+    vegetation_height: float
+    vegetation_drag_coefficient: float
 
 
 class CrossSectionLocations(Array[CrossSectionLocation]):
@@ -41,6 +45,8 @@ def apply_to_lines(self, lines, channels, extrapolate=False):
         - frict_type2: the friction type of the second cross section location
         - frict_value1: the friction value of the first cross section location
         - frict_value2: the friction value of the second cross section location
+        - veg_coef1: the product of vegetation properties of the first cross section location
+        - veg_coef2: the product of vegetation properties of the second cross section location
         - invert_level_start_point: 'reference_level' interpolated at the line end
         - invert_level_end_point: 'reference_level' interpolated at the line start
         - dpumax: the largest of the two invert levels
@@ -65,6 +71,22 @@ def apply_to_lines(self, lines, channels, extrapolate=False):
         lines.frict_value1 = self.friction_value[idx1]
         lines.frict_value2 = self.friction_value[idx2]
 
+        lines.veg_coef1 = (
+            self.vegetation_stem_density[idx1]
+            * self.vegetation_stem_diameter[idx1]
+            * self.vegetation_height[idx1]
+            * self.vegetation_drag_coefficient[idx1]
+        )
+        lines.veg_coef1[np.isnan(lines.veg_coef1)] = 0.0
+
+        lines.veg_coef2 = (
+            self.vegetation_stem_density[idx2]
+            * self.vegetation_stem_diameter[idx2]
+            * self.vegetation_height[idx2]
+            * self.vegetation_drag_coefficient[idx2]
+        )
+        lines.veg_coef2[np.isnan(lines.veg_coef2)] = 0.0
+
         # Compute invert levels and start and end
         lines.invert_level_start_point = compute_bottom_level(
             lines.content_pk,

threedigrid_builder/interface/db.py

@@ -468,6 +468,11 @@ def get_cross_section_definitions(self) -> CrossSectionDefinitions:
                     models.CrossSectionDefinition.shape,
                     models.CrossSectionDefinition.width,
                     models.CrossSectionDefinition.height,
+                    models.CrossSectionDefinition.friction_values,
+                    models.CrossSectionDefinition.vegetation_stem_densities,
+                    models.CrossSectionDefinition.vegetation_stem_diameters,
+                    models.CrossSectionDefinition.vegetation_heights,
+                    models.CrossSectionDefinition.vegetation_drag_coefficients,
                 )
                 .order_by(models.CrossSectionDefinition.id)
                 .as_structarray()
@@ -497,6 +502,10 @@ def get_cross_section_locations(self) -> CrossSectionLocations:
                     models.CrossSectionLocation.bank_level,
                     models.CrossSectionLocation.friction_type,
                     models.CrossSectionLocation.friction_value,
+                    models.CrossSectionLocation.vegetation_stem_density,
+                    models.CrossSectionLocation.vegetation_stem_diameter,
+                    models.CrossSectionLocation.vegetation_height,
+                    models.CrossSectionLocation.vegetation_drag_coefficient,
                 )
                 .order_by(models.CrossSectionLocation.id)
                 .as_structarray()

threedigrid_builder/interface/gridadmin.py

@@ -442,6 +442,8 @@ def write_lines(self, lines: Lines, cross_sections):
         self.write_dataset(group, "frict_type2", lines.frict_type2)
         self.write_dataset(group, "frict_value1", lines.frict_value1)
         self.write_dataset(group, "frict_value2", lines.frict_value2)
+        self.write_dataset(group, "veg_coef1", lines.veg_coef1)
+        self.write_dataset(group, "veg_coef2", lines.veg_coef2)
         self.write_dataset(group, "cross_weight", lines.cross_weight)
         self.write_dataset(
             group, "hydraulic_resistance_in", lines.hydraulic_resistance_in
@@ -690,6 +692,12 @@ def write_cross_sections(self, cross_sections: CrossSections):
         self.write_dataset(group, "offset", cross_sections.offset)
         self.write_dataset(group, "count", cross_sections.count)
         self.write_dataset(group, "tables", cross_sections.tables.T, insert_dummy=False)
+        if cross_sections.tables_yz is not None:
+            self.write_dataset(group, "offset_yz", cross_sections.offset_yz)
+            self.write_dataset(group, "count_yz", cross_sections.count_yz)
+            self.write_dataset(
+                group, "tables_yz", cross_sections.tables_yz.T, insert_dummy=False
+            )
 
     def write_obstacles(self, obstacles: Obstacles):
         """For backwards compat, the group is named 'levees'"""