Fix dommel network (#129)

Contains all work for #102, including:
- NetworkValidator to automated fixing of topological issues
- Misc new methods for the Ribasim-NL Model class:
  - `reverse_edge`: reverse an edge and it's `from` and `to`properties
- `remove_node` remove a node, possibility to remove connected_edges as
well
- `remove_edge`: remove an edge, possibility to remove
disconnected_nodes as well
- the modify_model.py script to work from Sweco 2024.6.3 model to
2024.8.0 available at https://deltares.thegood.cloud/f/117827

All results in a topologically correct network that passes the RIBASIM
network validator and simulation. Note, all parameters are implausible
still:

![image](https://github.com/user-attachments/assets/5314e68f-ce9b-45aa-bb12-eb7e7243aa38)

---------

Co-authored-by: Martijn Visser <[email protected]>

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)

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)

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)

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")
+# %%

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"]

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