Skip to content

Commit

Permalink
Update shape in yz & unit test, fric/veg function
Browse files Browse the repository at this point in the history
  • Loading branch information
Golnesa committed Jan 30, 2024
1 parent 322728f commit f48fe07
Show file tree
Hide file tree
Showing 6 changed files with 100 additions and 49 deletions.
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
4 changes: 2 additions & 2 deletions threedigrid_builder/base/lines.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
67 changes: 41 additions & 26 deletions threedigrid_builder/grid/cross_section_definitions.py
Original file line number Diff line number Diff line change
Expand Up @@ -56,24 +56,29 @@ 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],
result.width_1d[i],
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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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:
Expand All @@ -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 "
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 = {
Expand All @@ -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
71 changes: 53 additions & 18 deletions threedigrid_builder/tests/test_cross_section_definitions.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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(
Expand All @@ -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))
2 changes: 1 addition & 1 deletion threedigrid_builder/tests/test_db.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
3 changes: 2 additions & 1 deletion threedigrid_builder/tests/test_gridadmin.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"),
Expand Down

0 comments on commit f48fe07

Please sign in to comment.