diff --git a/notebooks/aaenmaas/get_model.py b/notebooks/aaenmaas/get_model.py new file mode 100644 index 0000000..29dd5a3 --- /dev/null +++ b/notebooks/aaenmaas/get_model.py @@ -0,0 +1,9 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +aaenmaas_url = cloud.joinurl("AaenMaas", "modellen", "AaenMaas_2024_6_3") + +# %% +cloud.download_content(aaenmaas_url) diff --git a/notebooks/de_dommel/get_model.py b/notebooks/de_dommel/get_model.py new file mode 100644 index 0000000..56538ad --- /dev/null +++ b/notebooks/de_dommel/get_model.py @@ -0,0 +1,9 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +dommel_url = cloud.joinurl("DeDommel", "modellen", "DeDommel_2024_6_3") + +# %% +cloud.download_content(dommel_url) diff --git a/notebooks/de_dommel/get_verwerkt.py b/notebooks/de_dommel/get_verwerkt.py new file mode 100644 index 0000000..f6d7d10 --- /dev/null +++ b/notebooks/de_dommel/get_verwerkt.py @@ -0,0 +1,9 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +dommel_url = cloud.joinurl("DeDommel", "verwerkt") + +# %% +cloud.download_content(dommel_url) diff --git a/notebooks/de_dommel/modify_model.py b/notebooks/de_dommel/modify_model.py new file mode 100644 index 0000000..8e636d5 --- /dev/null +++ b/notebooks/de_dommel/modify_model.py @@ -0,0 +1,210 @@ +# %% +import sqlite3 + +import geopandas as gpd +import pandas as pd +from ribasim import Node +from ribasim.nodes import basin, level_boundary, manning_resistance, outlet +from ribasim_nl import CloudStorage, Model, NetworkValidator + +cloud = CloudStorage() + +ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_2024_6_3", "model.toml") +database_gpkg = ribasim_toml.with_name("database.gpkg") + + +basin_data = [ + basin.Profile(level=[0.0, 1.0], area=[0.01, 1000.0]), + basin.Static( + drainage=[0.0], + potential_evaporation=[0.001 / 86400], + infiltration=[0.0], + precipitation=[0.005 / 86400], + ), + basin.State(level=[0]), +] + +# %% remove urban_runoff +# Connect to the SQLite database + +conn = sqlite3.connect(database_gpkg) + +# get table into DataFrame +table = "Basin / static" +df = pd.read_sql_query(f"SELECT * FROM '{table}'", conn) + +# drop urban runoff column if exists +if "urban_runoff" in df.columns: + df.drop(columns="urban_runoff", inplace=True) + + # Write the DataFrame back to the SQLite table + df.to_sql(table, conn, if_exists="replace", index=False) + +# # Close the connection +conn.close() + + +# %% read model +model = Model.read(ribasim_toml) + +network_validator = NetworkValidator(model) + +# %% verwijder duplicated edges + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2288780504 +model.edge.df = model.edge.df[~((model.edge.df.from_node_type == "Basin") & (model.edge.df.to_node_type == "Basin"))] + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291081244 +duplicated_fids = network_validator.edge_duplicated().index.to_list() +model.edge.df = model.edge.df.drop_duplicates(subset=["from_node_id", "to_node_id"]) + +if not network_validator.edge_duplicated().empty: + raise Exception("nog steeds duplicated edges") + +# %% toevoegen bovenstroomse knopen + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291091067 +edge_mask = model.edge.df.index.isin(duplicated_fids) +node_id = model.next_node_id +edge_fid = next(i for i in duplicated_fids if i in model.edge.df.index) +model.edge.df.loc[edge_mask, ["from_node_type"]] = "Basin" +model.edge.df.loc[edge_mask, ["from_node_id"]] = node_id + +node = Node(node_id, model.edge.df.at[edge_fid, "geometry"].boundary.geoms[0]) +model.basin.area.df.loc[model.basin.area.df.node_id == 1009, ["node_id"]] = node_id +area = basin.Area(geometry=model.basin.area[node_id].geometry.to_list()) +model.basin.add(node, basin_data + [area]) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291111647 + + +for row in network_validator.edge_incorrect_connectivity().itertuples(): + # drop edge from model + model.remove_edge(row.from_node_id, row.to_node_id, remove_disconnected_nodes=False) + + # add basin_node + area = basin.Area(geometry=model.basin.area[row.from_node_id].geometry.to_list()) + basin_node = Node(row.from_node_id, row.geometry.boundary.geoms[0]) + model.basin.add(basin_node, basin_data + [area]) + + # eindhovensch kanaal we need to add manning a 99% of the length + if row.to_node_id == 2: + geometry = row.geometry.interpolate(0.99, normalized=True) + name = "" + if row.to_node_id == 14: + gdf = gpd.read_file( + cloud.joinpath("DeDommel", "verwerkt", "1_ontvangen_data", "Geodata", "data_Q42018.gpkg"), + layer="DuikerSifonHevel", + engine="pyogrio", + fid_as_index=True, + ) + kdu = gdf.loc[5818] + geometry = kdu.geometry.interpolate(0.5, normalized=True) + name = kdu.objectId + + # add manning-node + manning_node_id = model.next_node_id + manning_data = manning_resistance.Static(length=[100], manning_n=[0.04], profile_width=[10], profile_slope=[1]) + model.manning_resistance.add( + Node(node_id=manning_node_id, geometry=geometry, name=name), + [manning_data], + ) + + # add edges + model.edge.add(model.basin[row.from_node_id], model.manning_resistance[manning_node_id]) + model.edge.add(model.manning_resistance[manning_node_id], model.level_boundary[row.to_node_id]) + + +if not network_validator.edge_incorrect_connectivity().empty: + raise Exception("nog steeds edges zonder knopen") + +# %% verwijderen internal basins + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291271525 +for row in network_validator.node_internal_basin().itertuples(): + if row.node_id not in model.basin.area.df.node_id.to_numpy(): # remove or change to level-boundary + edge_select_df = model.edge.df[model.edge.df.to_node_id == row.node_id] + if len(edge_select_df) == 1: + if edge_select_df.iloc[0]["from_node_type"] == "FlowBoundary": + model.remove_node(row.node_id) + model.remove_node(edge_select_df.iloc[0]["from_node_id"]) + model.edge.df.drop(index=edge_select_df.index[0], inplace=True) + +df = model.node_table().df[model.node_table().df.node_type == "Basin"] + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291876800 +boundary_node = Node(node_id=28, geometry=model.flow_boundary[28].geometry) + +for node_id in [29, 1828, 615, 28, 1329]: + model.remove_node(node_id, remove_edges=True) + +level_data = level_boundary.Static(level=[0]) +model.level_boundary.add(boundary_node, [level_data]) + +model.edge.add(model.tabulated_rating_curve[614], model.level_boundary[28]) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292014475 +model.remove_node(node_id=1898, remove_edges=True) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292017813 +for node_id in [1891, 989, 1058]: + model.remove_node(node_id, remove_edges=True) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291988317 +# for from_node_id, to_node_id in [799, 1580, 625, 1123, 597, 978]: +for from_node_id, to_node_id in [[616, 1032], [1030, 616], [393, 1242], [1852, 393], [353, 1700], [1253, 353]]: + model.reverse_edge(from_node_id, to_node_id) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292050862 + +gdf = gpd.read_file( + cloud.joinpath("DeDommel", "verwerkt", "1_ontvangen_data", "Geodata", "data_Q42018.gpkg"), + layer="HydroObject", + engine="pyogrio", + fid_as_index=True, +) + +geometry = gdf.loc[2751].geometry.interpolate(0.5, normalized=True) +node_id = model.next_node_id +data = manning_resistance.Static(length=[100], manning_n=[0.04], profile_width=[10], profile_slope=[1]) + +model.manning_resistance.add(Node(node_id=node_id, geometry=geometry), [manning_data]) + +model.edge.df = model.edge.df[~((model.edge.df.from_node_id == 611) & (model.edge.df.to_node_id == 1643))] +model.edge.add(model.basin[1643], model.manning_resistance[node_id]) +model.edge.add(model.manning_resistance[node_id], model.basin[1182]) +model.edge.add(model.tabulated_rating_curve[611], model.basin[1182]) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2293457160 +model.update_node(417, "Outlet", [outlet.Static(flow_rate=[0])]) + +if not network_validator.node_internal_basin().empty: + raise Exception("nog steeds interne basins") + +df = network_validator.edge_incorrect_type_connectivity( + from_node_type="ManningResistance", to_node_type="LevelBoundary" +) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2293486609 +for node_id in df.from_node_id: + model.update_node(node_id, "Outlet", [outlet.Static(flow_rate=[100])]) + +# for row in df.itertuples(): +# area_select_df = model.basin.area.df[model.basin.area.df.geometry.contains(row.geometry)] +# if len(area_select_df) == 1: +# area_id = area_select_df.iloc[0]["node_id"] +# if row.node_id != area_id: +# print(f"{row.node_id} == {area_id}") +# elif area_select_df.empty: +# raise Exception(f"Basin {row.node_id} not within Area") +# else: +# raise Exception(f"Basin {row.node_id} contained by multiple areas") + +# %% write model +ribasim_toml = ribasim_toml.parents[1].joinpath("DeDommel", ribasim_toml.name) +model.write(ribasim_toml) + +# %% upload model + +# cloud.upload_model("DeDommel", "DeDommel") +# %% diff --git a/src/ribasim_nl/ribasim_nl/__init__.py b/src/ribasim_nl/ribasim_nl/__init__.py index 593c2b2..ea8fb34 100644 --- a/src/ribasim_nl/ribasim_nl/__init__.py +++ b/src/ribasim_nl/ribasim_nl/__init__.py @@ -3,6 +3,7 @@ from ribasim_nl.cloud import CloudStorage from ribasim_nl.model import Model from ribasim_nl.network import Network +from ribasim_nl.network_validator import NetworkValidator from ribasim_nl.reset_index import reset_index -__all__ = ["CloudStorage", "Network", "reset_index", "Model"] +__all__ = ["CloudStorage", "Network", "reset_index", "Model", "NetworkValidator"] diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index 0bea2f8..3141aa2 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -92,6 +92,29 @@ def get_node(self, node_id: int): node_type = self.get_node_type(node_id) return getattr(self, pascal_to_snake_case(node_type))[node_id] + def remove_node(self, node_id: int, remove_edges: bool = False): + """Remove node from model""" + node_type = self.get_node_type(node_id) + + # read existing table + table = getattr(self, pascal_to_snake_case(node_type)) + + # remove node from all tables + for attr in table.model_fields.keys(): + df = getattr(table, attr).df + if df is not None: + getattr(table, attr).df = df[df.node_id != node_id] + + if remove_edges and (self.edge.df is not None): + for row in self.edge.df[self.edge.df.from_node_id == node_id].itertuples(): + self.remove_edge( + from_node_id=row.from_node_id, to_node_id=row.to_node_id, remove_disconnected_nodes=False + ) + for row in self.edge.df[self.edge.df.to_node_id == node_id].itertuples(): + self.remove_edge( + from_node_id=row.from_node_id, to_node_id=row.to_node_id, remove_disconnected_nodes=False + ) + def update_node(self, node_id, node_type, data, node_properties: dict = {}): """Update a node type and/or data""" # get existing network node_type @@ -102,6 +125,7 @@ def update_node(self, node_id, node_type, data, node_properties: dict = {}): # save node, so we can add it later node_dict = table.node.df[table.node.df["node_id"] == node_id].iloc[0].to_dict() + node_dict.pop("node_type") # remove node from all tables for attr in table.model_fields.keys(): @@ -178,6 +202,44 @@ def add_control_node( for _to_node_id in to_node_id: self.edge.add(table[node_id], self.get_node(_to_node_id)) + def reverse_edge(self, from_node_id: int, to_node_id: int): + """Reverse an edge""" + if self.edge.df is not None: + # get original edge-data + df = self.edge.df.copy() + df.loc[:, ["edge_id"]] = df.index + df = df.set_index(["from_node_id", "to_node_id"], drop=False) + edge_data = dict(df.loc[from_node_id, to_node_id]) + edge_id = edge_data["edge_id"] + + # revert node ids + self.edge.df.loc[edge_id, ["from_node_id"]] = edge_data["to_node_id"] + self.edge.df.loc[edge_id, ["to_node_id"]] = edge_data["from_node_id"] + + # revert node types + self.edge.df.loc[edge_id, ["from_node_type"]] = edge_data["to_node_type"] + self.edge.df.loc[edge_id, ["to_node_type"]] = edge_data["from_node_type"] + + # revert geometry + self.edge.df.loc[edge_id, ["geometry"]] = edge_data["geometry"].reverse() + + def remove_edge(self, from_node_id: int, to_node_id: int, remove_disconnected_nodes=True): + """Remove an edge and disconnected nodes""" + if self.edge.df is not None: + # get original edge-data + indices = self.edge.df[ + (self.edge.df.from_node_id == from_node_id) & (self.edge.df.to_node_id == to_node_id) + ].index + + # remove edge from edge-table + self.edge.df = self.edge.df[~self.edge.df.index.isin(indices)] + + # remove disconnected nodes + if remove_disconnected_nodes: + for node_id in [from_node_id, to_node_id]: + if node_id not in self.edge.df[["from_node_id", "to_node_id"]].to_numpy().ravel(): + self.remove_node(node_id) + def find_closest_basin(self, geometry: BaseGeometry, max_distance: float | None) -> NodeData: """Find the closest basin_node.""" # only works when basin area are defined diff --git a/src/ribasim_nl/ribasim_nl/network_validator.py b/src/ribasim_nl/ribasim_nl/network_validator.py new file mode 100644 index 0000000..ba95378 --- /dev/null +++ b/src/ribasim_nl/ribasim_nl/network_validator.py @@ -0,0 +1,116 @@ +from dataclasses import dataclass + +from ribasim import Model + + +def within_distance(row, gdf, tolerance=1.0) -> bool: + distance = gdf[gdf.index != row.name].distance(row.geometry) + return (distance < tolerance).any() + + +def check_node_connectivity(row, node_df, tolerance=1.0) -> bool: + invalid = True + + # check if from_node_id is valid + if row.from_node_id in node_df.index: + distance = row.geometry.boundary.geoms[0].distance(node_df.at[row.from_node_id, "geometry"]) + invalid = distance > tolerance + + # if valid, check if to_node_id is valid + if (not invalid) and (row.to_node_id in node_df.index): + distance = row.geometry.boundary.geoms[1].distance(node_df.at[row.to_node_id, "geometry"]) + invalid = distance > tolerance + + return invalid + + +def check_internal_basin(row, edge_df) -> bool: + if row.node_type == "Basin": + return row.node_id not in edge_df.from_node_id.to_numpy() + else: + return False + + +@dataclass +class NetworkValidator: + """Contains some methods for validating a RIBASIM-model + + Parameters + ---------- + model: ribasim.Model + Ribasim model + tolerance: float (default=1) + Tolerance to use for snapping. Should be in the units of the model.crs + + """ + + model: Model + tolerance: float = 1 + + @property + def node_df(self): + return self.model.node_table().df + + @property + def edge_df(self): + return self.model.edge.df + + def node_overlapping(self): + """Check if the node-geometry overlaps another node within tolerance (default=1m)""" + return self.node_df[self.node_df.apply(lambda row: within_distance(row, self.node_df, self.tolerance), axis=1)] + + def node_duplicated(self): + """Check if node_id is duplicated""" + return self.node_df[self.node_df.node_id.duplicated()] + + def node_internal_basin(self): + """Check if a Node with node_type Basin is not connected to another node""" + mask = self.node_df.apply(lambda row: check_internal_basin(row, self.edge_df), axis=1) + return self.node_df[mask] + + def edge_duplicated(self): + """Check if the `from_node_id` and `to_node_id` in the edge-table is duplicated""" + return self.edge_df[self.edge_df.duplicated(subset=["from_node_id", "to_node_id"], keep=False)] + + def edge_missing_nodes(self): + """Check if the `from_node_id` and `to_node_id` in the edge-table are both as node-id in the node-table""" + mask = ~( + self.edge_df.from_node_id.isin(self.node_df.node_id) & self.edge_df.to_node_id.isin(self.node_df.node_id) + ) + return self.edge_df[mask] + + def edge_incorrect_from_node(self): + """Check if the `from_node_type` in edge-table in matches the `node_type` of the corresponding node in the node-table""" + node_df = self.node_df.set_index("node_id") + mask = ~self.edge_df.apply( + lambda row: node_df.at[row["from_node_id"], "node_type"] == row["from_node_type"] + if row["from_node_id"] in node_df.index + else False, + axis=1, + ) + return self.edge_df[mask] + + def edge_incorrect_to_node(self): + """Check if the `to_node_type` in edge-table in matches the `node_type` of the corresponding node in the node-table""" + node_df = self.node_df.set_index("node_id") + mask = ~self.edge_df.apply( + lambda row: node_df.at[row["to_node_id"], "node_type"] == row["to_node_type"] + if row["to_node_id"] in node_df.index + else False, + axis=1, + ) + return self.edge_df[mask] + + def edge_incorrect_connectivity(self): + """Check if the geometries of the `from_node_id` and `to_node_id` are on the start and end vertices of the edge-geometry within tolerance (default=1m)""" + node_df = self.node_df.set_index("node_id") + mask = self.edge_df.apply( + lambda row: check_node_connectivity(row=row, node_df=node_df, tolerance=self.tolerance), axis=1 + ) + + return self.edge_df[mask] + + def edge_incorrect_type_connectivity(self, from_node_type="ManningResistance", to_node_type="LevelBoundary"): + """Check edges that contain wrong connectivity""" + mask = (self.edge_df.from_node_type == from_node_type) & (self.edge_df.to_node_type == to_node_type) + return self.edge_df[mask]