diff --git a/integration_tests/test_bergermeer.py b/integration_tests/test_bergermeer.py index d3c807e9..9843d32c 100644 --- a/integration_tests/test_bergermeer.py +++ b/integration_tests/test_bergermeer.py @@ -63,8 +63,6 @@ def test_integration(tmp_path, filename): NodeType.NODE_1D_BOUNDARIES: 4, } assert np.count_nonzero(f["nodes"]["is_manhole"][:] == 1) == 42 - is_manhole = f["nodes"]["is_manhole"][:] == 1 - assert np.all(f["nodes"]["zoom_category"][is_manhole] != -9999) assert np.count_nonzero(f["nodes"]["content_pk"][:] > 0) == 1360 assert ( np.count_nonzero(np.isfinite(f["nodes"]["initial_waterlevel"][:])) == 2531 diff --git a/setup.py b/setup.py index 246393da..c5847ad8 100644 --- a/setup.py +++ b/setup.py @@ -53,7 +53,7 @@ def get_version(): install_requires = [ "numpy>=1.15,<3.0", - "threedi-schema==0.227.*", + "threedi-schema==0.228.*", "shapely>=2", "pyproj>=3", "condenser[geo]>=0.1.1", diff --git a/threedigrid_builder/application.py b/threedigrid_builder/application.py index 64fa2651..fd2ae0ba 100644 --- a/threedigrid_builder/application.py +++ b/threedigrid_builder/application.py @@ -22,6 +22,7 @@ GridAdminOut, SQLite, ) +from threedigrid_builder.interface.db import map_cross_section_definition __all__ = ["make_grid", "make_gridadmin"] @@ -87,7 +88,7 @@ def _make_gridadmin( grid.meta.has_groundwater_flow, node_id_counter, line_id_counter ) - grid.set_obstacles(db.get_obstacles()) + grid.set_obstacles_2d(db.get_obstacles()) grid.set_boundary_conditions_2d( db.get_boundary_conditions_2d(), quadtree, @@ -101,6 +102,20 @@ def _make_gridadmin( connection_nodes = db.get_connection_nodes() if grid_settings.use_1d_flow and len(connection_nodes) > 0: progress_callback(0.8, "Constructing 1D computational grid...") + locations = db.get_cross_section_locations() + pipes = db.get_pipes() + weirs = db.get_weirs() + culverts = db.get_culverts() + orifices = db.get_orifices() + cross_section_definitions = db.get_cross_section_definitions() + ( + cross_section_definitions_unique, + mapping, + ) = cross_section_definitions.get_unique() + map_cross_section_definition( + [locations, orifices, pipes, culverts, weirs], mapping + ) + cn_grid = Grid.from_connection_nodes( connection_nodes=connection_nodes, node_id_counter=node_id_counter, @@ -109,6 +124,7 @@ def _make_gridadmin( grid += cn_grid channels = db.get_channels() + potential_breaches = db.get_potential_breaches() breach_points = potential_breaches.project_on_channels(channels).merge() @@ -125,13 +141,11 @@ def _make_gridadmin( connection_node_offset=connection_node_first_id, ) - locations = db.get_cross_section_locations() locations.apply_to_lines(channel_grid.lines, channels) windshieldings = db.get_windshieldings() windshieldings.apply_to_lines(channel_grid.lines) grid += channel_grid - pipes = db.get_pipes() grid += Grid.from_linear_objects( connection_nodes=connection_nodes, objects=pipes, @@ -145,7 +159,6 @@ def _make_gridadmin( connection_node_offset=connection_node_first_id, ) - culverts = db.get_culverts() grid += Grid.from_linear_objects( connection_nodes=connection_nodes, objects=culverts, @@ -159,8 +172,6 @@ def _make_gridadmin( connection_node_offset=connection_node_first_id, ) - weirs = db.get_weirs() - orifices = db.get_orifices() grid += grid.from_structures( connection_nodes=connection_nodes, weirs=weirs, @@ -168,6 +179,7 @@ def _make_gridadmin( line_id_counter=line_id_counter, connection_node_offset=connection_node_first_id, ) + grid.set_cross_sections(cross_section_definitions_unique) grid.set_calculation_types() grid.set_bottom_levels() @@ -179,7 +191,6 @@ def _make_gridadmin( culverts=culverts, ) grid.set_boundary_conditions_1d(db.get_boundary_conditions_1d()) - grid.set_cross_sections(db.get_cross_section_definitions()) grid.set_pumps(db.get_pumps()) if grid.nodes.has_1d and grid.nodes.has_2d: @@ -195,8 +206,8 @@ def _make_gridadmin( culverts=culverts, potential_breaches=potential_breaches, line_id_counter=line_id_counter, + node_open_water_detection=grid_settings.node_open_water_detection, ) - if grid.meta.has_groundwater: channels.apply_has_groundwater_exchange(grid.nodes, grid.lines) pipes.apply_has_groundwater_exchange(grid.nodes, grid.lines) diff --git a/threedigrid_builder/base/nodes.py b/threedigrid_builder/base/nodes.py index 5e822239..e1cb9875 100644 --- a/threedigrid_builder/base/nodes.py +++ b/threedigrid_builder/base/nodes.py @@ -48,7 +48,7 @@ class Node: breach_ids: Tuple[int, int] # referring to 1 or 2 ids of potential breaches has_dem_averaged: int # Boolean attribute to tell if dem is averaged for node. boundary_type: BoundaryType - manhole_id: int # referring to the id of a manhole + is_manhole: bool # node represents a manhole drain_level: float # Drain level entered for manholes initial_waterlevel: float zoom_category: int diff --git a/threedigrid_builder/base/settings.py b/threedigrid_builder/base/settings.py index 95ae3d49..26aaf159 100644 --- a/threedigrid_builder/base/settings.py +++ b/threedigrid_builder/base/settings.py @@ -38,6 +38,7 @@ class GridSettings: grid_space: float dist_calc_points: float kmax: int + node_open_water_detection: int embedded_cutoff_threshold: float = 0.05 max_angle_1d_advection: float = 0.4 * np.pi diff --git a/threedigrid_builder/grid/connection_nodes.py b/threedigrid_builder/grid/connection_nodes.py index 3019fbbf..c126786e 100644 --- a/threedigrid_builder/grid/connection_nodes.py +++ b/threedigrid_builder/grid/connection_nodes.py @@ -24,14 +24,12 @@ class ConnectionNode: the_geom: shapely.Geometry code: str storage_area: float - manhole_id: int + is_manhole: bool calculation_type: CalculationType manhole_indicator: int bottom_level: float drain_level: float surface_level: float - shape: str # enum with classes "00", "01", "02" - width: float initial_waterlevel: float display_name: str zoom_category: int @@ -59,7 +57,7 @@ def get_nodes(self, node_id_counter): - node_type: NODE_1D_STORAGE or NODE_1D_NO_STORAGE depending on storage_area - calculation_type: from calculation_type (which comes from manhole) - dmax: from bottom_level (which comes from manhole) - - manhole_id: id of associated manhole. + - is_manhole: node is a manhole - drain_level: drain_level of associated manhole. - storage_area: area of connection_node. """ @@ -77,7 +75,7 @@ def get_nodes(self, node_id_counter): node_type=node_type, calculation_type=self.calculation_type, dmax=self.bottom_level, - manhole_id=self.manhole_id, + is_manhole=~np.isnan(self.bottom_level), drain_level=self.drain_level, storage_area=self.storage_area, display_name=self.display_name, @@ -93,6 +91,14 @@ def is_closed(self, content_pk): """ return self.storage_area[self.id_to_index(content_pk)] >= 0 + def is_channel(self, content_pk, channels): + """Whether object is connected to a channel""" + has_channel = np.logical_or( + np.isin(self.id, channels.connection_node_start_id), + np.isin(self.id, channels.connection_node_end_id), + ) + return has_channel[content_pk] + def get_1d2d_exchange_levels(self, content_pk, channels, locations, **kwargs): """Compute the exchange level (dpumax) for 1D-2D flowlines. @@ -109,17 +115,16 @@ def get_1d2d_exchange_levels(self, content_pk, channels, locations, **kwargs): """ # get the corresponding connection_node ids and indexes connection_node_idx = self.id_to_index(content_pk) - is_manhole = self.manhole_id[connection_node_idx] != -9999 + is_manhole = ~np.isnan(self.bottom_level[connection_node_idx]) is_manhole_idx = connection_node_idx[is_manhole] - # Check if manhole has drain_level below bottom_level has_lower_drn_lvl = ( self.drain_level[is_manhole_idx] < self.bottom_level[is_manhole_idx] ) if np.any(has_lower_drn_lvl): - ids = self.manhole_id[is_manhole_idx[has_lower_drn_lvl]] + ids = self.id[is_manhole_idx[has_lower_drn_lvl]] logger.warning( - f"Manholes {sorted(ids.tolist())} have a " + f"ConnectionNode {sorted(ids.tolist())} have a " f"bottom_level that is above drain_level." ) @@ -143,7 +148,6 @@ def get_1d2d_exchange_levels(self, content_pk, channels, locations, **kwargs): _put_if_less(dpumax, cn_idx_with_channel, drain_level) # filter out connection nodes that were not in node_idx (non-connected ones) dpumax = dpumax[connection_node_idx] - # for manholes: put in the drain level dpumax[is_manhole] = self.drain_level[is_manhole_idx] return dpumax @@ -270,12 +274,13 @@ def set_bottom_levels(nodes: Nodes, lines: Lines): dmax = np.fmin(nodes.dmax[node_idx], dmax_per_node) # Check if the new node dmax is below the original manhole dmax - is_manhole = nodes.manhole_id[node_idx] != -9999 + is_manhole = nodes.is_manhole[node_idx] + has_lower_dmax = dmax[is_manhole] < nodes.dmax[node_idx[is_manhole]] if np.any(has_lower_dmax): - ids = nodes.manhole_id[node_idx[is_manhole][has_lower_dmax]] + ids = nodes.id[node_idx[is_manhole][has_lower_dmax]] logger.warning( - f"Manholes {sorted(ids.tolist())} have a " + f"Nodes {sorted(ids.tolist())} have a " f"bottom_level that is above one ore more of the following connected " f"objects: channel reference level, pipe/culvert invert level, " f"weir/orifice crest level." diff --git a/threedigrid_builder/grid/cross_section_definitions.py b/threedigrid_builder/grid/cross_section_definitions.py index 2f1566ab..dff7b308 100644 --- a/threedigrid_builder/grid/cross_section_definitions.py +++ b/threedigrid_builder/grid/cross_section_definitions.py @@ -1,5 +1,6 @@ import numpy as np import shapely +from threedi_schema import constants from threedigrid_builder.base import Array from threedigrid_builder.constants import CrossSectionShape @@ -15,16 +16,72 @@ class CrossSectionDefinition: id: int code: str 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 + height: float + width: float + friction_values: str # comma-separated list of floats + vegetation_stem_density: float + vegetation_stem_diameter: float + vegetation_height: float + vegetation_drag_coefficient: float + cross_section_table: str # csv table defining the shape of a tabulated shape + cross_section_vegetation_table: str # csv table with cross section vegetation table + origin_table: str # table definition originates from + origin_id: int # id in origin_table where definition originates from class CrossSectionDefinitions(Array[CrossSectionDefinition]): + def get_unique(self): + """ + Returns a tuple of unique cross section definitions and their mapping to original cross section definitions. + + Returns: + Tuple[CrossSectionDefinitions, Dict[str, Dict[int, int]]]: A tuple where the first element is + a dictionary of unique cross section definitions and the second element is a dictionary mapping + the original tables and rows where these definitions are used. + """ + # Map attributes used to define a cross section for each shape + cross_section_attributes = { + constants.CrossSectionShape.CLOSED_RECTANGLE.value: ["width", "height"], + constants.CrossSectionShape.RECTANGLE.value: ["width"], + constants.CrossSectionShape.CIRCLE.value: ["width"], + constants.CrossSectionShape.EGG.value: ["width"], + constants.CrossSectionShape.TABULATED_RECTANGLE.value: [ + "cross_section_table" + ], + constants.CrossSectionShape.TABULATED_TRAPEZIUM.value: [ + "cross_section_table" + ], + constants.CrossSectionShape.TABULATED_YZ.value: ["cross_section_table"], + constants.CrossSectionShape.INVERTED_EGG.value: ["width"], + } + definition_map = {name: {} for name in np.unique(self.origin_table)} + new_csd_dict = { + key: np.empty(0, dtype=data.dtype) for key, data in vars(self).items() + } + for shape in np.unique(self.shape): + mask = self.shape == shape + # collect cross section definition in array + attr_arr = np.column_stack( + [getattr(self, attr)[mask] for attr in cross_section_attributes[shape]] + ) + # Find unique rows and add these to the new_csd_dict with new id's + u_arr, u_idx = np.unique(attr_arr, return_index=True) + new_id = np.arange(len(u_idx)) + len(new_csd_dict["id"]) + for key in new_csd_dict: + if key == "id": + new_csd_dict[key] = np.concatenate([new_csd_dict[key], new_id]) + else: + new_csd_dict[key] = np.concatenate( + [new_csd_dict[key], getattr(self, key)[mask][u_idx]] + ) + # Map unique cross section definition to table and row of origin + for i, row in enumerate(u_arr): + for idx in np.where(attr_arr == row)[0]: + definition_map[self.origin_table[mask][idx]][ + self.origin_id[mask][idx] + ] = new_id[i] + return CrossSectionDefinitions(**new_csd_dict), definition_map + def convert(self, ids): """Convert to CrossSections. @@ -58,14 +115,24 @@ def convert(self, ids): for i, self_i in enumerate(idx): shape = self.shape[self_i] - tabulator = tabulators[shape] + if shape in [ + CrossSectionShape.TABULATED_YZ, + CrossSectionShape.TABULATED_RECTANGLE, + CrossSectionShape.TABULATED_TRAPEZIUM, + ]: + width, height = _parse_tabulated( + self.cross_section_table[self_i], shape + ) + else: + width = self.width[self_i] + height = self.height[self_i] ( result.shape[i], result.width_1d[i], result.height_1d[i], table, yz, - ) = tabulator(shape, self.width[self_i], self.height[self_i]) + ) = tabulators[shape](shape, width, height) if table is not None: result.count[i] = len(table) tables.append(table) @@ -74,10 +141,11 @@ def convert(self, ids): 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], + self.vegetation_stem_density[self_i], + self.vegetation_stem_diameter[self_i], + self.vegetation_height[self_i], + self.vegetation_drag_coefficient[self_i], + self.cross_section_vegetation_table[self_i], ) tables_yz.append(yz) @@ -118,36 +186,29 @@ class CrossSections(Array[CrossSection]): tables_yz = None -def tabulate_builtin(shape, width, height): +def tabulate_builtin(shape, width: float, height: float): """Tabulate built-in shapes (rectangle, circle) Built-in geometries only require a width to be fully specified. Args: shape (CrossSectionShape): returned as is - width (str): fully specifies the shape - height (str): ignored + width (float): fully specifies the shape + height (float): ignored Returns: tuple: shape, width_1d (float), None, None, None """ - try: - width = float(width) - except ValueError: - raise SchematisationError( - f"Unable to parse cross section definition width (got: '{width}')." - ) - return shape, width, None, None, None -def tabulate_egg(shape, width, height): +def tabulate_egg(shape, width: float, height: float): """Tabulate the egg shape. Args: shape (CrossSectionShape): ignored - width (str): the width of the egg, defines height and increment - height (str): ignored; height is set to 1.5 * width + width (float): the width of the egg, defines height and increment + height (float): ignored; height is set to 1.5 * width Returns: tuple: TABULATED_TRAPEZIUM, width_1d (float), @@ -155,14 +216,6 @@ def tabulate_egg(shape, width, height): """ NUM_INCREMENTS = 16 - # width is the only constant; height derives from width - try: - width = float(width) - except ValueError: - raise SchematisationError( - f"Unable to parse cross section definition width (got: '{width}')." - ) - height = width * 1.5 position = height / 3.0 # some parameter for the 'egg' curve heights = np.linspace(0, height, num=NUM_INCREMENTS, endpoint=True) @@ -190,50 +243,46 @@ def tabulate_inverted_egg(shape, width, height): return type_, width, height, table, None -def tabulate_closed_rectangle(shape, width, height): +def tabulate_closed_rectangle(shape, width: float, height: float): """Tabulate the closed rectangle shape. Args: shape (CrossSectionShape): ignored - width (str): the width of the rectangle - height (str): the height of the rectangle + width (float): the width of the rectangle + height (float): the height of the rectangle Returns: tuple: TABULATED_RECTANGLE, width_1d (float), height (float), table (ndarray of shape (M, 2)), None """ - try: - width = float(width) - height = float(height) - except ValueError: - raise SchematisationError( - f"Unable to parse cross section definition width and/or height " - f"(got: '{width}', '{height}')." - ) table = np.array([[0.0, width], [height, 0.0]], order="F") return CrossSectionShape.TABULATED_RECTANGLE, width, height, table, None -def _parse_tabulated(width, height): +def _parse_tabulated(cross_section_table, shape): try: - heights = np.array([float(x) for x in height.split(" ")]) - widths = np.array([float(x) for x in width.split(" ")]) + left, right = zip( + *[ + [float(item) for item in line.split(",")] + for line in cross_section_table.splitlines() + ] + ) except ValueError: raise SchematisationError( - f"Unable to parse cross section definition width and/or height " - f"(got: '{width}', '{height}')." + f"Unable to parse cross section definition table " + f"(got: '{cross_section_table}')." ) - if len(heights) == 0: + if len(left) == 0: raise SchematisationError( f"Cross section definitions of tabulated or profile type must have at least one " - f"height element (got: {height})." + f"height element (got: {cross_section_table})." ) - if len(heights) != len(widths): - raise SchematisationError( - f"Cross section definitions of tabulated or profile type must have equal number of " - f"height and width elements (got: {height}, {width})." - ) - return widths, heights + if shape == CrossSectionShape.TABULATED_YZ: + # for tabulated_yz, cross section table is y,z + return np.array(left), np.array(right) + else: + # for other tabulated, cross seciton table is height, width + return np.array(right), np.array(left) def tabulate_tabulated(shape, width, height): @@ -248,14 +297,13 @@ def tabulate_tabulated(shape, width, height): tuple: shape, width_1d (float), 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): + if len(height) > 1 and np.any(np.diff(height) < 0.0): raise SchematisationError( f"Cross section definitions of tabulated type must have increasing heights " f"(got: {height})." ) - return shape, np.max(widths), np.max(heights), np.array([heights, widths]).T, None + return shape, np.max(width), np.max(height), np.array([height, width]).T, None def tabulate_yz(shape, width, height): @@ -270,7 +318,8 @@ def tabulate_yz(shape, width, height): tuple: shape, width_1d (float), height_1d (float), table (ndarray of shape (M, 2)), yz (ndarray of shape (M, 4)) """ - ys, zs = _parse_tabulated(width, height) + ys = np.array(width) + zs = np.array(height) is_closed = ys[0] == ys[-1] and zs[0] == zs[-1] if is_closed and len(zs) < 4: raise SchematisationError( @@ -377,23 +426,33 @@ def tabulate_yz(shape, width, height): def set_friction_vegetation_values( yz, friction_values, - vegetation_stem_densities, - vegetation_stem_diameters, - vegetation_heights, - vegetation_drag_coefficients, + vegetation_stem_density, + vegetation_stem_diameter, + vegetation_height, + vegetation_drag_coefficient, + cross_section_vegetation_table, ): """Convert friction and vegetation properties from list into arrays, if available, and add to yz""" if friction_values: - fric = np.array([float(x) for x in friction_values.split(" ")]) + fric = np.array([float(x) for x in friction_values.split(",")]) yz[:-1, 2] = fric - - if vegetation_drag_coefficients: - 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_drag - yz[:-1, 4] = veg_hght - + if cross_section_vegetation_table: + parsed = np.array( + [ + np.fromstring(row, sep=",") + for row in cross_section_vegetation_table.splitlines() + ] + ) + vegetation_stem_density = parsed[:, 0] + vegetation_stem_diameter = parsed[:, 1] + vegetation_height = parsed[:, 2] + vegetation_drag_coefficient = parsed[:, 3] + if vegetation_drag_coefficient is not None: + yz[:-1, 3] = ( + vegetation_stem_density + * vegetation_stem_diameter + * vegetation_drag_coefficient + ) + yz[:-1, 4] = vegetation_height return yz diff --git a/threedigrid_builder/grid/grid.py b/threedigrid_builder/grid/grid.py index 370ceebb..770484de 100644 --- a/threedigrid_builder/grid/grid.py +++ b/threedigrid_builder/grid/grid.py @@ -20,7 +20,7 @@ from threedigrid_builder.base.settings import GridSettings, TablesSettings from threedigrid_builder.constants import ContentType, LineType, NodeType, WKT_VERSION from threedigrid_builder.exceptions import SchematisationError -from threedigrid_builder.grid import zero_d +from threedigrid_builder.grid import ConnectionNodes, zero_d from . import connection_nodes as connection_nodes_module from . import dem_average_area as dem_average_area_module @@ -30,7 +30,7 @@ from .cross_section_definitions import CrossSections from .linear import BaseLinear from .lines_1d2d import Lines1D2D -from .obstacles import Obstacles +from .obstacles import ObstacleAffectsType, Obstacles from .potential_breaches import PotentialBreaches, PotentialBreachPoints osr.UseExceptions() @@ -370,6 +370,10 @@ def from_connection_nodes(cls, connection_nodes, node_id_counter): - nodes.content_pk: the user-supplied id - nodes.node_type: NODE_1D_NO_STORAGE / NODE_1D_STORAGE - nodes.calculation_type: only if set on Manhole + - nodes.dmax: from bottom_level (which comes from manhole) + - nodes.is_manhole: node is a manhole + - nodes.drain_level: drain_level of associated manhole. + - nodes.storage_area: area of connection_node. """ return cls(connection_nodes.get_nodes(node_id_counter), Lines(id=[])) @@ -580,17 +584,43 @@ def set_initial_waterlevels(self, connection_nodes, channels, pipes, culverts): culverts=culverts, ) - def set_obstacles(self, obstacles: Obstacles): + def set_obstacles_2d( + self, + obstacles: Obstacles, + ): """Set obstacles on 2D lines by determining intersection between line_coords (these must be knows at this point) and obstacle geometry. Set kcu to LINE_2D_OBSTACLE and changes flod and flou to crest_level. Args: obstacles (Obstacles) + affects_type (ObstacleAffectsType): which affects attribute of obstacles are considered """ line_2d = [LineType.LINE_2D_U, LineType.LINE_2D_V] selection = np.where(np.isin(self.lines.kcu, line_2d))[0] - crest_level = obstacles.compute_dpumax(self.lines, where=selection)[0] + crest_level = obstacles.compute_dpumax( + self.lines, where=selection, affects_type=ObstacleAffectsType.AFFECTS_2D + )[0] + self.lines.set_2d_crest_levels(crest_level, where=selection) + self.obstacles = obstacles + + def set_obstacles_1d2d( + self, + obstacles: Obstacles, + ): + """Set obstacles on 2D lines by determining intersection between + line_coords (these must be knows at this point) and obstacle geometry. + Set kcu to LINE_2D_OBSTACLE and changes flod and flou to crest_level. + + Args: + obstacles (Obstacles) + affects_type (ObstacleAffectsType): which affects attribute of obstacles are considered + """ + line_2d = [LineType.LINE_2D_U, LineType.LINE_2D_V] + selection = np.where(np.isin(self.lines.kcu, line_2d))[0] + crest_level = obstacles.compute_dpumax( + self.lines, where=selection, affects_type=ObstacleAffectsType.AFFECTS_2D + )[0] self.lines.set_2d_crest_levels(crest_level, where=selection) self.obstacles = obstacles @@ -693,6 +723,7 @@ def add_1d2d_lines( culverts, potential_breaches, line_id_counter, + node_open_water_detection, ) -> Lines1D2D: """Connect 1D and 2D elements by computing 1D-2D lines. @@ -726,10 +757,15 @@ def add_1d2d_lines( for objects in (channels, connection_nodes, pipes, culverts): mask = self.nodes.content_type[node_idx] == objects.content_type content_pk = self.nodes.content_pk[node_idx[mask]] - lines_1d2d.assign_kcu(mask, objects.is_closed(content_pk)) + if node_open_water_detection == 0 and isinstance(objects, ConnectionNodes): + is_closed = objects.is_channel(content_pk, channels) + else: + is_closed = objects.is_closed(content_pk) + lines_1d2d.assign_kcu(mask, is_closed=is_closed) lines_1d2d.assign_dpumax_from_breaches(potential_breaches) lines_1d2d.assign_dpumax_from_exchange_lines(exchange_lines) - lines_1d2d.assign_dpumax_from_obstacles(self.obstacles) + lines_1d2d.assign_dpumax_from_obstacles_open(self.obstacles) + lines_1d2d.assign_dpumax_from_obstacles_closed(self.obstacles) # Go through objects and dispatch to get_1d2d_properties for objects in (channels, connection_nodes, pipes, culverts): mask = self.nodes.content_type[node_idx] == objects.content_type diff --git a/threedigrid_builder/grid/lines_1d2d.py b/threedigrid_builder/grid/lines_1d2d.py index d75cc0ed..1ba4cdba 100644 --- a/threedigrid_builder/grid/lines_1d2d.py +++ b/threedigrid_builder/grid/lines_1d2d.py @@ -1,6 +1,6 @@ import itertools import logging -from typing import Iterator +from typing import Iterator, List import numpy as np import shapely @@ -11,7 +11,7 @@ from .connection_nodes import ConnectionNodes from .exchange_lines import ExchangeLines -from .obstacles import Obstacles +from .obstacles import ObstacleAffectsType, Obstacles from .potential_breaches import PotentialBreaches logger = logging.getLogger(__name__) @@ -340,26 +340,46 @@ def assign_dpumax_from_exchange_lines(self, exchange_lines: ExchangeLines): ], ) - def assign_dpumax_from_obstacles(self, obstacles: Obstacles): - """Set the dpumax based on intersected obstacles - - Only for open connections. Requires line_geometries. - """ - is_open_water = np.isin( - self.kcu, - [ + def assign_dpumax_from_obstacles_open(self, obstacles: Obstacles): + self._assign_dpumax_from_obstacles( + obstacles, + affects_type=ObstacleAffectsType.AFFECTS_1D2D_OPEN_WATER, + line_types=[ LineType.LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER, LineType.LINE_1D2D_DOUBLE_CONNECTED_OPEN_WATER, ], ) - is_open_water_idx = np.where(is_open_water)[0] + + def assign_dpumax_from_obstacles_closed(self, obstacles: Obstacles): + self._assign_dpumax_from_obstacles( + obstacles, + affects_type=ObstacleAffectsType.AFFECTS_1D2D_CLOSED, + line_types=[ + LineType.LINE_1D2D_SINGLE_CONNECTED_CLOSED, + LineType.LINE_1D2D_DOUBLE_CONNECTED_CLOSED, + ], + ) + + def _assign_dpumax_from_obstacles( + self, + obstacles: Obstacles, + affects_type: ObstacleAffectsType, + line_types: List[LineType], + ): + """Set the dpumax based on intersected obstacles + + Only for open connections. Requires line_geometries. + """ + is_included = np.isin(self.kcu, line_types) + is_included_idx = np.where(is_included)[0] obstacle_crest_levels, obstacle_idx = obstacles.compute_dpumax( - self, where=is_open_water_idx + self, where=is_included_idx, affects_type=affects_type ) - self.assign_dpumax(is_open_water, obstacle_crest_levels) + self.assign_dpumax(is_included, obstacle_crest_levels) mask = obstacle_idx != -9999 + self.assign_ds1d_half_from_obstacles( - obstacles, is_open_water_idx[mask], obstacle_idx[mask] + obstacles, is_included_idx[mask], obstacle_idx[mask] ) def assign_ds1d_half_from_obstacles( diff --git a/threedigrid_builder/grid/obstacles.py b/threedigrid_builder/grid/obstacles.py index a8e8bb2b..c7382533 100644 --- a/threedigrid_builder/grid/obstacles.py +++ b/threedigrid_builder/grid/obstacles.py @@ -1,22 +1,36 @@ +from enum import Enum + import numpy as np import shapely from threedigrid_builder.base import Array, Lines +class ObstacleAffectsType(Enum): + AFFECTS_2D = "affects_2d" + AFFECTS_1D2D_OPEN_WATER = "affects_1d2d_open_water" + AFFECTS_1D2D_CLOSED = "affects_1d2d_closed" + + class Obstacle: id: int crest_level: float the_geom: shapely.Geometry + affects_1d2d_closed: bool + affects_1d2d_open_water: bool + affects_2d: bool class Obstacles(Array[Obstacle]): - def compute_dpumax(self, lines: Lines, where): + def compute_dpumax( + self, lines: Lines, where, affects_type: ObstacleAffectsType + ) -> Lines: """Compute the dpumax for lines that intersect with the obstacles. Args: lines (Lines) where (int array): indices into lines to compute dpumax for + affects_attr (str): name of Obstacles attribute used Returns: dpumax for each line (length equals length of 'where' array) @@ -33,6 +47,8 @@ def compute_dpumax(self, lines: Lines, where): lines_tree = shapely.STRtree(lines.line_geometries[where]) inscts = lines_tree.query(self.the_geom, predicate="intersects") for i in range(len(self)): + if not getattr(self, affects_type.value)[i]: + continue line_idx = inscts[1, np.where(inscts[0, :] == i)[0]] is_greater = ~(self.crest_level[i] <= exchange_level[line_idx]) exchange_level[line_idx[is_greater]] = self.crest_level[i] diff --git a/threedigrid_builder/interface/db.py b/threedigrid_builder/interface/db.py index 118072e0..d68c3c3f 100644 --- a/threedigrid_builder/interface/db.py +++ b/threedigrid_builder/interface/db.py @@ -4,14 +4,14 @@ from contextlib import contextmanager from enum import Enum from functools import lru_cache -from typing import Callable, ContextManager, Dict, Optional +from typing import Callable, ContextManager, Dict, List, Optional, Union import numpy as np import shapely from condenser import NumpyQuery from pyproj import Transformer from pyproj.crs import CRS -from sqlalchemy import cast, func, inspect, Integer +from sqlalchemy import case, cast, func, inspect, Integer, literal from sqlalchemy.orm import Session from threedi_schema import custom_types, models, ModelSchema, ThreediDatabase @@ -43,7 +43,7 @@ # hardcoded source projection SOURCE_EPSG = 4326 -MIN_SQLITE_VERSION = 227 +MIN_SQLITE_VERSION = 228 DAY_IN_SECONDS = 24.0 * 3600.0 @@ -108,6 +108,44 @@ def arr_to_attr_dict( } +def map_cross_section_definition( + objects: List[Union[CrossSectionLocations, Pipes, Weirs, Orifices, Culverts]], + definition_map: Dict[str, Dict[int, int]], +) -> None: + """ + Set cross section definition ids for cross_section_locations, + pipes, weirs, orifices and culverts to match the unique + cross section locations. + + Args: + objects: List of objects (CrossSectionLocations, Pipes, Weirs, Orifices, Culverts) to map the definition to + definition_map: Mapping of object names to their definition IDs + """ + # Map shapes to key in definition_map + object_map = { + CrossSectionLocations: "cross_section_location", + Pipes: "pipe", + Weirs: "weir", + Orifices: "orifice", + Culverts: "culvert", + } + for object in objects: + object_name = object_map.get(type(object)) + if object_name is None: + raise ValueError(f"Object of type {type(object)} cannot be mapped") + mapping = definition_map.get(object_name) + if mapping is None: + continue + idx = object.id_to_index(list(mapping.keys())) + # set correct cross section definition id's for mapped object + if isinstance(object, CrossSectionLocations): + object.definition_id[idx] = np.array(list(mapping.values()), dtype=int) + else: + object.cross_section_definition_id[idx] = np.array( + list(mapping.values()), dtype=int + ) + + class SQLite: def __init__(self, path: pathlib.Path, upgrade=False, convert_to_geopackage=False): if not path.exists(): @@ -129,9 +167,7 @@ def get_version(self) -> int: def upgrade(self, convert_to_geopackage=False): schema = ModelSchema(self.db) - schema.upgrade( - backup=False, set_views=False, convert_to_geopackage=convert_to_geopackage - ) + schema.upgrade(backup=False, convert_to_geopackage=convert_to_geopackage) @contextmanager def get_session(self) -> ContextManager[Session]: @@ -328,9 +364,9 @@ def get_surfaces(self) -> Surfaces: models.SurfaceParameter.min_infiltration_capacity, models.SurfaceParameter.infiltration_decay_constant, models.SurfaceParameter.infiltration_recovery_constant, - models.ConnectionNode.id.label("connection_node_id"), - models.ConnectionNode.the_geom.label("connection_node_the_geom"), + models.SurfaceMap.connection_node_id, models.SurfaceMap.percentage, + models.ConnectionNode.geom.label("connection_node_the_geom"), ) .select_from(models.Surface) .join(models.SurfaceParameter) @@ -347,9 +383,6 @@ def get_surfaces(self) -> Surfaces: # reproject arr["geom"] = self.reproject(arr["geom"]) - arr["connection_node_the_geom"] = self.reproject( - arr["connection_node_the_geom"] - ) return Surfaces( id=np.arange(0, len(arr["surface_id"] + 1), dtype=int), @@ -392,15 +425,14 @@ def get_boundary_conditions_2d(self) -> BoundaryConditions2D: def get_channels(self) -> Channels: """Return Channels""" cols = [ - models.Channel.the_geom, - models.Channel.dist_calc_points, + models.Channel.geom, + models.Channel.calculation_point_distance, models.Channel.id, models.Channel.code, - models.Channel.connection_node_start_id, - models.Channel.connection_node_end_id, - models.Channel.calculation_type, + models.Channel.connection_node_id_start, + models.Channel.connection_node_id_end, + models.Channel.exchange_type, models.Channel.display_name, - models.Channel.zoom_category, models.Channel.exchange_thickness, models.Channel.hydraulic_conductivity_out, models.Channel.hydraulic_conductivity_in, @@ -409,84 +441,132 @@ def get_channels(self) -> Channels: with self.get_session() as session: arr = session.query(*cols).order_by(models.Channel.id).as_structarray() - arr["the_geom"] = self.reproject(arr["the_geom"]) + arr["geom"] = self.reproject(arr["geom"]) # map "old" calculation types (100, 101, 102, 105) to (0, 1, 2, 5) - arr["calculation_type"][arr["calculation_type"] >= 100] -= 100 + arr["exchange_type"][arr["exchange_type"] >= 100] -= 100 arr["hydraulic_conductivity_out"] /= DAY_IN_SECONDS arr["hydraulic_conductivity_in"] /= DAY_IN_SECONDS + attr_dict = arr_to_attr_dict( + arr, + { + "geom": "the_geom", + "exchange_type": "calculation_type", + "calculation_point_distance": "dist_calc_points", + "connection_node_id_start": "connection_node_start_id", + "connection_node_id_end": "connection_node_end_id", + }, + ) + # transform to a Channels object - return Channels(**{name: arr[name] for name in arr.dtype.names}) + return Channels(**attr_dict) def get_connection_nodes(self) -> ConnectionNodes: - """Return ConnectionNodes (which are enriched using the manhole table)""" + """Return ConnectionNodes""" cols = [ - models.ConnectionNode.the_geom, + models.ConnectionNode.geom, models.ConnectionNode.id, models.ConnectionNode.code, models.ConnectionNode.storage_area, - models.ConnectionNode.initial_waterlevel, - models.Manhole.id.label("manhole_id"), - models.Manhole.calculation_type, - models.Manhole.bottom_level, - models.Manhole.drain_level, - models.Manhole.manhole_indicator, - models.Manhole.surface_level, - models.Manhole.shape, - models.Manhole.width, - models.Manhole.display_name, - models.Manhole.zoom_category, - models.Manhole.exchange_thickness, - models.Manhole.hydraulic_conductivity_out, - models.Manhole.hydraulic_conductivity_in, + models.ConnectionNode.initial_water_level, + models.ConnectionNode.exchange_type, + models.ConnectionNode.bottom_level, + models.ConnectionNode.exchange_level, + models.ConnectionNode.visualisation, + models.ConnectionNode.manhole_surface_level, + models.ConnectionNode.display_name, + models.ConnectionNode.exchange_thickness, + models.ConnectionNode.hydraulic_conductivity_out, + models.ConnectionNode.hydraulic_conductivity_in, ] with self.get_session() as session: arr = ( - session.query(*cols) - .join(models.ConnectionNode.manholes, isouter=True) - .order_by(models.ConnectionNode.id) - .as_structarray() + session.query(*cols).order_by(models.ConnectionNode.id).as_structarray() ) - arr["the_geom"] = self.reproject(arr["the_geom"]) + arr["geom"] = self.reproject(arr["geom"]) - # replace -9999.0 with NaN in initial_waterlevel - arr["initial_waterlevel"][arr["initial_waterlevel"] == -9999.0] = np.nan + # replace -9999.0 with NaN in initial_water_level + arr["initial_water_level"][arr["initial_water_level"] == -9999.0] = np.nan arr["hydraulic_conductivity_out"] /= DAY_IN_SECONDS arr["hydraulic_conductivity_in"] /= DAY_IN_SECONDS - return ConnectionNodes(**{name: arr[name] for name in arr.dtype.names}) + attr_dict = arr_to_attr_dict( + arr, + { + "geom": "the_geom", + "initial_water_level": "initial_waterlevel", + "exchange_type": "calculation_type", + "exchange_level": "drain_level", + "visualisation": "manhole_indicator", + "manhole_surface_level": "surface_level", + }, + ) - def get_cross_section_definitions(self) -> CrossSectionDefinitions: - """Return CrossSectionDefinitions""" - with self.get_session() as session: - arr = ( - session.query( - models.CrossSectionDefinition.id, - models.CrossSectionDefinition.code, - 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() - ) + return ConnectionNodes(**attr_dict) - # map shape 10 to 1 (circle) to match CrossSectionShape enum - arr["shape"][arr["shape"] == 10] = 1 - # map shape 11 to 5 (tabulated rectangle) to match CrossSectionShape enum - arr["shape"][arr["shape"] == 11] = 5 - # map shape 12 to 6 (tabulated trapezium) to match CrossSectionShape enum - arr["shape"][arr["shape"] == 12] = 6 + def get_cross_section_definition_for_table(self, table) -> np.ndarray: + with self.get_session() as session: + cols = [ + literal(table.__tablename__).label("origin_table"), + table.id.label("origin_id"), + table.cross_section_shape.label("shape"), + table.cross_section_width.label("width"), + table.cross_section_height.label("height"), + table.cross_section_table, + ] + if table == models.CrossSectionLocation: + cols += [ + table.cross_section_friction_values.label("friction_values"), + table.cross_section_vegetation_table, + table.vegetation_stem_density, + table.vegetation_stem_diameter, + table.vegetation_height, + table.vegetation_drag_coefficient, + ] + arr = session.query(*cols).select_from(table).as_structarray() + # map shape 10 to 1 (circle) to match CrossSectionShape enum + arr["shape"][arr["shape"] == 10] = 1 + # map shape 11 to 5 (tabulated rectangle) to match CrossSectionShape enum + arr["shape"][arr["shape"] == 11] = 5 + # map shape 12 to 6 (tabulated trapezium) to match CrossSectionShape enum + arr["shape"][arr["shape"] == 12] = 6 + return arr - # transform to a CrossSectionDefinitions object - return CrossSectionDefinitions(**{name: arr[name] for name in arr.dtype.names}) + def get_cross_section_definitions(self) -> CrossSectionDefinitions: + """Return CrossSectionDefinitions""" + attr_dict = { + "id": np.empty(0, dtype="i4"), + "origin_table": np.empty(0, dtype="O"), + "origin_id": np.empty(0, dtype="i4"), + "shape": np.empty(0, dtype="O"), + "width": np.empty(0, dtype="f8"), + "height": np.empty(0, dtype="f8"), + "cross_section_table": np.empty(0, dtype="O"), + "friction_values": np.empty(0, dtype="O"), + "cross_section_vegetation_table": np.empty(0, dtype="O"), + } + for table in [ + models.CrossSectionLocation, + models.Culvert, + models.Orifice, + models.Pipe, + models.Weir, + ]: + arr = self.get_cross_section_definition_for_table(table) + if len(arr) == 0: + continue + for name in attr_dict.keys(): + if name == "id": + data = len(attr_dict["id"]) + np.arange(len(arr)) + elif name in arr.dtype.names: + data = arr[name] + else: + data = np.empty(len(arr), attr_dict[name].dtype) + attr_dict[name] = np.concatenate((attr_dict[name], data)) + + return CrossSectionDefinitions(**attr_dict) def get_cross_section_locations(self) -> CrossSectionLocations: """Return CrossSectionLocations""" @@ -495,8 +575,7 @@ def get_cross_section_locations(self) -> CrossSectionLocations: session.query( models.CrossSectionLocation.id, models.CrossSectionLocation.code, - models.CrossSectionLocation.the_geom, - models.CrossSectionLocation.definition_id, + models.CrossSectionLocation.geom, models.CrossSectionLocation.channel_id, models.CrossSectionLocation.reference_level, models.CrossSectionLocation.bank_level, @@ -510,52 +589,83 @@ def get_cross_section_locations(self) -> CrossSectionLocations: .order_by(models.CrossSectionLocation.id) .as_structarray() ) + arr["geom"] = self.reproject(arr["geom"]) - arr["the_geom"] = self.reproject(arr["the_geom"]) + attr_dict = arr_to_attr_dict(arr, {"geom": "the_geom"}) # transform to a CrossSectionLocations object - return CrossSectionLocations(**{name: arr[name] for name in arr.dtype.names}) + return CrossSectionLocations(**attr_dict) def get_culverts(self) -> Culverts: """Return Culverts""" + cols = [ + models.Culvert.id, + models.Culvert.code, + models.Culvert.geom, + models.Culvert.calculation_point_distance, + models.Culvert.connection_node_id_start, + models.Culvert.connection_node_id_end, + models.Culvert.exchange_type, + models.Culvert.invert_level_start, + models.Culvert.invert_level_end, + models.Culvert.discharge_coefficient_negative, + models.Culvert.discharge_coefficient_positive, + models.Culvert.display_name, + case( + { + models.Culvert.friction_value.isnot(None) + & models.Culvert.friction_type.isnot( + None + ): models.Culvert.friction_value + }, + else_=models.Material.friction_coefficient, + ).label("friction_value"), + case( + { + models.Culvert.friction_value.isnot(None) + & models.Culvert.friction_type.isnot( + None + ): models.Culvert.friction_type + }, + else_=models.Material.friction_type, + ).label("friction_type"), + ] + with self.get_session() as session: arr = ( - session.query( - models.Culvert.id, - models.Culvert.code, - models.Culvert.the_geom, - models.Culvert.dist_calc_points, - models.Culvert.connection_node_start_id, - models.Culvert.connection_node_end_id, - models.Culvert.calculation_type, - models.Culvert.cross_section_definition_id, - models.Culvert.invert_level_start_point, - models.Culvert.invert_level_end_point, - models.Culvert.discharge_coefficient_negative, - models.Culvert.discharge_coefficient_positive, - models.Culvert.friction_type, - models.Culvert.friction_value, - models.Culvert.display_name, - models.Culvert.zoom_category, + session.query(*cols) + .outerjoin( + models.Material, models.Culvert.material_id == models.Material.id ) .order_by(models.Culvert.id) .as_structarray() ) - arr["the_geom"] = self.reproject(arr["the_geom"]) + arr["geom"] = self.reproject(arr["geom"]) # map friction_type 4 to friction_type 2 to match crosssectionlocation enum arr["friction_type"][arr["friction_type"] == 4] = 2 # When no calculation type is provides we default to isolated - arr["calculation_type"][ - arr["calculation_type"] == -9999 - ] = LineType.LINE_1D_ISOLATED + arr["exchange_type"][arr["exchange_type"] == -9999] = LineType.LINE_1D_ISOLATED # map "old" calculation types (100, 101, 102, 105) to (0, 1, 2, 5) - arr["calculation_type"][arr["calculation_type"] >= 100] -= 100 + arr["exchange_type"][arr["exchange_type"] >= 100] -= 100 + + attr_dict = arr_to_attr_dict( + arr, + { + "exchange_type": "calculation_type", + "calculation_point_distance": "dist_calc_points", + "invert_level_start": "invert_level_start_point", + "invert_level_end": "invert_level_end_point", + "geom": "the_geom", + "connection_node_id_start": "connection_node_start_id", + "connection_node_id_end": "connection_node_end_id", + }, + ) # transform to a CrossSectionLocations object - return Culverts(**{name: arr[name] for name in arr.dtype.names}) + return Culverts(**attr_dict) def get_exchange_lines(self) -> ExchangeLines: with self.get_session() as session: @@ -641,6 +751,9 @@ def get_obstacles(self) -> Obstacles: models.Obstacle.geom, models.Obstacle.id, models.Obstacle.crest_level, + models.Obstacle.affects_2d, + models.Obstacle.affects_1d2d_closed, + models.Obstacle.affects_1d2d_open_water, ) .order_by(models.Obstacle.id) .as_structarray() @@ -654,32 +767,58 @@ def get_obstacles(self) -> Obstacles: def get_orifices(self) -> Orifices: """Return Orifices""" + cols = [ + models.Orifice.id, + models.Orifice.code, + models.Orifice.connection_node_id_start, + models.Orifice.connection_node_id_end, + models.Orifice.crest_level, + models.Orifice.crest_type, + models.Orifice.discharge_coefficient_negative, + models.Orifice.discharge_coefficient_positive, + models.Orifice.display_name, + models.Orifice.sewerage, + case( + { + models.Orifice.friction_value.isnot(None) + & models.Orifice.friction_type.isnot( + None + ): models.Orifice.friction_value + }, + else_=models.Material.friction_coefficient, + ).label("friction_value"), + case( + { + models.Orifice.friction_value.isnot(None) + & models.Orifice.friction_type.isnot( + None + ): models.Orifice.friction_type + }, + else_=models.Material.friction_type, + ).label("friction_type"), + ] with self.get_session() as session: arr = ( - session.query( - models.Orifice.id, - models.Orifice.code, - models.Orifice.connection_node_start_id, - models.Orifice.connection_node_end_id, - models.Orifice.crest_level, - models.Orifice.crest_type, - models.Orifice.cross_section_definition_id, - models.Orifice.discharge_coefficient_negative, - models.Orifice.discharge_coefficient_positive, - models.Orifice.friction_type, - models.Orifice.friction_value, - models.Orifice.display_name, - models.Orifice.zoom_category, - models.Orifice.sewerage, + session.query(*cols) + .outerjoin( + models.Material, models.Orifice.material_id == models.Material.id ) .order_by(models.Orifice.id) .as_structarray() ) - # map friction_type 4 to friction_type 2 to match crosssectionlocation enum arr["friction_type"][arr["friction_type"] == 4] = 2 - - return Orifices(**{name: arr[name] for name in arr.dtype.names}) + attr_dict = arr_to_attr_dict( + arr, + { + "exchange_type": "calculation_type", + "calculation_point_distance": "dist_calc_points", + "connection_node_id_start": "connection_node_start_id", + "connection_node_id_end": "connection_node_end_id", + "geom": "the_geom", + }, + ) + return Orifices(**attr_dict) def get_pipes(self) -> Pipes: """Return Pipes""" @@ -687,88 +826,149 @@ def get_pipes(self) -> Pipes: models.Pipe.id, models.Pipe.code, models.Pipe.sewerage_type, - models.Pipe.calculation_type, - models.Pipe.invert_level_start_point, - models.Pipe.invert_level_end_point, - models.Pipe.friction_type, - models.Pipe.friction_value, - models.Pipe.dist_calc_points, - models.Pipe.connection_node_start_id, - models.Pipe.connection_node_end_id, - models.Pipe.cross_section_definition_id, + models.Pipe.exchange_type, + models.Pipe.invert_level_start, + models.Pipe.invert_level_end, + models.Pipe.calculation_point_distance, + models.Pipe.connection_node_id_start, + models.Pipe.connection_node_id_end, models.Pipe.display_name, - models.Pipe.zoom_category, - models.Pipe.material, models.Pipe.exchange_thickness, models.Pipe.hydraulic_conductivity_out, models.Pipe.hydraulic_conductivity_in, + models.Pipe.material_id, + case( + { + models.Pipe.friction_value.isnot(None) + & models.Pipe.friction_type.isnot(None): models.Pipe.friction_value + }, + else_=models.Material.friction_coefficient, + ).label("friction_value"), + case( + { + models.Pipe.friction_value.isnot(None) + & models.Pipe.friction_type.isnot(None): models.Pipe.friction_type + }, + else_=models.Material.friction_type, + ).label("friction_type"), ] - with self.get_session() as session: - arr = session.query(*cols).order_by(models.Pipe.id).as_structarray() + arr = ( + session.query(*cols) + .outerjoin( + models.Material, models.Pipe.material_id == models.Material.id + ) + .order_by(models.Pipe.id) + .as_structarray() + ) # map friction_type 4 to friction_type 2 to match crosssectionlocation enum arr["friction_type"][arr["friction_type"] == 4] = 2 arr["hydraulic_conductivity_out"] /= DAY_IN_SECONDS arr["hydraulic_conductivity_in"] /= DAY_IN_SECONDS + attr_dict = arr_to_attr_dict( + arr, + { + "exchange_type": "calculation_type", + "calculation_point_distance": "dist_calc_points", + "material_id": "material", + "invert_level_start": "invert_level_start_point", + "invert_level_end": "invert_level_end_point", + "connection_node_id_start": "connection_node_start_id", + "connection_node_id_end": "connection_node_end_id", + "geom": "the_geom", + }, + ) + # transform to a Pipes object - return Pipes(**{name: arr[name] for name in arr.dtype.names}) + return Pipes(**attr_dict) def get_pumps(self) -> Pumps: with self.get_session() as session: arr = ( session.query( - models.Pumpstation.id, - models.Pumpstation.code, - models.Pumpstation.capacity, - models.Pumpstation.connection_node_start_id, - models.Pumpstation.connection_node_end_id, - models.Pumpstation.type_, - models.Pumpstation.start_level, - models.Pumpstation.lower_stop_level, - models.Pumpstation.upper_stop_level, - models.Pumpstation.display_name, - models.Pumpstation.zoom_category, + models.Pump.id, + models.Pump.code, + models.Pump.capacity, + models.Pump.connection_node_id, + models.PumpMap.connection_node_id_end, + models.Pump.type_, + models.Pump.start_level, + models.Pump.lower_stop_level, + models.Pump.upper_stop_level, + models.Pump.display_name, ) - .order_by(models.Pumpstation.id) + .outerjoin(models.PumpMap, models.Pump.id == models.PumpMap.pump_id) + .order_by(models.Pump.id) .as_structarray() ) - # Pump capicity is entered as L/s but we need m3/s. + # Pump capacity is entered as L/s but we need m3/s. arr["capacity"] = arr["capacity"] / 1000 + attr_dict = arr_to_attr_dict( + arr, + { + "connection_node_id": "connection_node_start_id", + "connection_node_id_end": "connection_node_end_id", + "geom": "the_geom", + }, + ) + # transform to a Pumps object - return Pumps(**{name: arr[name] for name in arr.dtype.names}) + return Pumps(**attr_dict) def get_weirs(self) -> Weirs: """Return Weirs""" + cols = [ + models.Weir.id, + models.Weir.code, + models.Weir.connection_node_id_start, + models.Weir.connection_node_id_end, + models.Weir.crest_level, + models.Weir.crest_type, + models.Weir.discharge_coefficient_negative, + models.Weir.discharge_coefficient_positive, + models.Weir.display_name, + models.Weir.sewerage, + case( + { + models.Weir.friction_value.isnot(None) + & models.Weir.friction_type.isnot(None): models.Weir.friction_value + }, + else_=models.Material.friction_coefficient, + ).label("friction_value"), + case( + { + models.Weir.friction_value.isnot(None) + & models.Weir.friction_type.isnot(None): models.Weir.friction_type + }, + else_=models.Material.friction_type, + ).label("friction_type"), + ] with self.get_session() as session: arr = ( - session.query( - models.Weir.id, - models.Weir.code, - models.Weir.connection_node_start_id, - models.Weir.connection_node_end_id, - models.Weir.crest_level, - models.Weir.crest_type, - models.Weir.cross_section_definition_id, - models.Weir.discharge_coefficient_negative, - models.Weir.discharge_coefficient_positive, - models.Weir.friction_type, - models.Weir.friction_value, - models.Weir.display_name, - models.Weir.zoom_category, - models.Weir.sewerage, + session.query(*cols) + .outerjoin( + models.Material, models.Weir.material_id == models.Material.id ) .order_by(models.Weir.id) .as_structarray() ) - # map friction_type 4 to friction_type 2 to match crosssectionlocation enum arr["friction_type"][arr["friction_type"] == 4] = 2 - - return Weirs(**{name: arr[name] for name in arr.dtype.names}) + attr_dict = arr_to_attr_dict( + arr, + { + "exchange_type": "calculation_type", + "calculation_point_distance": "dist_calc_points", + "connection_node_id_start": "connection_node_start_id", + "connection_node_id_end": "connection_node_end_id", + "geom": "the_geom", + }, + ) + return Weirs(**attr_dict) def get_windshieldings(self) -> Windshieldings: with self.get_session() as session: diff --git a/threedigrid_builder/interface/dict_out.py b/threedigrid_builder/interface/dict_out.py index a98361f1..fef28cde 100644 --- a/threedigrid_builder/interface/dict_out.py +++ b/threedigrid_builder/interface/dict_out.py @@ -28,7 +28,7 @@ "storage_area", "boundary_id", # referring to the id of the boundary condition "boundary_type", - "manhole_id", # referring to the id of a manhole + "is_manhole", "initial_waterlevel", ) diff --git a/threedigrid_builder/interface/gridadmin.py b/threedigrid_builder/interface/gridadmin.py index be695a0c..139b955f 100644 --- a/threedigrid_builder/interface/gridadmin.py +++ b/threedigrid_builder/interface/gridadmin.py @@ -308,7 +308,6 @@ def write_nodes(self, nodes, group_name="nodes"): is_2d = np.isin( nodes.node_type, [NodeType.NODE_2D_OPEN_WATER, NodeType.NODE_2D_GROUNDWATER] ) - is_manhole = nodes.manhole_id != -9999 # Datasets that match directly to a nodes attribute: self.write_dataset(group, "id", nodes.id + 1) @@ -335,11 +334,11 @@ def write_nodes(self, nodes, group_name="nodes"): self.write_dataset( group, "is_manhole", - np.where(is_manhole, 1, -9999).astype("i4"), + np.where(nodes.is_manhole, 1, -9999).astype("i4"), ) self.write_dataset(group, "dimp", nodes.dimp) self.write_dataset( - group, "bottom_level", np.where(is_manhole, nodes.dmax, np.nan) + group, "bottom_level", np.where(nodes.is_manhole, nodes.dmax, np.nan) ) self.write_dataset(group, "z_coordinate", nodes.dmax) @@ -358,7 +357,10 @@ def write_nodes(self, nodes, group_name="nodes"): group, "display_name", to_bytes_array(nodes.display_name, 64) ) self.write_dataset(group, "zoom_category", nodes.zoom_category) - self.write_dataset(group, "manhole_id", nodes.manhole_id) + # Set manhole_id to match nodes.id when there is a manhole + self.write_dataset( + group, "manhole_id", np.where(nodes.is_manhole, nodes.id + 1, -9999) + ) self.write_dataset(group, "manhole_indicator", nodes.manhole_indicator) self.write_dataset(group, "shape", to_bytes_array(nodes.shape, 4)) self.write_dataset(group, "drain_level", nodes.drain_level) diff --git a/threedigrid_builder/tests/conftest.py b/threedigrid_builder/tests/conftest.py index 0946ab0c..1807b30e 100644 --- a/threedigrid_builder/tests/conftest.py +++ b/threedigrid_builder/tests/conftest.py @@ -99,6 +99,7 @@ def grid_all(): grid_space=20.0, dist_calc_points=25.0, kmax=4, + node_open_water_detection=1, ), tables_settings=TablesSettings( table_step_size=0.05, @@ -137,6 +138,9 @@ def grid_all(): shapely.linestrings([[1, 1], [2, 2], [4, 4]]), shapely.linestrings([[1, 1], [2, 2], [3, 3]]), ], + affects_2d=[True, True], + affects_1d2d_open_water=[True, True], + affects_1d2d_closed=[True, True], ) breaches = PotentialBreaches( id=[0, 1], diff --git a/threedigrid_builder/tests/test_connection_nodes.py b/threedigrid_builder/tests/test_connection_nodes.py index d4fde138..aaecff00 100644 --- a/threedigrid_builder/tests/test_connection_nodes.py +++ b/threedigrid_builder/tests/test_connection_nodes.py @@ -35,9 +35,8 @@ def connection_nodes(): id=np.array([1, 3, 4, 9, 10]), code=np.array(["one", "two", "", "", ""]), storage_area=np.array([15, np.nan, np.nan, np.nan, 0]), - manhole_id=np.array([42, 11, -9999, -9999, -9999]), calculation_type=np.array([1, 2, 5, 2, 2]), - bottom_level=np.array([2.1, np.nan, np.nan, np.nan, 2.1]), + bottom_level=np.array([2.1, 2.1, np.nan, np.nan, 2.1]), drain_level=[1.2, np.nan, np.nan, np.nan, np.nan], ) @@ -200,7 +199,7 @@ def test_bottom_levels_above_invert_level(structure_type, caplog): content_type=ContentType.TYPE_V2_CONNECTION_NODES, dmax=[10.0, 20.0], content_pk=[52, 22], - manhole_id=[3, 6], + is_manhole=[True, True], ) lines = Lines( id=[1], @@ -213,7 +212,7 @@ def test_bottom_levels_above_invert_level(structure_type, caplog): # setting the bottom levels emits a warning set_bottom_levels(nodes, lines) - assert caplog.messages[0].startswith("Manholes [3] have a bottom_level") + assert caplog.messages[0].startswith("Nodes [1] have a bottom_level") # assert the resulting value of dmax assert_almost_equal(nodes.dmax, [3.0, 20.0]) @@ -267,10 +266,9 @@ def test_get_1d2d_exchange_levels(connection_nodes: ConnectionNodes, caplog): ) assert_array_equal(actual, [1.2, np.nan, np.nan, 1.0]) - assert ( caplog.record_tuples[0][2] - == "Manholes [42] have a bottom_level that is above drain_level." + == "ConnectionNode [1] have a bottom_level that is above drain_level." ) @@ -279,6 +277,23 @@ def test_is_closed(connection_nodes): assert_array_equal(actual, [True, False, False, False, True]) +@pytest.mark.parametrize( + "start_id, end_id, content_pk, expected", + [ + ([1, 3], [9, 10], np.arange(5, dtype=int), [True, True, False, True, True]), + ([1, 3], [9, 10], [0, 1], [True, True]), + ], +) +def test_has_channel(connection_nodes, start_id, end_id, content_pk, expected): + channels = Channels( + id=np.arange(len(start_id), dtype=int), + connection_node_start_id=start_id, + connection_node_end_id=end_id, + ) + has_channel = connection_nodes.is_channel(content_pk, channels) + np.testing.assert_array_equal(has_channel, expected) + + @pytest.mark.parametrize( "area,thickness,hc_out,hc_in,expected", [ diff --git a/threedigrid_builder/tests/test_cross_section_definitions.py b/threedigrid_builder/tests/test_cross_section_definitions.py index a1ac28fe..80abf0ae 100644 --- a/threedigrid_builder/tests/test_cross_section_definitions.py +++ b/threedigrid_builder/tests/test_cross_section_definitions.py @@ -3,11 +3,13 @@ import numpy as np import pytest from numpy.testing import assert_almost_equal, assert_array_equal +from threedi_schema import constants from threedigrid_builder.constants import CrossSectionShape from threedigrid_builder.exceptions import SchematisationError from threedigrid_builder.grid import CrossSectionDefinitions, CrossSections from threedigrid_builder.grid.cross_section_definitions import ( + _parse_tabulated, set_friction_vegetation_values, tabulate_builtin, tabulate_closed_rectangle, @@ -26,8 +28,9 @@ def cross_section_definitions(): id=[1, 3, 9], code=["1", "3", "9"], shape=[SHP.CIRCLE, SHP.TABULATED_TRAPEZIUM, SHP.TABULATED_RECTANGLE], - width=["1.22", "3.7 5.0", "5.0 6.0"], - height=["", "1 2", "1 2"], + width=[1.22, None, None], + height=[None, None, None], + cross_section_table=[None, "1,3.7\n2,5.0", "1,5\n2,6"], ) @@ -89,13 +92,13 @@ def test_convert_empty(cross_section_definitions): def test_tabulate_builtin(): - actual = tabulate_builtin("my-shape", "1.52", "1.33") + actual = tabulate_builtin("my-shape", 1.52, 1.33) assert actual == ("my-shape", 1.52, None, None, None) def test_tabulate_closed_rectangle(): shape, width_1d, height_1d, table, yz = tabulate_closed_rectangle( - "my-shape", "1.52", "5.2" + "my-shape", 1.52, 5.2 ) assert shape == CrossSectionShape.TABULATED_RECTANGLE @@ -106,7 +109,7 @@ def test_tabulate_closed_rectangle(): def test_tabulate_egg(): - shape, width_1d, height_1d, table, yz = tabulate_egg("my-shape", "1.52", "ignored") + shape, width_1d, height_1d, table, yz = tabulate_egg("my-shape", 1.52, 0) assert shape == CrossSectionShape.TABULATED_TRAPEZIUM assert width_1d == 1.52 @@ -138,10 +141,12 @@ def test_tabulate_egg(): def test_tabulate_tabulated(): + width, height = _parse_tabulated( + "0,1\n1,2\n2,3", CrossSectionShape.TABULATED_RECTANGLE + ) shape, width_1d, height_1d, table, yz = tabulate_tabulated( - "my-shape", "1 2 3", "0 1 2" + "my-shape", width, height ) - assert shape == "my-shape" assert width_1d == 3.0 # the max assert height_1d == 2.0 @@ -149,25 +154,16 @@ def test_tabulate_tabulated(): assert yz is None -@pytest.mark.parametrize( - "width,height,match", - [ - ("", "1", r"Unable to parse cross section definition.*"), - ("1", "", r"Unable to parse cross section definition.*"), - ("", "", r"Unable to parse cross section definition.*"), - ("1", "1 2", r".*of tabulated or profile type must have equal number.*"), - ("1 2", "1", r".*of tabulated or profile type must have equal number.*"), - ("1 1", "2 1", r".*of tabulated type must have increasing heights.*"), - ], -) -def test_tabulate_tabulated_err(width, height, match): - with pytest.raises(SchematisationError, match=match): - tabulate_tabulated(CrossSectionShape.TABULATED_RECTANGLE, width, height) +def test_tabulate_tabulated_err(): + with pytest.raises( + SchematisationError, match=r".*of tabulated type must have increasing heights.*" + ): + tabulate_tabulated(CrossSectionShape.TABULATED_RECTANGLE, [1, 1], [2, 1]) def test_tabulate_inverted_egg(): shape, width_1d, height_1d, table, yz = tabulate_inverted_egg( - "my-shape", "1.52", "ignored" + "my-shape", 1.52, "ignored" ) assert shape == CrossSectionShape.TABULATED_TRAPEZIUM @@ -200,11 +196,29 @@ def test_tabulate_inverted_egg(): @pytest.mark.parametrize( - "width,height,friction_values,vegetation_stem_densities,vegetation_stem_diameters,vegetation_heights,vegetation_drag_coefficients,exp_width,exp_height,exp_table,exp_yz", + "shape, expected_width, expected_height", + [ + (CrossSectionShape.TABULATED_YZ, [1, 3, 5], [2, 4, 6]), + (CrossSectionShape.TABULATED_RECTANGLE, [2, 4, 6], [1, 3, 5]), + ], +) +def test_parse_tabulated(shape, expected_width, expected_height): + width, height = _parse_tabulated("1,2\n3,4\n5,6", shape) + np.testing.assert_array_equal(width, expected_width) + np.testing.assert_array_equal(height, expected_height) + + +@pytest.mark.parametrize("cross_section_table", ["", "1,2\n3", "1,2\n3,"]) +def test_parse_tabulated_invalid(cross_section_table): + with pytest.raises(SchematisationError): + _parse_tabulated(cross_section_table, CrossSectionShape.TABULATED_RECTANGLE) + + +@pytest.mark.parametrize( + "cross_section_table,friction_values,vegetation_stem_densities,vegetation_stem_diameters,vegetation_heights,vegetation_drag_coefficients,exp_width,exp_height,exp_table,exp_yz", [ ( - "0 0.5 1 1.5", - "0.5 0 0 0.5", + "0,0.5\n0.5,0\n1,0\n1.5,0.5", None, None, None, @@ -221,8 +235,7 @@ def test_tabulate_inverted_egg(): ], ), ( - "0 0.5 1 1.5", - "0.5 0 0 0.25", + "0,0.5\n0.5,0\n1,0\n1.5,0.25", "1 1 1", "0.5 1 1", "1 1 1", @@ -239,8 +252,7 @@ def test_tabulate_inverted_egg(): ], ), ( - "0 1 2 3 4 5", - "1 0 0.5 0.5 0 1", + "0,1\n1,0\n2,0.5\n3,0.5\n4,0\n5,1", "1 1 1 1 1", "1 1 1 1 1", "1 1 1 1 1", @@ -259,8 +271,7 @@ def test_tabulate_inverted_egg(): ], ), ( - "0 1 2 2 0 0", - "0.5 0 0.5 1.5 1.5 0.5", + "0,0.5\n1,0\n2,0.5\n2,1.5\n0,1.5\n0,0.5", "1 1 1 1", "1 1 1 1", "1 1 1 1", @@ -272,8 +283,7 @@ def test_tabulate_inverted_egg(): None, ), ( - "0 0.5 0.75 1.0 1.5", - "0.5 0 0 0 0.5", + "0,0.5\n0.5,0\n0.75,0\n1.0,0\n1.5,0.5", "1 1 1 1", "1 0.1 1 1", "2 1 1 1", @@ -291,8 +301,7 @@ def test_tabulate_inverted_egg(): ], ), ( - "0 1 0 1 0", - "0 1 1 0 0", + "0,0\n1,1\n0,1\n1,0\n0,0", "1 1 1 1", "1 1 1 1", "1 1 1 1", @@ -304,8 +313,7 @@ def test_tabulate_inverted_egg(): None, ), ( # Open profile left side higher than right side - "0 1 2 3 4", - "1 0 0 0 2", + "0,1\n1,0\n2,0\n3,0\n4,2", None, None, None, @@ -330,8 +338,7 @@ def test_tabulate_inverted_egg(): ), ), ( # Open profile right side higher than left side - "1 2 3 4 5", - "3 1 0 1 2", + "1,3\n2,1\n3,0\n4,1\n5,2", None, None, None, @@ -356,8 +363,7 @@ def test_tabulate_inverted_egg(): ), ), ( # Open profile same height left and right - "1 2 3 4 5", - "3 1 0 1 3", + "1,3\n2,1\n3,0\n4,1\n5,3", None, None, None, @@ -382,8 +388,7 @@ def test_tabulate_inverted_egg(): ), ), ( # Open profile, left side rises then falls - "0 5 10 15 20 25 30 35 40", - "5.2 5.26 5.25 5.2 5.1 5.1 0 5 5", + "0,5.2\n5,5.26\n10,5.25\n15,5.2\n20,5.1\n25,5.1\n30,0\n35,5\n40,5", None, None, None, @@ -410,8 +415,7 @@ def test_tabulate_inverted_egg(): ], ) def test_tabulate_yz( - width, - height, + cross_section_table, friction_values, vegetation_stem_densities, vegetation_stem_diameters, @@ -422,6 +426,9 @@ def test_tabulate_yz( exp_table, exp_yz, ): + width, height = _parse_tabulated( + cross_section_table, CrossSectionShape.TABULATED_YZ + ) shape, width_1d, height_1d, table, yz = tabulate_yz( "my-shape", width, @@ -432,39 +439,167 @@ def test_tabulate_yz( assert height_1d == exp_height assert_almost_equal(table, np.array(exp_table, dtype=float)) if yz is not None: - 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 +class TestSetFrictionVegetationValues: + def test_set_friction_values(self): + friction_values = "0.5, 0.6" + yz = np.zeros((3, 5), dtype=float) + result = set_friction_vegetation_values( + yz, friction_values, None, None, None, None, None + ) + expected = np.zeros_like(yz) + expected[:, 2] = [0.5, 0.6, 0] + np.testing.assert_array_equal(result, expected) + + def test_set_with_cross_section_vegetation_table(self): + vegetation_table = "1,2,3,4\n10,20,30,40" + yz = np.zeros((3, 5), dtype=float) + result = set_friction_vegetation_values( + yz, None, None, None, None, None, vegetation_table + ) + expected = np.zeros_like(yz) + expected[:, 3] = [8, 8000, 0] + expected[:, 4] = [3, 30, 0] + np.testing.assert_array_equal(result, expected) + + def test_with_cross_section_vegetation_values(self): + vegetation_stem_density = 2.0 + vegetation_stem_diameter = 3.0 + vegetation_drag_coefficient = 0.5 + vegetation_height = 1.0 + yz = np.zeros((2, 5), dtype=float) + result = set_friction_vegetation_values( + yz, + None, + vegetation_stem_density, + vegetation_stem_diameter, + vegetation_height, + vegetation_drag_coefficient, + None, + ) + expected = [[0, 0, 0, 3, 1], [0, 0, 0, 0, 0]] + np.testing.assert_array_equal(result, expected) + + @pytest.mark.parametrize( "width,height,match", [ - ("", "1", r"Unable to parse cross section definition.*"), - ("1", "", r"Unable to parse cross section definition.*"), - ("", "", r"Unable to parse cross section definition.*"), - ("1", "1 2", r".*of tabulated or profile type must have equal number.*"), - ("1 2", "1", r".*of tabulated or profile type must have equal number.*"), ( - "0 0.5 0", - "0 1 0", + [0, 0.5, 0], + [0, 1, 0], r".*of closed profiles must have at least 4 coordinates.*", ), - ("0 0.5", "0 1", r".*of open profiles must have at least 3 coordinates.*"), - ("0 0.5 1.0", "0 1 -1", r".*cannot have negative height coordinate.*"), - ("0 0.5 1.0", "1 2 1", r".*must have at least one height coordinate at 0.*"), - ("0 0.5 1.0 0.5", "0 1 2 3", r".*should be closed or have increasing widths.*"), + ([0, 0.5], [0, 1], r".*of open profiles must have at least 3 coordinates.*"), + ([0, 0.5, 1.0], [0, 1, -1], r".*cannot have negative height coordinate.*"), + ( + [0, 0.5, 1.0], + [1, 2, 1], + r".*must have at least one height coordinate at 0.*", + ), + ( + [0, 0.5, 1.0, 0.5], + [0, 1, 2, 3], + r".*should be closed or have increasing widths.*", + ), ], ) def test_tabulate_yz_err(width, height, match): with pytest.raises(SchematisationError, match=match): tabulate_yz("my-shape", width, height) + + +class TestCrossSectionDefinitionGetUnique: + def test_for_closed_rectangle(self): + csd_in = CrossSectionDefinitions( + id=[100, 200], + shape=[ + constants.CrossSectionShape.CLOSED_RECTANGLE.value, + constants.CrossSectionShape.CLOSED_RECTANGLE.value, + ], + width=[1, 1], + height=[1, 1], + cross_section_table=["foo", "bar"], + origin_table=["pipe", "weir"], + origin_id=[10, 20], + ) + unique_definition, _ = csd_in.get_unique() + assert unique_definition.id == [ + 0, + ] + assert unique_definition.width == [ + 1, + ] + assert unique_definition.height == [ + 1, + ] + + @pytest.mark.parametrize( + "shape", + [ + constants.CrossSectionShape.RECTANGLE.value, + constants.CrossSectionShape.CIRCLE.value, + constants.CrossSectionShape.EGG.value, + constants.CrossSectionShape.INVERTED_EGG.value, + ], + ) + def test_for_tabulated_shape(self, shape): + csd_in = CrossSectionDefinitions( + id=[100, 200], + shape=[shape, shape], + width=[1, 1], + height=[1, 2], + cross_section_table=["foo", "bar"], + origin_table=["pipe", "weir"], + origin_id=[10, 20], + ) + unique_definition, _ = csd_in.get_unique() + assert unique_definition.id == [ + 0, + ] + assert unique_definition.width == [ + 1, + ] + + @pytest.mark.parametrize( + "shape", + [ + constants.CrossSectionShape.TABULATED_YZ.value, + constants.CrossSectionShape.TABULATED_RECTANGLE.value, + constants.CrossSectionShape.TABULATED_TRAPEZIUM.value, + ], + ) + def test_for_other_shapes(self, shape): + csd_in = CrossSectionDefinitions( + id=[100, 200], + shape=[shape, shape], + width=[1, 1], + height=[10, 21], + cross_section_table=["foo", "foo"], + origin_table=["pipe", "weir"], + origin_id=[10, 20], + ) + unique_definition, _ = csd_in.get_unique() + assert unique_definition.id == [ + 0, + ] + assert unique_definition.cross_section_table == ["foo"] + + def test_mapping(self): + csd_in = CrossSectionDefinitions( + id=[100, 200, 300], + shape=[ + constants.CrossSectionShape.CLOSED_RECTANGLE.value, + constants.CrossSectionShape.CLOSED_RECTANGLE.value, + constants.CrossSectionShape.RECTANGLE.value, + ], + width=[1, 1, 3], + height=[1, 1, 10], + origin_table=["pipe", "weir", "pipe"], + origin_id=[10, 20, 20], + ) + _, mapping = csd_in.get_unique() + assert mapping == {"pipe": {10: 0, 20: 1}, "weir": {20: 0}} diff --git a/threedigrid_builder/tests/test_db.py b/threedigrid_builder/tests/test_db.py index 3ee81830..7b03980d 100644 --- a/threedigrid_builder/tests/test_db.py +++ b/threedigrid_builder/tests/test_db.py @@ -4,6 +4,7 @@ import pytest import shapely from shapely.testing import assert_geometries_equal +from threedi_schema import constants from threedigrid_builder.base import GridSettings, Pumps, TablesSettings from threedigrid_builder.constants import ( @@ -28,6 +29,7 @@ Weirs, ) from threedigrid_builder.interface import SQLite +from threedigrid_builder.interface.db import map_cross_section_definition def test_init(tmp_path): @@ -37,7 +39,7 @@ def test_init(tmp_path): with mock.patch( "threedigrid_builder.interface.db.ThreediDatabase" ) as db, mock.patch.object(SQLite, "get_version") as get_version: - get_version.return_value = 227 + get_version.return_value = 228 sqlite = SQLite(path) db.assert_called_with(path) @@ -65,7 +67,7 @@ def test_init_bad_version(tmp_path): def test_get_version(db): - assert db.get_version() == 227 + assert db.get_version() == 228 def test_get_boundary_conditions_1d(db): @@ -86,7 +88,6 @@ def test_get_boundary_conditions_2d(db): assert len(boundary_conditions_2d) == 0 -# TODO: fix this, no clue why this fails def test_get_channels(db): channels = db.get_channels() assert isinstance(channels, Channels) @@ -127,15 +128,12 @@ def test_get_connection_nodes(db): assert connection_nodes.code[494] == "" assert connection_nodes.initial_waterlevel[169] == -0.4 # manhole fields - assert connection_nodes.manhole_id[10] == 11 - assert connection_nodes.manhole_id[100] == -9999 assert connection_nodes.calculation_type[1] == CalculationType.CONNECTED assert connection_nodes.manhole_indicator[6] == 1 assert connection_nodes.bottom_level[9] == -3.51 + assert np.isnan(connection_nodes.bottom_level[100]) assert connection_nodes.drain_level[1] == -0.82 assert connection_nodes.surface_level[35] == -0.54 - assert connection_nodes.shape[40] == "00" - assert connection_nodes.width[32] == 0.8 assert connection_nodes.display_name[33] == "71512" @@ -144,11 +142,23 @@ def test_get_cross_section_definitions(db): assert isinstance(definitions, CrossSectionDefinitions) # some test samples - assert len(definitions.id) == 13 - assert definitions.id[8] == 97 + assert len(definitions.id) == 1365 + assert definitions.id[8] == 8 assert definitions.shape[7] == CrossSectionShape.RECTANGLE - assert definitions.height[10] == "0.4" - assert definitions.width[2] == "0.315" + assert definitions.height[10] == 0.4 + assert definitions.width[2] == 3 + mask_tabulated_yz = ( + definitions.shape == constants.CrossSectionShape.TABULATED_YZ.value + ) + idx_tabulated_yz = definitions.id_to_index(definitions.id[mask_tabulated_yz])[0] + assert ( + definitions.cross_section_table[idx_tabulated_yz] + == "0,2\n1,1\n2,0.5\n3,0\n4,1\n5,1.5\n6,2" + ) + assert ( + definitions.cross_section_vegetation_table[idx_tabulated_yz] + == "10,0.2,1,7\n5,0.1,0.5,7\n5,0.1,0.5,7\n5,0.3,0.5,7\n15,0.4,1,7\n20,0.1,1.5,7" + ) def test_get_cross_section_locations(db): @@ -161,7 +171,6 @@ def test_get_cross_section_locations(db): locations.the_geom[96], shapely.from_wkt("POINT (111104 521655)"), tolerance=1 ) assert locations.id[11] == 12 - assert locations.definition_id[365] == 98 assert locations.channel_id[448] == 452 assert locations.reference_level[691] == -3.0 assert locations.bank_level[995] == -1.7 @@ -219,14 +228,12 @@ def test_get_pipes(db): assert pipes.calculation_type[3] == CalculationType.ISOLATED assert pipes.connection_node_start_id[4] == 37 assert pipes.connection_node_end_id[5] == 35 - assert pipes.cross_section_definition_id[9] == 7 assert pipes.invert_level_start_point[16] == -3.91 assert pipes.invert_level_end_point[19] == -3.62 assert pipes.sewerage_type[24] == 2 assert pipes.friction_type[28] == FrictionType.MANNING assert pipes.friction_value[36] == 0.0145 assert pipes.display_name[33] == "71518_71517" - assert pipes.zoom_category[15] == 3 def test_get_settings(db): @@ -303,7 +310,6 @@ def test_get_settings(db): def test_get_pumps(db): pumps = db.get_pumps() assert isinstance(pumps, Pumps) - # some test samples assert len(pumps) == 19 assert pumps.id[11] == 13 @@ -316,7 +322,6 @@ def test_get_pumps(db): assert pumps.start_level[0] == -4.0 assert pumps.lower_stop_level[18] == -1.9 assert np.isnan(pumps.upper_stop_level[15]) - assert pumps.zoom_category[15] == 5 assert pumps.display_name[10] == "KGM-JL-18" @@ -338,7 +343,6 @@ def test_get_culverts(db): assert culverts.connection_node_start_id[45] == 1220 assert culverts.connection_node_end_id[61] == 1620 assert culverts.calculation_type[4] == CalculationType.ISOLATED - assert culverts.cross_section_definition_id[0] == 97 assert culverts.invert_level_start_point[21] == -2.39 assert culverts.invert_level_end_point[83] == -3.28 assert culverts.discharge_coefficient_negative[0] == 0.8 @@ -346,7 +350,6 @@ def test_get_culverts(db): assert culverts.friction_type[28] == FrictionType.MANNING assert culverts.friction_value[36] == 0.03 assert culverts.display_name[32] == "KDU-Q-18482" - assert culverts.zoom_category[15] == 2 def test_get_orifices(db): @@ -370,13 +373,11 @@ def test_get_weirs(db): assert weirs.connection_node_end_id[7] == 394 assert weirs.crest_level[26] == -2.508 assert weirs.crest_type[0] == CalculationType.SHORT_CRESTED - assert weirs.cross_section_definition_id[0] == 8 assert weirs.discharge_coefficient_negative[0] == 0.0 assert weirs.discharge_coefficient_positive[0] == 0.8 assert weirs.friction_type[28] == FrictionType.MANNING assert weirs.friction_value[36] == 0.03 assert weirs.display_name[33] == "KST-JL-76" - assert weirs.zoom_category[15] == 2 assert weirs.sewerage[0] == 1 @@ -407,3 +408,15 @@ def test_get_exchange_lines(db): exchange_lines = db.get_exchange_lines() # No exchange lines in test dataset assert len(exchange_lines) == 0 + + +def test_map_cross_section_definition(): + mapping = {"pipe": {0: 0, 1: 1}, "weir": {1: 0}, "cross_section_location": {0: 10}} + pipes = Pipes(id=[0, 1], cross_section_definition_id=[10, 20]) + weirs = Weirs(id=[0, 1], cross_section_definition_id=[10, 20]) + cross_section_locations = CrossSectionLocations(id=[0], definition_id=[100]) + objects = [pipes, weirs, cross_section_locations] + map_cross_section_definition(objects, mapping) + np.testing.assert_array_equal(pipes.cross_section_definition_id, [0, 1]) + np.testing.assert_array_equal(weirs.cross_section_definition_id, [10, 0]) + np.testing.assert_array_equal(cross_section_locations.definition_id, [10]) diff --git a/threedigrid_builder/tests/test_grid.py b/threedigrid_builder/tests/test_grid.py index 26fcc89a..78301941 100644 --- a/threedigrid_builder/tests/test_grid.py +++ b/threedigrid_builder/tests/test_grid.py @@ -48,6 +48,7 @@ def meta(): grid_space=20.0, dist_calc_points=25.0, kmax=4, + node_open_water_detection=1, ), tables_settings=TablesSettings( table_step_size=0.05, diff --git a/threedigrid_builder/tests/test_gridadmin.py b/threedigrid_builder/tests/test_gridadmin.py index ff8bbc45..bb1cf1cc 100644 --- a/threedigrid_builder/tests/test_gridadmin.py +++ b/threedigrid_builder/tests/test_gridadmin.py @@ -50,6 +50,7 @@ def h5_out_1d(tmpdir_factory, grid_all): ("id", (4,), "int32"), ("initial_waterlevel", (4,), "float64"), ("is_manhole", (4,), "int32"), + ("manhole_id", (4,), "int32"), ("manhole_indicator", (4,), "int32"), ("node_type", (4,), "int32"), ("pixel_coords", (4, 4), "int32"), @@ -89,6 +90,7 @@ def test_write_nodes(h5_out, dataset, shape, dtype): ("id", (3,), "int32"), ("initial_waterlevel", (3,), "float64"), ("is_manhole", (3,), "int32"), + ("manhole_id", (3,), "int32"), ("manhole_indicator", (3,), "int32"), ("node_type", (3,), "int32"), ("pixel_coords", (4, 3), "int32"), diff --git a/threedigrid_builder/tests/test_lines_1d2d.py b/threedigrid_builder/tests/test_lines_1d2d.py index 52b9137b..a57aaf54 100644 --- a/threedigrid_builder/tests/test_lines_1d2d.py +++ b/threedigrid_builder/tests/test_lines_1d2d.py @@ -303,31 +303,47 @@ def test_assign_dpumax_from_breaches(assign_dpumax): @mock.patch.object(Lines1D2D, "assign_dpumax") @mock.patch.object(Lines1D2D, "assign_ds1d_half_from_obstacles") +@pytest.mark.parametrize("test_case", ["open", "closed"]) def test_assign_dpumax_from_obstacles( assign_ds1d_half_from_obstacles, assign_dpumax, + test_case, ): obstacles = mock.Mock() + obstacles.compute_dpumax.return_value = ( np.array([1.2, np.nan]), np.array([1, -9999]), ) - - lines = Lines1D2D( - id=range(3), - kcu=[ - LineType.LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER, - LineType.LINE_1D2D_SINGLE_CONNECTED_CLOSED, - LineType.LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER, - ], - ) - - lines.assign_dpumax_from_obstacles(obstacles) + if test_case == "open": + obstacles.affects_1d2d_open_water = np.array([True, True, True]) + obstacles.affects_1d2d_closed = np.array([False, False, False]) + lines = Lines1D2D( + id=range(3), + kcu=[ + LineType.LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER, + LineType.LINE_1D2D_SINGLE_CONNECTED_CLOSED, + LineType.LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER, + ], + ) + lines.assign_dpumax_from_obstacles_open(obstacles) + elif test_case == "closed": + obstacles.affects_1d2d_closed = np.array([True, True, True]) + obstacles.affects_1d2d_open_water = np.array([False, False, False]) + lines = Lines1D2D( + id=range(3), + kcu=[ + LineType.LINE_1D2D_SINGLE_CONNECTED_CLOSED, + LineType.LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER, + LineType.LINE_1D2D_SINGLE_CONNECTED_CLOSED, + ], + ) + lines.assign_dpumax_from_obstacles_closed(obstacles) args, kwargs = obstacles.compute_dpumax.call_args assert len(args) == 1 assert args[0] is lines - assert len(kwargs) == 1 + assert len(kwargs) == 2 assert_array_equal(kwargs["where"], [0, 2]) (actual_mask, actual_dpumax), _ = assign_dpumax.call_args diff --git a/threedigrid_builder/tests/test_obstacles.py b/threedigrid_builder/tests/test_obstacles.py index 067312bc..c47fb7a1 100644 --- a/threedigrid_builder/tests/test_obstacles.py +++ b/threedigrid_builder/tests/test_obstacles.py @@ -5,7 +5,7 @@ from threedigrid_builder.base import Lines, Nodes from threedigrid_builder.constants import LineType -from threedigrid_builder.grid import Obstacles +from threedigrid_builder.grid import ObstacleAffectsType, Obstacles from threedigrid_builder.grid.grid import Grid @@ -44,6 +44,9 @@ def obstacles(): the_geom=shapely.linestrings( [[[12.0, 11.0], [17.5, 17.0]], [[13.0, 11.0], [13.0, 17.0]]] ), + affects_2d=[True, True], + affects_1d2d_open_water=[True, True], + affects_1d2d_closed=[True, True], ) @@ -55,11 +58,14 @@ def obstacles_no_intersect(): the_geom=shapely.linestrings( [[[11, 11], [11, 18], [18, 18], [18, 11], [11, 11]]] ), + affects_2d=[True], + affects_1d2d_open_water=[True], + affects_1d2d_closed=[True], ) def test_set_obstacles(grid, obstacles): - grid.set_obstacles(obstacles) + grid.set_obstacles_2d(obstacles) expected_kcu = np.array( [102, 102, 102, 102, 98, 99, 103, 99, 99, 103], dtype=np.int32 @@ -74,7 +80,7 @@ def test_set_obstacles(grid, obstacles): def test_set_obstacles_no_intersect(grid, obstacles_no_intersect): - grid.set_obstacles(obstacles_no_intersect) + grid.set_obstacles_2d(obstacles_no_intersect) assert np.all(grid.lines.kcu != LineType.LINE_2D_OBSTACLE) assert_equal(grid.lines.flod, np.nan) @@ -83,8 +89,9 @@ def test_set_obstacles_no_intersect(grid, obstacles_no_intersect): def test_compute_dpumax(lines: Lines, obstacles: Obstacles): expected = [-0.1, -0.3, -0.3, -0.1, np.nan, np.nan, -0.1, np.nan, np.nan, -0.1] - actual, idx = obstacles.compute_dpumax(lines, where=np.arange(len(lines))) - + actual, idx = obstacles.compute_dpumax( + lines, where=np.arange(len(lines)), affects_type=ObstacleAffectsType.AFFECTS_2D + ) assert_equal(idx, [0, 1, 1, 0, -9999, -9999, 0, -9999, -9999, 0]) assert_almost_equal(actual, expected) @@ -93,7 +100,7 @@ def test_compute_dpumax_no_intersections( lines: Lines, obstacles_no_intersect: Obstacles ): actual, idx = obstacles_no_intersect.compute_dpumax( - lines, where=np.arange(len(lines)) + lines, where=np.arange(len(lines)), affects_type=ObstacleAffectsType.AFFECTS_2D ) assert len(actual) == len(lines) @@ -102,16 +109,43 @@ def test_compute_dpumax_no_intersections( def test_compute_dpumax_no_obstacles(lines: Lines): - actual, idx = Obstacles(id=[]).compute_dpumax(lines, where=np.arange(len(lines))) + actual, idx = Obstacles(id=[]).compute_dpumax( + lines, where=np.arange(len(lines)), affects_type=ObstacleAffectsType.AFFECTS_2D + ) assert len(actual) == len(lines) assert_equal(idx, -9999) assert_almost_equal(actual, np.nan) +@pytest.mark.parametrize("affects_type", list(ObstacleAffectsType)) +def test_compute_dpumax_affects_some(obstacles: Obstacles, lines: Lines, affects_type): + affects_args = { + obstacle_type.value: np.full(2, True, dtype=bool) + for obstacle_type in ObstacleAffectsType + } + affects_args[affects_type.value] = [True, False] + obstacles = Obstacles( + id=[1, 2], + crest_level=[-0.1, -0.3], + the_geom=shapely.linestrings( + [[[12.0, 11.0], [17.5, 17.0]], [[13.0, 11.0], [13.0, 17.0]]] + ), + **affects_args + ) + actual, idx = obstacles.compute_dpumax( + lines, where=[0, 1], affects_type=affects_type + ) + assert len(actual) == 2 + assert_equal(idx, [0, -9999]) + assert_almost_equal(actual, [-0.1, np.nan]) + + def test_compute_dpumax_where(lines: Lines, obstacles: Obstacles): expected = [-0.1, -0.1, np.nan, -0.1] - actual, idx = obstacles.compute_dpumax(lines, where=[0, 3, 4, 6]) + actual, idx = obstacles.compute_dpumax( + lines, where=[0, 3, 4, 6], affects_type=ObstacleAffectsType.AFFECTS_2D + ) assert_equal(idx, [0, 0, -9999, 0]) assert_almost_equal(actual, expected)