diff --git a/setup.py b/setup.py index d301f23e..575e9999 100644 --- a/setup.py +++ b/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.218.*", "shapely>=2", "pyproj>=3", "condenser[geo]>=0.1.1", diff --git a/threedigrid_builder/base/lines.py b/threedigrid_builder/base/lines.py index d52aa25a..8e0123c2 100644 --- a/threedigrid_builder/base/lines.py +++ b/threedigrid_builder/base/lines.py @@ -38,8 +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 - veg_coef2: float + 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 diff --git a/threedigrid_builder/grid/cross_section_definitions.py b/threedigrid_builder/grid/cross_section_definitions.py index 79e1db94..8fa63640 100644 --- a/threedigrid_builder/grid/cross_section_definitions.py +++ b/threedigrid_builder/grid/cross_section_definitions.py @@ -56,11 +56,8 @@ def convert(self, ids): 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], @@ -68,12 +65,20 @@ def convert(self, ids): result.height_1d[i], table, yz, - ) = tabulator(shape, width_idx[i], height_idx[i]) + ) = tabulator(shape, self.width[self_i], self.height[self_i]) if table is not None: result.count[i] = 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) @@ -180,7 +185,7 @@ def tabulate_inverted_egg(shape, width, height): See tabulate_egg. """ - type_, width, height, table, yz = tabulate_egg(shape, width, height) + type_, width, height, table, _ = tabulate_egg(shape, width, height) table[:, 1] = table[::-1, 1] return type_, width, height, table, None @@ -253,16 +258,7 @@ def tabulate_tabulated(shape, width, height): return shape, np.max(widths), np.max(heights), np.array([heights, widths]).T, None -def tabulate_yz( - shape, - width, - height, - friction_values, - vegetation_stem_densities, - vegetation_stem_diameters, - vegetation_heights, - vegetation_drag_coefficients, -): +def tabulate_yz(shape, width, height): """Tabulate an (open or closed) YZ profile Args: @@ -272,10 +268,11 @@ def tabulate_yz( 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] + shape_return = CrossSectionShape.TABULATED_TRAPEZIUM if is_closed and len(zs) < 4: raise SchematisationError( f"Cross section definitions of closed profiles must have at least " @@ -310,13 +307,7 @@ def tabulate_yz( yz = np.zeros((len(ys), 4), dtype=float) yz[:, 0] = ys yz[:, 1] = zs - fric = np.array([float(x) for x in friction_values.split(" ")]) - 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, 2] = fric - yz[:-1, 3] = veg_stemden * veg_stemdia * veg_hght * veg_drag + 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 @@ -360,7 +351,7 @@ def tabulate_yz( table = table[1:] height_1d, width_1d = table.max(axis=0).tolist() - return CrossSectionShape.TABULATED_TRAPEZIUM, width_1d, height_1d, table, yz + return shape_return, width_1d, height_1d, table, yz tabulators = { @@ -373,3 +364,27 @@ def tabulate_yz( 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 len(friction_values) > 0: + fric = np.array([float(x) for x in friction_values.split(" ")]) + yz[:-1, 2] = fric + + if len(vegetation_drag_coefficients) > 0: + 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 diff --git a/threedigrid_builder/tests/test_cross_section_definitions.py b/threedigrid_builder/tests/test_cross_section_definitions.py index 99efe7fa..8b8ac2c0 100644 --- a/threedigrid_builder/tests/test_cross_section_definitions.py +++ b/threedigrid_builder/tests/test_cross_section_definitions.py @@ -8,6 +8,7 @@ from threedigrid_builder.exceptions import SchematisationError from threedigrid_builder.grid import CrossSectionDefinitions, CrossSections from threedigrid_builder.grid.cross_section_definitions import ( + set_friction_vegetation_values, tabulate_builtin, tabulate_closed_rectangle, tabulate_egg, @@ -204,15 +205,15 @@ def test_tabulate_inverted_egg(): ( "0 0.5 1 1.5", "0.5 0 0 0.5", - "1 1 1", - "1 1 1", - "1 1 1", - "1 1 1", - "1 1 1", + "", + "", + "", + "", + "", 1.5, 0.5, [[0, 0.5], [0.5, 1.5]], - [[0, 0.5, 1, 1], [0.5, 0, 1, 1], [1, 0, 1, 1], [1.5, 0.5, 0, 0]], + [[0, 0.5, 0, 0], [0.5, 0, 0, 0], [1, 0, 0, 0], [1.5, 0.5, 0, 0]], ), ( "0 0.5 1 1.5", @@ -264,16 +265,16 @@ def test_tabulate_inverted_egg(): "0 0.5 0.75 1.0 1.5", "0.5 0 0 0 0.5", "1 1 1 1", - "1 1 1 1", - "1 1 1 1", + "1 0.1 1 1", + "2 1 1 1", "1 1 1 1", "1 1 1 1", 1.5, 0.5, [[0, 0.5], [0.5, 1.5]], [ - [0, 0.5, 1, 1], - [0.5, 0, 1, 1], + [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], @@ -311,19 +312,24 @@ def test_tabulate_yz( "my-shape", width, height, - friction_values, - vegetation_stem_densities, - vegetation_stem_diameters, - vegetation_heights, - vegetation_drag_coefficients, ) - assert shape == CrossSectionShape.TABULATED_TRAPEZIUM assert width_1d == exp_width assert height_1d == exp_height assert_almost_equal(table, np.array(exp_table, dtype=float)) if yz is not None: - assert_almost_equal(yz, np.array(exp_yz, dtype=float)) + yz_1d = set_friction_vegetation_values( + yz, + friction_values, + vegetation_stem_densities, + vegetation_stem_diameters, + vegetation_heights, + vegetation_drag_coefficients, + ) + assert_almost_equal(yz_1d, np.array(exp_yz, dtype=float)) + assert shape == CrossSectionShape.TABULATED_YZ + else: + assert shape == CrossSectionShape.TABULATED_TRAPEZIUM @pytest.mark.parametrize( @@ -347,4 +353,33 @@ def test_tabulate_yz( ) def test_tabulate_yz_err(width, height, match): with pytest.raises(SchematisationError, match=match): - tabulate_yz("my-shape", width, height, None, None, None, None, None) + tabulate_yz("my-shape", width, height) + + +# @pytest.mark.parametrize( +# "yz,friction_values,vegetation_stem_densities,vegetation_stem_diameters,vegetation_heights,vegetation_drag_coefficients,exp_yz", +# [ +# ( +# [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], +# "0.2 0.3 0.1", +# "1 2 3", +# "0.5 0.1 0.3", +# "2 3 1", +# "1 1 1", +# [ +# [0, 0, 0.2, 1], +# [0, 0, 0.3, 0.6], +# [0, 0, 0.1, 0.9], +# [0, 0, 0, 0], +# ], +# ), +# ], +# ) +# def test_set_friction_vegetation_values(yz,friction_values,vegetation_stem_densities,vegetation_stem_diameters,vegetation_heights,vegetation_drag_coefficients,exp_yz): +# yz_1d = set_friction_vegetation_values(yz, +# friction_values, +# vegetation_stem_densities, +# vegetation_stem_diameters, +# vegetation_heights, +# vegetation_drag_coefficients) +# assert_almost_equal(yz_1d, np.array(exp_yz, dtype=float)) diff --git a/threedigrid_builder/tests/test_db.py b/threedigrid_builder/tests/test_db.py index 7f329c9d..435d2202 100644 --- a/threedigrid_builder/tests/test_db.py +++ b/threedigrid_builder/tests/test_db.py @@ -65,7 +65,7 @@ def test_init_bad_version(tmp_path): def test_get_version(db): - assert db.get_version() == 217 + assert db.get_version() == 218 def test_get_boundary_conditions_1d(db): diff --git a/threedigrid_builder/tests/test_gridadmin.py b/threedigrid_builder/tests/test_gridadmin.py index 2f1dc334..ff8bbc45 100644 --- a/threedigrid_builder/tests/test_gridadmin.py +++ b/threedigrid_builder/tests/test_gridadmin.py @@ -140,7 +140,8 @@ def test_write_nodes_embedded(h5_out, dataset, shape, dtype): ("dist_calc_points", (6,), "float64"), ("friction_type", (6,), "int32"), ("friction_value", (6,), "float64"), - # ("vegetation_stem_density", (6,), "float64"), + ("veg_coef1", (6,), "float64"), + ("veg_coef2", (6,), "float64"), ("id", (6,), "int32"), ("invert_level_end_point", (6,), "float64"), ("invert_level_start_point", (6,), "float64"),