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

Adding vegetation height to table_yz and single

threedigrid_builder/base/lines.py

@@ -38,8 +38,10 @@ 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
+    veg_coef1: float  # the product of vegetation properties (except height)
+    veg_coef2: float  # the product of vegetation properties (except height)
+    veg_height1: float  # the vegetaion height
+    veg_height2: float  # the vegetaion height
     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

@@ -94,7 +94,7 @@ def convert(self, ids):
         if len(tables_yz) > 0:
             result.tables_yz = np.concatenate(tables_yz, axis=0)
         else:
-            result.tables_yz = np.empty((0, 4))
+            result.tables_yz = np.empty((0, 5))
 
         return result
 
@@ -304,7 +304,7 @@ def tabulate_yz(shape, width, height):
         yz = None
         shape_return = CrossSectionShape.TABULATED_TRAPEZIUM
     else:
-        yz = np.zeros((len(ys), 4), dtype=float)
+        yz = np.zeros((len(ys), 5), dtype=float)
         yz[:, 0] = ys
         yz[:, 1] = zs
         shape_return = CrossSectionShape.TABULATED_YZ
@@ -385,6 +385,7 @@ def set_friction_vegetation_values(
         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
+        yz[:-1, 3] = veg_stemden * veg_stemdia * veg_drag
+        yz[:-1, 4] = veg_hght
 
     return yz

threedigrid_builder/grid/cross_section_locations.py

@@ -45,8 +45,10 @@ 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
+        - veg_coef1: the product of vegetation properties of the first cross section location (except height)
+        - veg_coef2: the product of vegetation properties of the second cross section location (except height)
+        - veg_height1: the vegetation height of the first cross section
+        - veg_height2: the vegetation height of the second cross section
         - 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
@@ -74,18 +76,20 @@ def apply_to_lines(self, lines, channels, extrapolate=False):
         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
+        lines.veg_height1 = self.vegetation_height[idx1]
+        lines.veg_height1[np.isnan(lines.veg_height1)] = 0.0
+        lines.veg_height2 = self.vegetation_height[idx2]
+        lines.veg_height2[np.isnan(lines.veg_height2)] = 0.0
 
         # Compute invert levels and start and end
         lines.invert_level_start_point = compute_bottom_level(

threedigrid_builder/interface/gridadmin.py

@@ -444,6 +444,8 @@ def write_lines(self, lines: Lines, cross_sections):
         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, "veg_height1", lines.veg_height1)
+        self.write_dataset(group, "veg_height2", lines.veg_height2)
         self.write_dataset(group, "cross_weight", lines.cross_weight)
         self.write_dataset(
             group, "hydraulic_resistance_in", lines.hydraulic_resistance_in

threedigrid_builder/tests/test_cross_section_definitions.py

@@ -213,7 +213,12 @@ def test_tabulate_inverted_egg():
             1.5,
             0.5,
             [[0, 0.5], [0.5, 1.5]],
-            [[0, 0.5, 0, 0], [0.5, 0, 0, 0], [1, 0, 0, 0], [1.5, 0.5, 0, 0]],
+            [
+                [0, 0.5, 0, 0, 0],
+                [0.5, 0, 0, 0, 0],
+                [1, 0, 0, 0, 0],
+                [1.5, 0.5, 0, 0, 0],
+            ],
         ),
         (
             "0 0.5 1 1.5",
@@ -226,7 +231,12 @@ def test_tabulate_inverted_egg():
             1.5,
             0.5,
             [[0, 0.5], [0.25, 1.25], [0.5, 1.5]],
-            [[0, 0.5, 1, 3], [0.5, 0, 1, 1], [1, 0, 1, 1], [1.5, 0.25, 0, 0]],
+            [
+                [0, 0.5, 1, 1, 3],
+                [0.5, 0, 1, 1, 1],
+                [1, 0, 1, 1, 1],
+                [1.5, 0.25, 0, 0, 0],
+            ],
         ),
         (
             "0 1 2 3 4 5",
@@ -240,12 +250,12 @@ def test_tabulate_inverted_egg():
             1,
             [[0, 0], [0.5, 3], [0.5, 4], [1, 5]],
             [
-                [0, 1, 1, 1],
-                [1, 0, 1, 1],
-                [2, 0.5, 1, 1],
-                [3, 0.5, 1, 1],
-                [4, 0, 1, 1],
-                [5, 1, 0, 0],
+                [0, 1, 1, 1, 1],
+                [1, 0, 1, 1, 1],
+                [2, 0.5, 1, 1, 1],
+                [3, 0.5, 1, 1, 1],
+                [4, 0, 1, 1, 1],
+                [5, 1, 0, 0, 0],
             ],
         ),
         (
@@ -273,11 +283,11 @@ def test_tabulate_inverted_egg():
             0.5,
             [[0, 0.5], [0.5, 1.5]],
             [
-                [0, 0.5, 1, 2],
-                [0.5, 0, 1, 0.1],
-                [0.75, 0, 1, 1],
-                [1, 0, 1, 1],
-                [1.5, 0.5, 0, 0],
+                [0, 0.5, 1, 2, 1],
+                [0.5, 0, 1, 0.1, 1],
+                [0.75, 0, 1, 1, 1],
+                [1, 0, 1, 1, 1],
+                [1.5, 0.5, 0, 0, 0],
             ],
         ),
         (

threedigrid_builder/tests/test_cross_section_locations.py

@@ -169,8 +169,10 @@ def test_apply_to_lines(channels, channel_lines, locations):
     assert_equal(channel_lines.frict_type2, [1, 2, 1, 1, 1, 2, 2])
     assert_equal(channel_lines.frict_value1, [30.0, 0.02, 40.0, 40.0, 40.0, 0.03, 0.03])
     assert_equal(channel_lines.frict_value2, [30.0, 0.02, 35.0, 35.0, 35.0, 0.03, 0.03])
-    assert_almost_equal(channel_lines.veg_coef1, [0.1, 0.0, 0.2, 0.2, 0.2, 0.0, 0.0])
-    assert_almost_equal(channel_lines.veg_coef2, [0.1, 0.0, 0.15, 0.15, 0.15, 0.0, 0.0])
+    assert_almost_equal(channel_lines.veg_coef1, [0.5, 0.0, 0.2, 0.2, 0.2, 0.0, 0.0])
+    assert_almost_equal(channel_lines.veg_coef2, [0.5, 0.0, 0.15, 0.15, 0.15, 0.0, 0.0])
+    assert_almost_equal(channel_lines.veg_height1, [0.2, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0])
+    assert_almost_equal(channel_lines.veg_height2, [0.2, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0])
     assert_almost_equal(
         channel_lines.cross_weight, [1.0, 1.0, 1.0, 0.65, 0.0, 1.0, 1.0]
     )