From 71b67e6eecbf783ea25fc91c12c9744db78ad3e0 Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Fri, 11 Oct 2024 12:36:26 +0200 Subject: [PATCH] Afbouwen wdod vechtstromen (#163) Alle verbeteringen aan model-klasse voor het debuggen van Vechtstromen en WDOD --------- Co-authored-by: Martijn Visser --- notebooks/basin_area_nodes.py | 28 + notebooks/de_dommel/00_get_model.py | 4 +- notebooks/de_dommel/00_get_verwerkt.py | 4 +- .../drents_overijsselse_delta/00_get_model.py | 15 + .../00_get_verwerkt.py | 9 + .../01_fix_model_network.py | 403 +++++++ notebooks/vechtstromen/00_get_model.py | 15 + notebooks/vechtstromen/00_get_verwerkt.py | 9 + .../vechtstromen/01_fix_model_network.py | 1028 +++++++++++++++++ notebooks/vechtstromen/99_upload_model.py | 6 + src/ribasim_nl/ribasim_nl/geometry.py | 49 + src/ribasim_nl/ribasim_nl/model.py | 161 ++- .../ribasim_nl/network_validator.py | 26 + 13 files changed, 1744 insertions(+), 13 deletions(-) create mode 100644 notebooks/basin_area_nodes.py create mode 100644 notebooks/drents_overijsselse_delta/00_get_model.py create mode 100644 notebooks/drents_overijsselse_delta/00_get_verwerkt.py create mode 100644 notebooks/drents_overijsselse_delta/01_fix_model_network.py create mode 100644 notebooks/vechtstromen/00_get_model.py create mode 100644 notebooks/vechtstromen/00_get_verwerkt.py create mode 100644 notebooks/vechtstromen/01_fix_model_network.py create mode 100644 notebooks/vechtstromen/99_upload_model.py diff --git a/notebooks/basin_area_nodes.py b/notebooks/basin_area_nodes.py new file mode 100644 index 0000000..9f809b0 --- /dev/null +++ b/notebooks/basin_area_nodes.py @@ -0,0 +1,28 @@ +# %% + +import pandas as pd +from ribasim_nl import CloudStorage, Model + +cloud = CloudStorage() + + +# %% +data = [] +for authority in cloud.water_authorities: + ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_2024_6_3", "model.toml") + if ribasim_toml.exists(): + model = Model.read(ribasim_toml) + data += [ + { + "waterschap": authority, + "basin_nodes": len(model.basin.node.df), + "basin_areas": len(model.basin.area.df), + "basin_verschil": abs(len(model.basin.node.df) - len(model.basin.area.df)), + "basin_area_lt_5000m2": len(model.basin.area.df[model.basin.area.df.area < 5000]), + } + ] + +df = pd.DataFrame(data) + + +df.to_excel(cloud.joinpath("verschil_basins.xlsx"), index=False) diff --git a/notebooks/de_dommel/00_get_model.py b/notebooks/de_dommel/00_get_model.py index 56538ad..12e044a 100644 --- a/notebooks/de_dommel/00_get_model.py +++ b/notebooks/de_dommel/00_get_model.py @@ -3,7 +3,7 @@ cloud = CloudStorage() -dommel_url = cloud.joinurl("DeDommel", "modellen", "DeDommel_2024_6_3") +model_url = cloud.joinurl("DeDommel", "modellen", "DeDommel_2024_6_3") # %% -cloud.download_content(dommel_url) +cloud.download_content(model_url) diff --git a/notebooks/de_dommel/00_get_verwerkt.py b/notebooks/de_dommel/00_get_verwerkt.py index f6d7d10..5cdc8b3 100644 --- a/notebooks/de_dommel/00_get_verwerkt.py +++ b/notebooks/de_dommel/00_get_verwerkt.py @@ -3,7 +3,7 @@ cloud = CloudStorage() -dommel_url = cloud.joinurl("DeDommel", "verwerkt") +data_url = cloud.joinurl("DeDommel", "verwerkt") # %% -cloud.download_content(dommel_url) +cloud.download_content(data_url) diff --git a/notebooks/drents_overijsselse_delta/00_get_model.py b/notebooks/drents_overijsselse_delta/00_get_model.py new file mode 100644 index 0000000..2cf928b --- /dev/null +++ b/notebooks/drents_overijsselse_delta/00_get_model.py @@ -0,0 +1,15 @@ +# %% +from ribasim_nl import CloudStorage + +authority = "DrentsOverijsselseDelta" +short_name = "dod" + +cloud = CloudStorage() + +model_url = cloud.joinurl(authority, "modellen", f"{authority}_2024_6_3") + +cloud.download_content(model_url) + +ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_2024_6_3", "model.toml") +if ribasim_toml.exists(): + ribasim_toml.rename(ribasim_toml.with_name(f"{short_name}.toml")) diff --git a/notebooks/drents_overijsselse_delta/00_get_verwerkt.py b/notebooks/drents_overijsselse_delta/00_get_verwerkt.py new file mode 100644 index 0000000..5cdc8b3 --- /dev/null +++ b/notebooks/drents_overijsselse_delta/00_get_verwerkt.py @@ -0,0 +1,9 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +data_url = cloud.joinurl("DeDommel", "verwerkt") + +# %% +cloud.download_content(data_url) diff --git a/notebooks/drents_overijsselse_delta/01_fix_model_network.py b/notebooks/drents_overijsselse_delta/01_fix_model_network.py new file mode 100644 index 0000000..fd1aa83 --- /dev/null +++ b/notebooks/drents_overijsselse_delta/01_fix_model_network.py @@ -0,0 +1,403 @@ +# %% +import geopandas as gpd +import numpy as np +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 +from ribasim_nl.geometry import split_basin_multi_polygon + +cloud = CloudStorage() + +authority = "DrentsOverijsselseDelta" +short_name = "dod" + +ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_2024_6_3", f"{short_name}.toml") +database_gpkg = ribasim_toml.with_name("database.gpkg") +hydroobject_gdf = gpd.read_file( + cloud.joinpath(authority, "verwerkt", "4_ribasim", "hydamo.gpkg"), layer="hydroobject", fid_as_index=True +) +duikersifonhevel_gdf = gpd.read_file( + cloud.joinpath(authority, "aangeleverd", "Aanlevering_202311", "HyDAMO_WM_20231117.gpkg"), + fid_as_index=True, + layer="duikersifonhevel", +) + +split_line_gdf = gpd.read_file( + cloud.joinpath(authority, "verwerkt", "fix_user_data.gpkg"), layer="split_basins", fid_as_index=True +) + +# level_boundary_gdf = gpd.read_file( +# cloud.joinpath(authority, "verwerkt", "fix_user_data.gpkg"), layer="level_boundary", fid_as_index=True +# ) + +# %% read model +model = Model.read(ribasim_toml) +ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_fix_model_network", f"{short_name}.toml") +network_validator = NetworkValidator(model) + +# %% some stuff we'll need again +manning_data = manning_resistance.Static(length=[100], manning_n=[0.04], profile_width=[10], profile_slope=[1]) +level_data = level_boundary.Static(level=[0]) + +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]), +] +outlet_data = outlet.Static(flow_rate=[100]) + + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2393424844 +# Verwijderen duplicate edges + +model.edge.df.drop_duplicates(inplace=True) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2393458802 + +# Toevoegen ontbrekende basins (oplossen topologie) +model.remove_node(7, remove_edges=True) +model.remove_node(84, remove_edges=True) +basin_edges_df = network_validator.edge_incorrect_connectivity() +basin_nodes_df = network_validator.node_invalid_connectivity() + +for row in basin_nodes_df.itertuples(): + # maak basin-node + basin_node = model.basin.add(Node(geometry=row.geometry), tables=basin_data) + + # update edge_table + model.edge.df.loc[basin_edges_df[basin_edges_df.from_node_id == row.node_id].index, ["from_node_id"]] = ( + basin_node.node_id + ) + model.edge.df.loc[basin_edges_df[basin_edges_df.to_node_id == row.node_id].index, ["to_node_id"]] = ( + basin_node.node_id + ) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2393672367 + +# Omdraaien edge-richting rondom outlets (inlaten/uitlaten) +# for edge_id in [2282, ] + +# https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2393731749 + +# Opruimen Reeve + +# basin 2484 wordt LevelBoundary (IJssel) +model.update_node(2484, "LevelBoundary", data=[level_data]) + +# nodes 1536, 762, 761, 1486 + aangesloten edges gooien we weg +for node_id in [1536, 762, 761, 1486]: + model.remove_node(node_id, remove_edges=True) + +# edges 2841, 2842, 2843, 2846 gooien we weg +model.remove_edges([2841, 2842, 2843, 2846]) + +# duiker 286309 voegen we toe +kdu = duikersifonhevel_gdf.loc[9063] +outlet_node = model.outlet.add( + Node(geometry=kdu.geometry.interpolate(0.5, normalized=True), name=f"duikersifonhevel.{kdu.objectid}"), + tables=[outlet_data], +) + +model.edge.add(model.level_boundary[10], outlet_node) +model.edge.add(outlet_node, model.basin[2240]) +model.edge.add(model.manning_resistance[849], model.basin[2240]) +model.edge.add(model.manning_resistance[760], model.basin[2240]) +model.edge.add(model.tabulated_rating_curve[187], model.basin[2240]) +model.edge.add(model.basin[2240], model.pump[1100]) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2393871075 + +# Ramsgeul bij Ramspol +for node_id in [81, 839]: + model.remove_node(node_id, remove_edges=True) + +model.update_node(83, "Basin", data=basin_data) +model.move_node(83, hydroobject_gdf.at[21045, "geometry"].boundary.geoms[0]) +model.reverse_edge(edge_id=3013) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399104407 + +# Deelstroomgebied Frysland + +# Nieuwe Kanaal / Tussen Linde voorzien van basin nabij gemaal.20 (node_id 701) +basin_node = model.basin.add(Node(geometry=hydroobject_gdf.at[19671, "geometry"].boundary.geoms[1]), tables=basin_data) +outlet_node = model.outlet.add( + Node(geometry=hydroobject_gdf.at[19671, "geometry"].interpolate(0.5, normalized=True)), + tables=[outlet_data], +) + +# basin 1623 verbinden met inlaatduikers (3x) en gemaal 1623; overige verbindingen verwijderen. +model.remove_edges([3038, 3040, 3037, 3041, 3039]) + +# nw basin verbinden met gemaal 20, level boundary 94 en alle inlaatduikers +model.reverse_edge(edge_id=2282) +model.edge.add(outlet_node, model.level_boundary[94]) +model.edge.add(basin_node, outlet_node) +model.edge.add(basin_node, model.manning_resistance[1182]) +model.edge.add(basin_node, model.manning_resistance[969]) +model.edge.add(basin_node, model.manning_resistance[1050]) +model.edge.add(basin_node, model.outlet[539]) +model.edge.add(model.pump[701], basin_node) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399209787 + +# Aansluiten NW boezem op Fryslan + +# basin /area 1681 op te knippen nabij basin 1717 (rode lijn) +model.split_basin(split_line_gdf.at[14, "geometry"]) +model.basin.area.df = model.basin.area.df[model.basin.area.df.node_id != 1717] + +# basin 1682 te veranderen in een LevelBoundary +model.update_node(1682, "LevelBoundary", [level_data]) + +# Alle edges die nu naar basin 1717 lopen naar LevelBoundary 1682 of opheffen +model.remove_node(27, remove_edges=True) +model.remove_node(1556, remove_edges=True) +model.remove_edges([945, 2537, 2536]) + +boundary_node = model.level_boundary.add(Node(geometry=hydroobject_gdf.at[7778, "geometry"].boundary.geoms[0])) + +model.edge.add(model.pump[642], boundary_node) +model.update_node(1202, "Outlet", data=[outlet_data]) +model.edge.add(boundary_node, model.outlet[1202]) +model.update_node(1203, "Outlet", data=[outlet_data]) +model.edge.add(boundary_node, model.outlet[1203]) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399328441 + +# Misc pump benedenstroomse edges +for edge_id in [2862, 3006, 3049]: + model.reverse_edge(edge_id=edge_id) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399355028 + +# Misc tabulated_rating_curve (stuwen) stroomrichting +for edge_id in [1884, 2197]: + model.reverse_edge(edge_id=edge_id) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399382478 + +# Misc manning_resistance (duikers) stroomrichting +for edge_id in [1081, 518]: + model.reverse_edge(edge_id=edge_id) + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399425885 +# Opknippen NW boezem + +poly1, poly2 = split_basin_multi_polygon(model.basin.area.df.at[598, "geometry"], split_line_gdf.at[15, "geometry"]) +model.basin.area.df.loc[model.basin.area.df.node_id == 1681, ["geometry"]] = poly1 + +poly1, poly2 = split_basin_multi_polygon(poly2, split_line_gdf.at[16, "geometry"]) +model.basin.area.df.loc[598] = {"node_id": 1686, "geometry": poly1} + +poly1, poly2 = split_basin_multi_polygon(poly2, split_line_gdf.at[17, "geometry"]) +model.basin.area.df.loc[model.basin.area.df.index.max() + 1] = {"geometry": poly1, "node_id": 1695} +model.basin.area.df.crs = model.crs + +tables = basin_data + [basin.Area(geometry=[poly2])] +basin_node = model.basin.add(Node(geometry=hydroobject_gdf.at[19608, "geometry"].boundary.geoms[1]), tables=tables) + + +model.move_node(1686, hydroobject_gdf.at[19566, "geometry"].boundary.geoms[1]) +model.merge_basins(basin_id=2426, to_basin_id=1696, are_connected=True) +model.merge_basins(basin_id=2460, to_basin_id=1696, are_connected=True) +model.merge_basins(basin_id=1648, to_basin_id=1696, are_connected=True) + +model.merge_basins(basin_id=1696, to_basin_id=2453, are_connected=True) + +model.merge_basins(basin_id=2453, to_basin_id=1686, are_connected=True) +model.merge_basins(basin_id=1719, to_basin_id=1686, are_connected=True) +model.merge_basins(basin_id=1858, to_basin_id=1686, are_connected=True) + +model.remove_node(1532, remove_edges=True) +model.remove_node(722, remove_edges=True) +model.remove_node(536, remove_edges=True) +model.remove_node(2506, remove_edges=True) + +edge_ids = [ + 2866, + 2867, + 2868, + 2869, + 2870, + 2871, + 2872, + 2873, + 2875, + 2876, + 2877, + 2878, + 2879, + 2880, + 2881, + 2883, + 2885, + 2886, + 2889, + 2890, + 2891, + 2895, + 2897, + 2899, + 2901, + 2902, + 2903, + 2905, + 2906, + 2907, + 2908, + 2910, + 2911, + 2912, + 2913, + 2915, + 2916, + 2918, +] + +for edge_id in edge_ids: + model.redirect_edge(edge_id, to_node_id=basin_node.node_id) + +model.remove_edges([2887, 2892]) +model.edge.add(basin_node, model.pump[547]) +model.edge.add(basin_node, model.outlet[540]) + +for edge_id in [2914, 2894]: + model.redirect_edge(edge_id, to_node_id=31) + +model.redirect_edge(461, to_node_id=1585) + +model.basin.area.df.loc[model.basin.area.df.node_id == 1545, ["node_id"]] = 1585 + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2400050483 + +# Ontbrekende basin beneden-Vecht +basin_node = model.basin.add(Node(geometry=hydroobject_gdf.at[12057, "geometry"].boundary.geoms[0]), tables=basin_data) +outlet_node = model.outlet.add( + Node(geometry=hydroobject_gdf.at[12057, "geometry"].interpolate(0.5, normalized=True)), tables=[outlet_data] +) + +for edge_id in [2956, 2957, 2958, 2959, 2960, 2961]: + model.redirect_edge(edge_id, to_node_id=basin_node.node_id) + +model.remove_node(76, remove_edges=True) +model.edge.add(basin_node, model.pump[598]) +model.edge.add(basin_node, outlet_node) +model.edge.add(outlet_node, model.level_boundary[50]) + + +# %% https://github.com/Deltares/Ribasim-NL/issues/147#issuecomment-2399931763 + +# Samenvoegen Westerveldse Aa +model.merge_basins(basin_id=1592, to_basin_id=1645, are_connected=True) +model.merge_basins(basin_id=1593, to_basin_id=1645, are_connected=True) + +model.merge_basins(basin_id=1645, to_basin_id=1585, are_connected=True) +model.merge_basins(basin_id=2567, to_basin_id=1585, are_connected=True) +model.merge_basins(basin_id=2303, to_basin_id=1585, are_connected=True) +model.merge_basins(basin_id=2549, to_basin_id=1585, are_connected=True) +model.merge_basins(basin_id=2568, to_basin_id=1585, are_connected=True) +model.merge_basins(basin_id=2572, to_basin_id=1585, are_connected=True) +model.merge_basins(basin_id=2374, to_basin_id=1585, are_connected=True) + +model.merge_basins(basin_id=2559, to_basin_id=2337, are_connected=False) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2382572457 + +# Administratie basin node_id in node_table en Basin / Area correct maken +# model.fix_unassigned_basin_area() +# model.fix_unassigned_basin_area(method="closest", distance=100) +# model.fix_unassigned_basin_area() + +# model.unassigned_basin_area.to_file("unassigned_basins.gpkg") +# model.basin.area.df = model.basin.area.df[~model.basin.area.df.node_id.isin(model.unassigned_basin_area.node_id)] + +# %% +# corrigeren knoop-topologie + +# ManningResistance bovenstrooms LevelBoundary naar Outlet +# for row in network_validator.edge_incorrect_type_connectivity().itertuples(): +# model.update_node(row.from_node_id, "Outlet", data=[outlet_data]) + +# # Inlaten van ManningResistance naar Outlet +# for row in network_validator.edge_incorrect_type_connectivity( +# from_node_type="LevelBoundary", to_node_type="ManningResistance" +# ).itertuples(): +# model.update_node(row.to_node_id, "Outlet", data=[outlet_data]) + + +# # buffer out small slivers +# model.basin.area.df.loc[:, ["geometry"]] = ( +# model.basin.area.df.buffer(0.1) +# .buffer(-0.1) +# .apply(lambda x: x if x.geom_type == "MultiPolygon" else MultiPolygon([x])) +# ) +# %% +# basin-profielen updaten + +df = pd.DataFrame( + { + "node_id": np.repeat(model.basin.node.df.index.to_numpy(), 2), + "level": [0.0, 1.0] * len(model.basin.node.df), + "area": [0.01, 1000.0] * len(model.basin.node.df), + } +) +df.index.name = "fid" +model.basin.profile.df = df + +df = model.basin.profile.df.groupby("node_id")[["level"]].max().reset_index() +df.index.name = "fid" +model.basin.state.df = df + +# %% +# tabulated_rating_curves updaten +df = pd.DataFrame( + { + "node_id": np.repeat(model.tabulated_rating_curve.node.df.index.to_numpy(), 2), + "level": [0.0, 5] * len(model.tabulated_rating_curve.node.df), + "flow_rate": [0, 0.1] * len(model.tabulated_rating_curve.node.df), + } +) +df.index.name = "fid" +model.tabulated_rating_curve.static.df = df + + +# %% + +# level_boundaries updaten +df = pd.DataFrame( + { + "node_id": model.level_boundary.node.df.index.to_list(), + "level": [0.0] * len(model.level_boundary.node.df), + } +) +df.index.name = "fid" +model.level_boundary.static.df = df + +# %% +# manning_resistance updaten +length = len(model.manning_resistance.node.df) +df = pd.DataFrame( + { + "node_id": model.manning_resistance.node.df.index.to_list(), + "length": [100.0] * length, + "manning_n": [100.0] * length, + "profile_width": [100.0] * length, + "profile_slope": [100.0] * length, + } +) +df.index.name = "fid" +model.manning_resistance.static.df = df + + +# %% write model +model.use_validation = False +model.write(ribasim_toml) + +# %% diff --git a/notebooks/vechtstromen/00_get_model.py b/notebooks/vechtstromen/00_get_model.py new file mode 100644 index 0000000..6964055 --- /dev/null +++ b/notebooks/vechtstromen/00_get_model.py @@ -0,0 +1,15 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +model_url = cloud.joinurl("Vechtstromen", "modellen", "Vechtstromen_2024_6_3") + +# %% +cloud.download_content(model_url) + + +# %% rename, so we can seperate in QGIS +ribasim_toml = cloud.joinpath("Vechtstromen", "modellen", "Vechtstromen_2024_6_3", "model.toml") +if ribasim_toml.exists(): + ribasim_toml.rename(ribasim_toml.with_name("vechtstromen.toml")) diff --git a/notebooks/vechtstromen/00_get_verwerkt.py b/notebooks/vechtstromen/00_get_verwerkt.py new file mode 100644 index 0000000..5cdc8b3 --- /dev/null +++ b/notebooks/vechtstromen/00_get_verwerkt.py @@ -0,0 +1,9 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +data_url = cloud.joinurl("DeDommel", "verwerkt") + +# %% +cloud.download_content(data_url) diff --git a/notebooks/vechtstromen/01_fix_model_network.py b/notebooks/vechtstromen/01_fix_model_network.py new file mode 100644 index 0000000..7d94912 --- /dev/null +++ b/notebooks/vechtstromen/01_fix_model_network.py @@ -0,0 +1,1028 @@ +# %% +import geopandas as gpd +import numpy as np +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 +from ribasim_nl.geometry import edge, split_basin, split_basin_multi_polygon +from shapely.geometry import LineString, MultiPolygon, Point, Polygon +from shapely.ops import nearest_points + +cloud = CloudStorage() + +ribasim_toml = cloud.joinpath("Vechtstromen", "modellen", "Vechtstromen_2024_6_3", "vechtstromen.toml") +database_gpkg = ribasim_toml.with_name("database.gpkg") +hydroobject_gdf = gpd.read_file( + cloud.joinpath("Vechtstromen", "verwerkt", "4_ribasim", "hydamo.gpkg"), layer="hydroobject", fid_as_index=True +) + +split_line_gdf = gpd.read_file( + cloud.joinpath("Vechtstromen", "verwerkt", "fix_user_data.gpkg"), layer="split_basins", fid_as_index=True +) + +level_boundary_gdf = gpd.read_file( + cloud.joinpath("Vechtstromen", "verwerkt", "fix_user_data.gpkg"), layer="level_boundary", fid_as_index=True +) + +# %% read model +model = Model.read(ribasim_toml) +ribasim_toml = cloud.joinpath("Vechtstromen", "modellen", "Vechtstromen_fix_model_network", "vechtstromen.toml") +network_validator = NetworkValidator(model) + +# %% some stuff we'll need again +manning_data = manning_resistance.Static(length=[100], manning_n=[0.04], profile_width=[10], profile_slope=[1]) +level_data = level_boundary.Static(level=[0]) + +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]), +] +outlet_data = outlet.Static(flow_rate=[100]) + + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2385111465 + +# Verwijderen duplicate edges + +model.edge.df.drop_duplicates(inplace=True) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2352686763 + +# Toevoegen benedenstroomse randvoorwaarden Beneden Dinkel + +# verander basin met node_id 2250 naar type level_boundary +model.update_node(2250, "LevelBoundary", data=[level_data]) + + +# verplaats basin 1375 naar het hydroobject +node_id = 1375 + +model.basin.node.df.loc[node_id, "geometry"] = hydroobject_gdf.at[3135, "geometry"].interpolate(0.5, normalized=True) +edge_ids = model.edge.df[ + (model.edge.df.from_node_id == node_id) | (model.edge.df.to_node_id == node_id) +].index.to_list() +model.reset_edge_geometry(edge_ids=edge_ids) + +# verplaats basin 1375 naar het hydroobject + + +# verbind basins met level_boundaries +for fid, node_id in [(1, 1375), (2, 1624)]: + boundary_node_geometry = level_boundary_gdf.at[fid, "geometry"] + + # line for interpolation + basin_node_geometry = Point( + model.basin.node.df.at[node_id, "geometry"].x, model.basin.node.df.at[node_id, "geometry"].y + ) + line_geometry = LineString((basin_node_geometry, boundary_node_geometry)) + + # define level_boundary_node + boundary_node = model.level_boundary.add(Node(geometry=boundary_node_geometry), tables=[level_data]) + level_node = model.level_boundary.add(Node(geometry=boundary_node_geometry), tables=[level_data]) + + # define manning_resistance_node + outlet_node_geometry = line_geometry.interpolate(line_geometry.length - 20) + outlet_node = model.outlet.add(Node(geometry=outlet_node_geometry), tables=[outlet_data]) + + from_node_id = model.basin[node_id].node_id + to_node_id = outlet_node.node_id + + # draw edges + # FIXME: we force edges to be z-less untill this is solved: https://github.com/Deltares/Ribasim/issues/1854 + model.edge.add( + model.basin[node_id], outlet_node, geometry=edge(model.basin[node_id].geometry, outlet_node.geometry) + ) + model.edge.add(outlet_node, boundary_node, geometry=edge(outlet_node.geometry, boundary_node.geometry)) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2382565944 + +# Verwijderen Twentekanaal (zit al bij RWS-HWS) +remove_node_ids = [1562, 1568, 1801, 1804, 1810, 1900, 2114, 2118, 2119, 32] + +# remove by edge so we also remove all resistance nodes in between +edge_df = model.edge.df[ + model.edge.df.from_node_id.isin(remove_node_ids) | model.edge.df.to_node_id.isin(remove_node_ids) +][["from_node_id", "to_node_id"]] + +for row in edge_df.itertuples(): + model.remove_edge(from_node_id=row.from_node_id, to_node_id=row.to_node_id, remove_disconnected_nodes=True) + +# add level_boundaries at twentekanaal for later coupling +hws_model = Model.read(cloud.joinpath("Rijkswaterstaat", "modellen", "hws", "hws.toml")) +basin_ids = hws_model.node_table().df[hws_model.node_table().df.name.str.contains("Twentekanaal")].index.to_list() +twentekanaal_poly = hws_model.basin.area.df[hws_model.basin.area.df.node_id.isin(basin_ids)].union_all() + +connect_node_ids = [ + i for i in set(edge_df[["from_node_id", "to_node_id"]].to_numpy().flatten()) if i in model._used_node_ids +] + +for node_id in connect_node_ids: + node = model.get_node(node_id=node_id) + + # update node to Outlet if it's a manning resistance + if node.node_type == "ManningResistance": + model.update_node(node.node_id, "Outlet", data=[outlet_data]) + node = model.get_node(node_id=node_id) + + _, boundary_node_geometry = nearest_points(node.geometry, twentekanaal_poly.boundary) + + boundary_node = model.level_boundary.add(Node(geometry=boundary_node_geometry), tables=[level_data]) + + # draw edge in the correct direction + if model.edge.df.from_node_id.isin([node_id]).any(): # supply + model.edge.add(boundary_node, node, geometry=edge(boundary_node.geometry, node.geometry)) + else: + model.edge.add(node, boundary_node, geometry=edge(node.geometry, boundary_node.geometry)) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2385525533 + +# Opruimen situatie rondom gemaal Oude Drostendiep + +# pumps met node_id 639, 608 en 603 op te heffen (1 gemaal ipv 3) +remove_node_ids = [639, 608, 603] + +for node_id in remove_node_ids: + model.remove_node(node_id, remove_edges=True) + +# remove by edge so we also remove all resistance nodes in between +edge_df = model.edge.df[ + model.edge.df.from_node_id.isin(remove_node_ids) | model.edge.df.to_node_id.isin(remove_node_ids) +][["from_node_id", "to_node_id"]] + +for row in edge_df.itertuples(): + model.remove_edge(from_node_id=row.from_node_id, to_node_id=row.to_node_id, remove_disconnected_nodes=True) + +# basin met node_id 1436 te verplaatsen naar locatie basin node_id 2259 +basin_id = 1436 +model.basin.node.df.loc[basin_id, "geometry"] = model.basin[2259].geometry +edge_ids = model.edge.df[ + (model.edge.df.from_node_id == basin_id) | (model.edge.df.to_node_id == basin_id) +].index.to_list() + +model.reset_edge_geometry(edge_ids=edge_ids) + +# basin met node_id 2259 opheffen (klein niets-zeggend bakje) +model.remove_node(2259, remove_edges=True) + +# stuw ST05005 (node_id 361) verbinden met basin met node_id 1436 +model.edge.add(model.tabulated_rating_curve[361], model.basin[basin_id]) +model.edge.add(model.basin[basin_id], model.pump[635]) + +# basin met node_id 2250 verplaatsen naar logische plek bovenstrooms ST05005 en bendenstrooms ST02886 op hydroobjec +basin_id = 2255 +model.basin.node.df.loc[basin_id, ["geometry"]] = hydroobject_gdf.at[6444, "geometry"].interpolate(0.5, normalized=True) + +edge_ids = model.edge.df[ + (model.edge.df.from_node_id == basin_id) | (model.edge.df.to_node_id == basin_id) +].index.to_list() + +model.reset_edge_geometry(edge_ids=edge_ids) + +model.split_basin(split_line_gdf.at[9, "geometry"]) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2385409772 + +incorrect_edges_df = network_validator.edge_incorrect_connectivity() +false_basin_ids = [1356, 1357, 1358, 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368, 1369, 1370] + +for false_basin_id in false_basin_ids: + basin_geom = ( + incorrect_edges_df[incorrect_edges_df.from_node_id == false_basin_id].iloc[0].geometry.boundary.geoms[0] + ) + basin_node = model.basin.add(Node(geometry=basin_geom), tables=basin_data) + + # fix edge topology + model.edge.df.loc[ + incorrect_edges_df[incorrect_edges_df.from_node_id == false_basin_id].index.to_list(), ["from_node_id"] + ] = basin_node.node_id + + model.edge.df.loc[ + incorrect_edges_df[incorrect_edges_df.to_node_id == false_basin_id].index.to_list(), ["to_node_id"] + ] = basin_node.node_id + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2386671759 + + +# basin 2224 en manning_resistance 898 (DK28491) opheffen +# tabulated_rating_cuves 336 (ST03745) en 238 (ST03744) opheffen +remove_node_ids = [2224, 898, 336, 238] + +for node_id in remove_node_ids: + model.remove_node(node_id, remove_edges=True) + +# pump 667 (GM00088) verbinden met basin 1495 +model.edge.add(model.pump[667], model.basin[1495]) + +# model.basin.area.df = model.basin.area.df[model.basin.area.df.node_id.isin(model.unassigned_basin_area.node_id)] + + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2387026622 + +# opruimen basin at Amsterdamscheveld +model.remove_node(1683, remove_edges=True) + +# verbinden basin node_id 1680 met tabulated_rating_curve node_id 101 en 125 +model.edge.add(model.basin[1680], model.tabulated_rating_curve[101]) +model.edge.add(model.basin[1680], model.tabulated_rating_curve[125]) + +# verbinden pump node_id 622 met basin node_id 1680 +model.edge.add(model.pump[622], model.basin[1680]) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2387056481 + +# Fix stieltjeskanaal + +# split basin_area bij manning_resistance node_id 1365 + +line = split_line_gdf.at[1, "geometry"] + +basin_polygon = model.basin.area.df.at[8, "geometry"].geoms[0] +basin_polygons = split_basin(basin_polygon, line) + +model.basin.area.df.loc[8, ["geometry"]] = MultiPolygon([basin_polygons.geoms[0]]) +model.basin.area.df.loc[model.basin.area.df.index.max() + 1, ["geometry"]] = MultiPolygon([basin_polygons.geoms[1]]) + +# hef basin node_id 1901 op +model.remove_node(1901, remove_edges=True) + +# hef pump node_id 574 (GM00246) en node_id 638 (GM00249) op +model.remove_node(574, remove_edges=True) +model.remove_node(638, remove_edges=True) + +# verbind basin node_id 1876 met pump node_ids 626 (GM00248) en 654 (GM00247) +model.edge.add(model.basin[1876], model.pump[626]) +model.edge.add(model.basin[1876], model.pump[654]) +model.edge.add(model.tabulated_rating_curve[113], model.basin[1876]) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2387815013 + +# opruimen Zwinderskanaal + + +# split basin_area bij rode lijn +line = split_line_gdf.at[2, "geometry"] +basin_polygon = model.basin.area.df.at[65, "geometry"].geoms[0] +basin_polygons = split_basin(basin_polygon, line) + +model.basin.area.df.loc[65, ["geometry"]] = MultiPolygon([basin_polygons.geoms[0]]) +model.basin.area.df.loc[model.basin.area.df.index.max() + 1, ["geometry"]] = MultiPolygon([basin_polygons.geoms[1]]) + +# verwijderen basin 2226, 2258, 2242 en manning_resistance 1366 +for node_id in [2226, 2258, 2242, 1366, 1350]: + model.remove_node(node_id, remove_edges=True) + +# verbinden tabulated_rating_curves 327 (ST03499) en 510 (ST03198) met basin 1897 +model.edge.add(model.basin[1897], model.tabulated_rating_curve[327]) +model.edge.add(model.tabulated_rating_curve[510], model.basin[1897]) + +# verbinden basin 1897 met tabulated_rating_curve 279 (ST03138) +model.edge.add(model.basin[1897], model.tabulated_rating_curve[279]) + +# verbinden basin 1897 met manning_resistance 1351 en 1352 +model.edge.add(model.basin[1897], model.manning_resistance[1351]) +model.edge.add(model.basin[1897], model.manning_resistance[1352]) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2387888742 + +# Oplossen situatie Van Echtenskanaal/Scholtenskanaal Klazinaveen + +# verplaatsen level_boundary 33 naar splitsing scholtenskanaal/echtenskanaal en omzetten naar basin +outlet_node_geometry = model.level_boundary[33].geometry +model.update_node(33, "Basin", data=basin_data) +model.move_node(node_id=33, geometry=hydroobject_gdf.loc[2679].geometry.boundary.geoms[1]) + +# plaatsen outlet bij oorspronkelijke plaats level_boundary 33 +outlet_node = model.outlet.add(Node(geometry=outlet_node_geometry), tables=[outlet_data]) + +# plaatsen nieuwe level_boundary op scholtenskanaal aan H&A zijde scholtenskanaal van outlet +boundary_node = model.level_boundary.add(Node(geometry=level_boundary_gdf.at[3, "geometry"]), tables=[level_data]) + +# toevoegen edges vanaf nieuwe basin 33 naar nieuwe outlet naar nieuwe boundary +model.edge.add(model.basin[33], outlet_node, geometry=edge(model.basin[33].geometry, outlet_node.geometry)) +model.edge.add(outlet_node, boundary_node) + +# opheffen manning_resistance 1330 bij GM00213 +model.remove_node(1330, remove_edges=True) + +# verbinden nieuwe basin met outlet en oorspronkijke manning_knopen en pompen in oorspronkelijke richting +for edge_id in [2711, 2712, 2713, 2714, 2708]: + model.reverse_edge(edge_id=edge_id) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2388009499 + +# basin met node_id 1873 gaat richting geesburg +model.move_node(node_id=1873, geometry=hydroobject_gdf.loc[6616].geometry.boundary.geoms[1]) +model.basin.area.df.loc[model.basin.area.df.node_id == 1873, ["node_id"]] = pd.NA +# ege 2700, 2701, 2702 worden opgeheven +model.edge.df = model.edge.df[~model.edge.df.index.isin([2700, 2701, 2702])] + +# basin 1873 wordt verbonden met manning_resistance 1054 +model.edge.add( + model.basin[1873], + model.manning_resistance[1054], + geometry=edge(model.basin[1873].geometry, model.manning_resistance[1054].geometry), +) + +# manning_resistance 1308 en 1331 worden verbonden met basin 1873 +model.edge.add( + model.manning_resistance[1308], + model.basin[1873], + geometry=edge(model.manning_resistance[1308].geometry, model.basin[1873].geometry), +) +model.edge.add( + model.basin[1873], + model.manning_resistance[1331], + geometry=edge(model.basin[1873].geometry, model.manning_resistance[1331].geometry), +) + +# level_boundary 26 wordt een outlet +model.update_node(26, "Outlet", data=[outlet_data]) + +# nieuwe level_boundary benedenstrooms nieuwe outlet 26 +boundary_node = model.level_boundary.add(Node(geometry=level_boundary_gdf.at[4, "geometry"]), tables=[level_data]) + + +# basin 1873 wordt verbonden met outlet en outlet met level_boundary +model.edge.add( + model.outlet[26], + model.basin[1873], + geometry=edge(model.outlet[26].geometry, model.basin[1873].geometry), +) + +model.edge.add(boundary_node, model.outlet[26]) + + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2388334544 + +# Kruising Dinkel/Kanaal Almelo Nordhorn corrigeren + +# ege 2700, 2701, 2702 worden opgeheven +model.edge.df = model.edge.df[~model.edge.df.index.isin([2690, 2691, 2692, 2693, 2694, 2695, 2696])] + +# basin / area splitten bij rode lijn in twee vlakken +line = split_line_gdf.at[3, "geometry"] + +total_basin_polygon = model.basin.area.df.at[544, "geometry"] +basin_polygon = [i for i in model.basin.area.df.at[544, "geometry"].geoms if i.intersects(line)][0] +basin_polygons = split_basin(basin_polygon, line) +model.basin.area.df.loc[544, ["geometry"]] = MultiPolygon( + [i for i in model.basin.area.df.at[544, "geometry"].geoms if not i.intersects(line)] + [basin_polygons.geoms[0]] +) +model.basin.area.df.loc[model.basin.area.df.index.max() + 1, ["geometry"]] = MultiPolygon([basin_polygons.geoms[1]]) + +# basin op dinkel bovenstrooms kanaal +dinkel_basin_node = model.basin.add( + Node(geometry=hydroobject_gdf.loc[2966].geometry.boundary.geoms[1]), tables=basin_data +) + +# basin in kanaal +kanaal_basin_node = model.basin.add( + Node(geometry=hydroobject_gdf.loc[7720].geometry.boundary.geoms[1]), tables=basin_data +) + +# edges v.a. tabulated_rating_curve 298 (ST01865) en 448 (ST01666) naar dinkel-basin +model.edge.add( + model.tabulated_rating_curve[298], + dinkel_basin_node, + geometry=edge(model.tabulated_rating_curve[298].geometry, dinkel_basin_node.geometry), +) + +model.edge.add( + model.tabulated_rating_curve[448], + dinkel_basin_node, + geometry=edge(model.tabulated_rating_curve[448].geometry, dinkel_basin_node.geometry), +) + +# edge v.a. manning_resistance 915 naar dinkel basin +model.edge.add( + model.manning_resistance[915], + dinkel_basin_node, + geometry=edge(model.manning_resistance[915].geometry, dinkel_basin_node.geometry), +) + +# edges v.a. dinkel basin naar tabulate_rating_curves 132 (ST02129) en 474 (ST02130) +model.edge.add( + dinkel_basin_node, + model.tabulated_rating_curve[132], + geometry=edge(dinkel_basin_node.geometry, model.tabulated_rating_curve[132].geometry), +) + +model.edge.add( + dinkel_basin_node, + model.tabulated_rating_curve[474], + geometry=edge(dinkel_basin_node.geometry, model.tabulated_rating_curve[474].geometry), +) + +# nieuwe manning_resistance in nieuwe dinkel-basin bovenstrooms kanaal +manning_node = model.manning_resistance.add( + Node(geometry=hydroobject_gdf.at[7721, "geometry"].interpolate(0.5, normalized=True)), tables=[manning_data] +) + +# nieuwe basin verbinden met nieuwe manning_resistance en nieuw kanaal basin +model.edge.add( + dinkel_basin_node, + manning_node, + geometry=edge(dinkel_basin_node.geometry, manning_node.geometry), +) + +model.edge.add( + manning_node, + kanaal_basin_node, + geometry=edge(manning_node.geometry, kanaal_basin_node.geometry), +) + +# nieuw kanaal-basin vervinden met tabulated_rating_curve 471 (ST01051) +model.edge.add( + kanaal_basin_node, + model.tabulated_rating_curve[471], + geometry=edge(kanaal_basin_node.geometry, model.tabulated_rating_curve[471].geometry), +) + +# nieuw kanaal-basin vervinden met manning_resistance 1346 +model.edge.add( + kanaal_basin_node, + model.manning_resistance[1346], + geometry=edge(kanaal_basin_node.geometry, model.manning_resistance[1346].geometry), +) + +# nieuwe outletlet bij grensduiker kanaal +outlet_node = model.outlet.add( + Node(geometry=hydroobject_gdf.at[7746, "geometry"].boundary.geoms[0]), tables=[outlet_data] +) + +# nieuwe basin verbinden met outlet verbinden met level_boundary 21 +model.edge.add( + outlet_node, + kanaal_basin_node, + geometry=edge(outlet_node.geometry, kanaal_basin_node.geometry), +) + +model.edge.add( + model.level_boundary[21], + outlet_node, + geometry=edge(model.level_boundary[21].geometry, outlet_node.geometry), +) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2389192454 +model.reverse_edge(edge_id=2685) +model.remove_node(node_id=2229, remove_edges=True) +model.edge.add( + model.basin[1778], + model.outlet[1080], + geometry=edge(model.basin[1778].geometry, model.outlet[1080].geometry), +) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2389198178 +model.reverse_edge(edge_id=2715) +model.reverse_edge(edge_id=2720) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2390712613 + +# Oplossen toplogische situatie kanaal Coevorden + +# opheffen basin 2243 en basin 2182 +model.remove_node(2243, remove_edges=True) +model.remove_node(2182, remove_edges=True) +model.remove_node(1351, remove_edges=True) +model.remove_node(1268, remove_edges=True) +model.remove_node(1265, remove_edges=True) + +# onknippen basin bij rode lijn +line = split_line_gdf.at[4, "geometry"] +basin_area_row = model.basin.area.df[model.basin.area.df.contains(line.centroid)].iloc[0] +basin_area_index = basin_area_row.name +basin_polygon = basin_area_row.geometry.geoms[0] +basin_polygons = split_basin(basin_polygon, line) + +model.basin.area.df.loc[basin_area_index, ["geometry"]] = MultiPolygon([basin_polygons.geoms[1]]) +model.basin.area.df.loc[model.basin.area.df.index.max() + 1, ["geometry"]] = MultiPolygon([basin_polygons.geoms[0]]) + +# # verplaatsen basin 1678 naar kruising waterlopen +model.move_node(node_id=1678, geometry=hydroobject_gdf.loc[6594].geometry.boundary.geoms[1]) + +# verwijderen edges 809, 814, 807, 810, 1293, 2772 +model.edge.df = model.edge.df[~model.edge.df.index.isin([809, 814, 807, 810, 887])] + + +# verbinden manning 1270, 1127 en pumps 644, 579 en 649 met basin 1678 +for node_id in [1270, 1127]: + model.edge.add( + model.manning_resistance[node_id], + model.basin[1678], + geometry=edge(model.manning_resistance[node_id].geometry, model.basin[1678].geometry), + ) + +for node_id in [644, 579, 649]: + model.edge.add( + model.pump[node_id], + model.basin[1678], + geometry=edge(model.pump[node_id].geometry, model.basin[1678].geometry), + ) + +# verplaatsen manning 1267 naar basin-edge tussen 1678 en 1678 +model.move_node(node_id=1267, geometry=hydroobject_gdf.loc[6609].geometry.boundary.geoms[1]) + +# maak nieuwe manning-node tussen 1678 en 1897 +manning_node = model.manning_resistance.add( + Node(geometry=hydroobject_gdf.loc[6596].geometry.interpolate(0.5, normalized=True)), tables=[manning_data] +) + +# verbinden basin 1897 met manning-node +model.edge.add( + model.basin[1897], + manning_node, + geometry=edge(model.basin[1897].geometry, manning_node.geometry), +) + +model.edge.add( + manning_node, + model.basin[1678], + geometry=edge(manning_node.geometry, model.basin[1678].geometry), +) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2390952469 + +# Schoonebekerdiep v.a. Twist Bült + +# verplaatsen basin 1909 nabij tabulated_rating_curve 383 (ST03607) +model.move_node(1909, geometry=hydroobject_gdf.loc[6865].geometry.boundary.geoms[1]) + +# verwijderen edges 780 en 778 +model.edge.df = model.edge.df[~model.edge.df.index.isin([780, 778])] + +# toevoegen edge tussen tabulated_rating_curve 383 en basin 1909 +model.edge.add( + model.tabulated_rating_curve[383], + model.basin[1909], + geometry=edge(model.tabulated_rating_curve[383].geometry, model.basin[1909].geometry), +) + +# toevoegen edge tussen manning_resistance 851 en basin 1909 +model.edge.add( + model.manning_resistance[851], + model.basin[1909], + geometry=edge(model.manning_resistance[851].geometry, model.basin[1909].geometry), +) + +# opknippen basin 1538 nabij 1909 en verbinden basin 1909 met 1539 via nieuwe manning_knoop +line = split_line_gdf.at[5, "geometry"] +model.split_basin(line=line) +manning_node = model.manning_resistance.add( + Node(geometry=line.intersection(hydroobject_gdf.at[6866, "geometry"])), tables=[manning_data] +) + +model.edge.add(model.basin[1909], manning_node) +model.edge.add(manning_node, model.basin[1539], geometry=edge(manning_node.geometry, model.basin[1539].geometry)) + +# verwijderen edge 2716,2718,2718,2719 +model.edge.df = model.edge.df[~model.edge.df.index.isin([2716, 2717, 2718, 2719])] + +# opknippen basin 2181 nabij 1881 en verbinden basin 1881 met 2181 via nieuwe manning_knoop +model.move_node(1881, geometry=hydroobject_gdf.loc[6919].geometry.boundary.geoms[1]) +line = split_line_gdf.at[6, "geometry"] +model.split_basin(line=line) +manning_node = model.manning_resistance.add( + Node(geometry=line.intersection(hydroobject_gdf.at[6879, "geometry"])), tables=[manning_data] +) + +model.edge.add(model.basin[1881], manning_node) +model.edge.add(manning_node, model.basin[2181], geometry=edge(manning_node.geometry, model.basin[2181].geometry)) + +for node_id in [139, 251, 267, 205]: + model.edge.add( + model.tabulated_rating_curve[node_id], + model.basin[1881], + geometry=edge(model.tabulated_rating_curve[node_id].geometry, model.basin[1881].geometry), + ) + +model.move_node(1269, hydroobject_gdf.at[7749, "geometry"].interpolate(0.5, normalized=True)) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391168839 + +# Molengoot Hardenberg + +# opheffen basin 1903 +model.remove_node(1903, remove_edges=True) + +# verbinden basin 1433 met pump 621 +model.edge.add( + model.basin[1433], + model.pump[621], + geometry=edge(model.basin[1433].geometry, model.pump[621].geometry), +) + +# verbinden tabulated_rating_curves 99 en 283 met basin 1433 +for node_id in [99, 283]: + model.edge.add( + model.tabulated_rating_curve[node_id], + model.basin[1433], + geometry=edge(model.tabulated_rating_curve[node_id].geometry, model.basin[1433].geometry), + ) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2390898004 + +model.remove_node(1131, remove_edges=True) +model.remove_node(1757, remove_edges=True) +model.edge.add( + model.basin[1588], + model.tabulated_rating_curve[112], + geometry=edge(model.basin[1588].geometry, model.tabulated_rating_curve[112].geometry), +) +model.edge.add( + model.basin[1588], + model.manning_resistance[57], + geometry=edge(model.basin[1588].geometry, model.manning_resistance[57].geometry), +) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391191673 + +# verwijderen basin 1905 +model.remove_node(1905, remove_edges=True) + +# verbinden manning_resistance 995 met basin 2148 +model.edge.add( + model.basin[2148], + model.manning_resistance[995], +) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391460899 + +# Samenvoegen basin-knopen Overijsselse Vecht & Coevorden Vechtkanaal +for basin_id in [1845, 2244, 2006, 1846]: + model.merge_basins(basin_id=basin_id, to_basin_id=2222) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391666745 + +# Opruimen basins Nieuw-Amsterdam + +# opknippen basin 1611 bij rode lijn, area mergen met basin 1879 +basin_polygon = model.basin.area.df.set_index("node_id").at[1611, "geometry"] +left_poly, right_poly = split_basin_multi_polygon(basin_polygon, split_line_gdf.at[8, "geometry"]) +model.basin.area.df.loc[model.basin.area.df.node_id == 1611, ["geometry"]] = right_poly + +left_poly = model.basin.area.df.set_index("node_id").at[1879, "geometry"].union(left_poly) +model.basin.area.df.loc[model.basin.area.df.node_id == 1879, ["geometry"]] = MultiPolygon([left_poly]) + +# merge basins 2186, 2173, 2022, 1611, 2185 in basin 1902 +for basin_id in [2186, 2173, 2022, 1611, 2185]: + model.merge_basins(basin_id=basin_id, to_basin_id=1902) + +# verplaats 1902 iets bovenstrooms +model.move_node(1902, hydroobject_gdf.at[6615, "geometry"].boundary.geoms[1]) + + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391686269 +model.remove_node(2198, remove_edges=True) +model.remove_node(2200, remove_edges=True) +model.edge.add(model.basin[2111], model.pump[671]) +model.edge.add(model.tabulated_rating_curve[542], model.basin[2111]) +model.edge.add(model.pump[671], model.basin[2316]) +model.edge.add(model.basin[2316], model.tabulated_rating_curve[542]) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391710413 + +model.remove_node(2202, remove_edges=True) +model.edge.add(model.basin[1590], model.pump[657]) +model.edge.add(model.manning_resistance[1058], model.basin[1590]) + + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391672700 + +# Merge basin 2176 in 1605 +model.merge_basins(basin_id=2176, to_basin_id=1605) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391726774 +# Merge basins 2206 in 1518 +model.merge_basins(basin_id=2206, to_basin_id=1518, are_connected=False) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391734144 +# dood takje uit Overijsselse Vecht +model.remove_node(2210, remove_edges=True) +model.remove_node(1294, remove_edges=True) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391740603 + +# Merge basin 2225 met 2304 +model.merge_basins(2225, 2304) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391815016 + +# Wetteringe als laterale inflow +model.merge_basins(2231, 1853) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391750536 + +# Rondom SL00010 opruimen +model.remove_node(2230, remove_edges=True) +model.remove_node(2251, remove_edges=True) +model.edge.add(model.outlet[41], model.level_boundary[15]) +model.edge.add(model.basin[1442], model.pump[664]) +model.edge.add(model.basin[1442], model.pump[665]) + + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391820198 + +# Merge basin 2232 in 1591 +model.merge_basins(basin_id=2232, to_basin_id=1591) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391825301 + +# Basin 2236 naar LevelBoundary +model.update_node(2236, "LevelBoundary", data=[level_data]) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391829471 + +# Merge basin 2246 en 1419 +model.merge_basins(basin_id=2246, to_basin_id=1419) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391946915 + +# Opruimen Elsbeek + +# Basin 2256 verplaatsen naar punt +model.move_node(node_id=2256, geometry=hydroobject_gdf.loc[2896].geometry.boundary.geoms[1]) + +# Basin knippen over lijn +model.basin.area.df.loc[model.basin.area.df.node_id == 2256, ["node_id"]] = pd.NA +model.split_basin(split_line_gdf.at[10, "geometry"]) + +# Edges 446, 1516, 443 en 444 verwijderen +model.edge.df = model.edge.df[~model.edge.df.index.isin([446, 1516, 443, 444])] + +# tabulated_rating_curves 202 en 230 verbinden met basin 2256 +model.edge.add( + model.tabulated_rating_curve[202], + model.basin[2256], + geometry=edge(model.tabulated_rating_curve[202].geometry, model.basin[2256].geometry), +) +model.edge.add( + model.tabulated_rating_curve[230], + model.basin[2256], + geometry=edge(model.tabulated_rating_curve[230].geometry, model.basin[2256].geometry), +) + +# resistance 954 verbinden met basin 2256 +model.edge.add( + model.manning_resistance[954], + model.basin[2256], + geometry=edge(model.manning_resistance[954].geometry, model.basin[2256].geometry), +) + +# basin 2256 verbinden met resistance 1106 +model.edge.add( + model.basin[2256], + model.manning_resistance[1106], + geometry=edge(model.basin[2256].geometry, model.manning_resistance[1106].geometry), +) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391984234 + +# Merge basin 2261 in basin 1698 +model.merge_basins(2261, 1698) +model.remove_node(390, remove_edges=True) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2391995841 + +# Merge basin 2260 met basin 1645 +model.merge_basins(2260, 1645) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392010526 +# Merge basin 2220 met basin 1371 +model.merge_basins(2220, 1371, are_connected=False) + + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392017041 + +# Kanaal Almelo Nordhorn bij Almelo +model.merge_basins(2219, 1583, are_connected=False) +model.merge_basins(2209, 1583, are_connected=False) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392022887 + +# Merge basin 2203 met 2227 +model.merge_basins(2203, 2227, are_connected=False) +model.remove_node(1219, remove_edges=True) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392026739 + +# Merge basin 2014 met 2144 +model.merge_basins(2014, 2144) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392030268 + +# Merge basin 1696 met 1411 +model.merge_basins(1696, 1411) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392037263 + +# Merge basin 2264 met 1459 +model.merge_basins(2264, 1459) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392043973 + +# Merge basin 2212 en 2310 +model.merge_basins(2212, 2310) +poly = model.basin.area.df.at[59, "geometry"].union(model.basin.area.df.set_index("node_id").at[2310, "geometry"]) +model.basin.area.df.loc[model.basin.area.df.node_id == 2310, ["geometry"]] = MultiPolygon([poly]) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392048684 + +# Merge basin 2253 in basin 2228 +model.merge_basins(2253, 2228) + + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392052379 + +# Merge basin 2221 in basin 1634 +model.merge_basins(2221, 1634) + +# %% https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2392076634 + +# Verbinding rondwaterleiding / Lennelwaterleiding herstellen +model.merge_basins(1859, 2235, are_connected=False) + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2382572457 + +# Administratie basin node_id in node_table en Basin / Area correct maken +model.fix_unassigned_basin_area() +model.fix_unassigned_basin_area(method="closest", distance=100) +model.fix_unassigned_basin_area() + +model.unassigned_basin_area.to_file("unassigned_basins.gpkg") +model.basin.area.df = model.basin.area.df[~model.basin.area.df.node_id.isin(model.unassigned_basin_area.node_id)] + +# %% +# corrigeren knoop-topologie + +# ManningResistance bovenstrooms LevelBoundary naar Outlet +for row in network_validator.edge_incorrect_type_connectivity().itertuples(): + model.update_node(row.from_node_id, "Outlet", data=[outlet_data]) + +# Inlaten van ManningResistance naar Outlet +for row in network_validator.edge_incorrect_type_connectivity( + from_node_type="LevelBoundary", to_node_type="ManningResistance" +).itertuples(): + model.update_node(row.to_node_id, "Outlet", data=[outlet_data]) + + +# %% see: https://github.com/Deltares/Ribasim-NL/issues/146#issuecomment-2382578661 + +# opvullen gaten +basin_polygon = model.basin.area.df.union_all() +# holes = [Polygon(interior) for polygon in basin_polygon.buffer(10).buffer(-10).geoms for interior in polygon.interiors] +holes = [Polygon(interior) for interior in basin_polygon.buffer(10).buffer(-10).interiors] +holes_df = gpd.GeoSeries(holes, crs=28992) +holes_df.index = holes_df.index + 1 +holes_df.to_file( + "holes.gpkg", + index=True, + fid="fid", +) +# splitsen Alemelo-Nordhorn / Overijsselskanaal. Overijsselskanaal zit in HWS +line = split_line_gdf.at[12, "geometry"] +idx = holes_df[holes_df.intersects(line)].index[0] +poly = split_basin(holes_df[holes_df.intersects(line)].iloc[0], line).geoms[0] +poly = model.basin.area.df.set_index("node_id").at[1583, "geometry"].union(poly) +model.basin.area.df.loc[model.basin.area.df.node_id == 1583, ["geometry"]] = MultiPolygon([poly]) + +# Split Overijsselskanaal bij Zwolsekanaal +line = split_line_gdf.at[13, "geometry"] +poly1, poly2 = split_basin(holes_df[holes_df.intersects(line)].iloc[0], line).geoms + +poly1 = model.basin.area.df.set_index("node_id").at[2116, "geometry"].union(poly1) +poly1 = MultiPolygon([i for i in poly1.geoms if i.geom_type == "Polygon"]) +model.basin.area.df.loc[model.basin.area.df.node_id == 2116, ["geometry"]] = poly1 + +poly2 = model.basin.area.df.set_index("node_id").at[2115, "geometry"].union(poly2) +poly2 = MultiPolygon([i for i in poly2.geoms if i.geom_type == "Polygon"] + [holes_df.loc[38], holes_df.loc[29]]) +model.basin.area.df.loc[model.basin.area.df.node_id == 2115, ["geometry"]] = poly2 + +# de rest gaan we automatisch vullen +holes_df = holes_df[~holes_df.index.isin([10, 22, 29, 32, 38, 39, 41])] + +holes_df.to_file( + "holes.gpkg", + index=True, + fid="fid", +) + +drainage_areas_df = gpd.read_file( + cloud.joinpath("Vechtstromen", "verwerkt", "4_ribasim", "areas.gpkg"), layer="drainage_areas" +) + +drainage_areas_df = drainage_areas_df[drainage_areas_df.buffer(-10).intersects(basin_polygon)] + +for idx, geometry in enumerate(holes_df): + # select drainage-area + drainage_area_select = drainage_areas_df[drainage_areas_df.contains(geometry.buffer(-10))] + if not drainage_area_select.empty: + if not len(drainage_area_select) == 1: + raise ValueError("hole contained by multiple drainage areas, can't fix that yet") + + drainage_area = drainage_area_select.iloc[0].geometry + + # find basin_id to merge to + selected_basins_df = model.basin.area.df[ + model.basin.area.df.node_id.isin(model.basin.node.df[model.basin.node.df.within(drainage_area)].index) + ].set_index("node_id") + if selected_basins_df.empty: + selected_basins_df = model.basin.area.df[ + model.basin.area.df.buffer(-10).intersects(drainage_area) + ].set_index("node_id") + + assigned_basin_id = selected_basins_df.intersection(geometry.buffer(10)).area.idxmax() + + # clip and merge geometry + geometry = geometry.buffer(10).difference(basin_polygon) + geometry = ( + model.basin.area.df.set_index("node_id") + .at[assigned_basin_id, "geometry"] + .union(geometry) + .buffer(0.1) + .buffer(-0.1) + ) + + if isinstance(geometry, Polygon): + geometry = MultiPolygon([geometry]) + model.basin.area.df.loc[model.basin.area.df.node_id == assigned_basin_id, "geometry"] = geometry + +# buffer out small slivers +model.basin.area.df.loc[:, ["geometry"]] = ( + model.basin.area.df.buffer(0.1) + .buffer(-0.1) + .apply(lambda x: x if x.geom_type == "MultiPolygon" else MultiPolygon([x])) +) +# %% +# basin-profielen updaten + +df = pd.DataFrame( + { + "node_id": np.repeat(model.basin.node.df.index.to_numpy(), 2), + "level": [0.0, 1.0] * len(model.basin.node.df), + "area": [0.01, 1000.0] * len(model.basin.node.df), + } +) +df.index.name = "fid" +model.basin.profile.df = df + +df = model.basin.profile.df.groupby("node_id")[["level"]].max().reset_index() +df.index.name = "fid" +model.basin.state.df = df + +# %% +# tabulated_rating_curves updaten +df = pd.DataFrame( + { + "node_id": np.repeat(model.tabulated_rating_curve.node.df.index.to_numpy(), 2), + "level": [0.0, 5] * len(model.tabulated_rating_curve.node.df), + "flow_rate": [0, 0.1] * len(model.tabulated_rating_curve.node.df), + } +) +df.index.name = "fid" +model.tabulated_rating_curve.static.df = df + + +# %% + +# level_boundaries updaten +df = pd.DataFrame( + { + "node_id": model.level_boundary.node.df.index.to_list(), + "level": [0.0] * len(model.level_boundary.node.df), + } +) +df.index.name = "fid" +model.level_boundary.static.df = df + +# %% +# manning_resistance updaten +length = len(model.manning_resistance.node.df) +df = pd.DataFrame( + { + "node_id": model.manning_resistance.node.df.index.to_list(), + "length": [100.0] * length, + "manning_n": [100.0] * length, + "profile_width": [100.0] * length, + "profile_slope": [100.0] * length, + } +) +df.index.name = "fid" +model.manning_resistance.static.df = df + +# %% write model + +model.basin.area.df.loc[:, ["meta_area"]] = model.basin.area.df.area +model.basin.node.df[~model.basin.node.df.index.isin(model.basin.area.df.node_id)].to_file("missing_areas.gpkg") + + +# model.use_validation = True +model.write(ribasim_toml) + +# %% diff --git a/notebooks/vechtstromen/99_upload_model.py b/notebooks/vechtstromen/99_upload_model.py new file mode 100644 index 0000000..31aa38b --- /dev/null +++ b/notebooks/vechtstromen/99_upload_model.py @@ -0,0 +1,6 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +cloud.upload_model("Vechtstromen", "Vechtstromen", include_results=False, include_plots=False) diff --git a/src/ribasim_nl/ribasim_nl/geometry.py b/src/ribasim_nl/ribasim_nl/geometry.py index 379333d..1457da7 100644 --- a/src/ribasim_nl/ribasim_nl/geometry.py +++ b/src/ribasim_nl/ribasim_nl/geometry.py @@ -104,6 +104,37 @@ def split_basin(basin_polygon: Polygon, line: LineString) -> MultiPolygon: return MultiPolygon(sort_basins(keep_polys)) +def split_basin_multi_polygon(basin_polygon: MultiPolygon, line: LineString) -> tuple[MultiPolygon, MultiPolygon]: + """Split a MultiPolygon into two given a LineString.""" + line_centre = line.interpolate(0.5, normalized=True) + + # get the polygon to cut + basin_geoms = list(basin_polygon.geoms) + if len(basin_geoms) == 0: + cut_idx = 0 + else: + try: + cut_idx = next(idx for idx, i in enumerate(basin_geoms) if i.contains(line_centre)) + except StopIteration: + cut_idx = next(idx for idx, i in enumerate(basin_geoms) if line.intersects(i)) + + # split it + right_basin_poly, left_basin_poly = split_basin(basin_geoms[cut_idx], line).geoms + + # concat left-over polygons to the right-side + right_basin_poly = [right_basin_poly] + left_basin_poly = [left_basin_poly] + + for idx, geom in enumerate(basin_geoms): + if idx != cut_idx: + if geom.distance(right_basin_poly[0]) < geom.distance(left_basin_poly[0]): + right_basin_poly += [geom] + else: + left_basin_poly += [geom] + + return MultiPolygon(right_basin_poly), MultiPolygon(left_basin_poly) + + def drop_z(geometry: LineString | MultiPolygon | Point | Polygon) -> Point | Polygon | MultiPolygon: """Drop the z-coordinate of a geometry if it has. @@ -145,3 +176,21 @@ def drop_z(geometry: LineString | MultiPolygon | Point | Polygon) -> Point | Pol ) return geometry + + +def edge(point_from: Point, point_to: Point) -> LineString: + """Create a LineString geometry between two Point geometries, dropping z-coordinate if any + + Args: + point_from (Point): _description_ + point_to (Point): _description_ + + Returns + ------- + LineString: LineString without z-coordinate + """ + if point_from.has_z: + point_from = Point(point_from.x, point_from.y) + if point_to.has_z: + point_to = Point(point_to.x, point_to.y) + return LineString((point_from, point_to)) diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index f8a4b51..0a9a8f6 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -1,14 +1,16 @@ from pathlib import Path from typing import Literal +import networkx as nx import pandas as pd from pydantic import BaseModel from ribasim import Model, Node from ribasim.geometry.edge import NodeData -from shapely.geometry import LineString, Point +from shapely.geometry import LineString, MultiPolygon, Point, Polygon from shapely.geometry.base import BaseGeometry from ribasim_nl.case_conversions import pascal_to_snake_case +from ribasim_nl.geometry import split_basin def read_arrow(filepath: Path) -> pd.DataFrame: @@ -46,6 +48,10 @@ def basin_results(self): self._basin_results = BasinResults(filepath=filepath) return self._basin_results + @property + def graph(self): + return nx.from_pandas_edgelist(self.edge.df[["from_node_id", "to_node_id"]], "from_node_id", "to_node_id") + @property def next_node_id(self): return self.node_table().df.index.max() + 1 @@ -260,15 +266,18 @@ 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): + def reverse_edge(self, from_node_id: int | None = None, to_node_id: int | None = None, edge_id: int | None = None): """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"] + if edge_id is 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"] + else: + edge_data = dict(self.edge.df.loc[edge_id]) # revert node ids self.edge.df.loc[edge_id, ["from_node_id"]] = edge_data["to_node_id"] @@ -294,6 +303,25 @@ def remove_edge(self, from_node_id: int, to_node_id: int, remove_disconnected_no if node_id not in self.edge.df[["from_node_id", "to_node_id"]].to_numpy().ravel(): self.remove_node(node_id) + def remove_edges(self, edge_ids: list[int]): + if self.edge.df is not None: + self.edge.df = self.edge.df[~self.edge.df.index.isin(edge_ids)] + + def move_node(self, node_id: int, geometry: Point): + node_type = self.node_table().df.at[node_id, "node_type"] + + # read existing table + table = getattr(self, pascal_to_snake_case(node_type)) + + # update geometry + table.node.df.loc[node_id, ["geometry"]] = geometry + + # reset all edges + edge_ids = self.edge.df[ + (self.edge.df.from_node_id == node_id) | (self.edge.df.to_node_id == node_id) + ].index.to_list() + self.reset_edge_geometry(edge_ids=edge_ids) + def find_closest_basin(self, geometry: BaseGeometry, max_distance: float | None) -> NodeData: """Find the closest basin_node.""" # only works when basin area are defined @@ -361,7 +389,9 @@ def reset_edge_geometry(self, edge_ids: list | None = None): df = self.edge.df for row in df.itertuples(): - geometry = LineString([node_df.at[row.from_node_id, "geometry"], node_df.at[row.to_node_id, "geometry"]]) + from_point = Point(node_df.at[row.from_node_id, "geometry"].x, node_df.at[row.from_node_id, "geometry"].y) + to_point = Point(node_df.at[row.to_node_id, "geometry"].x, node_df.at[row.to_node_id, "geometry"].y) + geometry = LineString([from_point, to_point]) self.edge.df.loc[row.Index, ["geometry"]] = geometry @property @@ -373,3 +403,116 @@ def edge_from_node_type(self): def edge_to_node_type(self): node_df = self.node_table().df return self.edge.df.to_node_id.apply(lambda x: node_df.at[x, "node_type"] if x in node_df.index else None) + + def split_basin(self, line: LineString): + if self.basin.area.df is None: + raise ValueError("provide basin / area table first") + + line_centre = line.interpolate(0.5, normalized=True) + basin_area_df = self.basin.area.df[self.basin.area.df.contains(line_centre)] + + if len(basin_area_df) != 1: + raise ValueError("Overlapping basin-areas at cut_line") + + # get all we need + basin_fid = int(basin_area_df.iloc[0].name) + basin_geometry = basin_area_df.iloc[0].geometry + self.basin.area.df = self.basin.area.df[self.basin.area.df.index != basin_fid] + + # get the polygon to cut + basin_geoms = list(basin_geometry.geoms) + cut_idx = next(idx for idx, i in enumerate(basin_geoms) if i.contains(line_centre)) + + # split it + right_basin_poly, left_basin_poly = split_basin(basin_geoms[cut_idx], line).geoms + + # concat left-over polygons to the right-side + right_basin_poly = [right_basin_poly] + left_basin_poly = [left_basin_poly] + + for idx, geom in enumerate(basin_geoms): + if idx != cut_idx: + if geom.distance(right_basin_poly[0]) < geom.distance(left_basin_poly[0]): + right_basin_poly += [geom] + else: + left_basin_poly += [geom] + + right_basin_poly = MultiPolygon(right_basin_poly) + left_basin_poly = MultiPolygon(left_basin_poly) + + for poly in [right_basin_poly, left_basin_poly]: + if self.basin.node.df.geometry.within(poly).any(): + node_ids = self.basin.node.df[self.basin.node.df.geometry.within(poly)].index.to_list() + if len(node_ids) == 1: + if node_ids[0] in self.basin.area.df.node_id.to_numpy(): + self.basin.area.df.loc[self.basin.area.df.node_id.isin(node_ids), ["geometry"]] = poly + else: + self.basin.area.df.loc[self.basin.area.df.index.max() + 1] = { + "node_id": node_ids[0], + "geometry": poly, + } + else: + self.basin.area.df.loc[self.basin.area.df.index.max() + 1, ["geometry"]] = poly + else: + self.basin.area.df.loc[self.basin.area.df.index.max() + 1, ["geometry"]] = poly + + if self.basin.area.df.crs is None: + self.basin.area.df.crs = self.crs + + def redirect_edge(self, edge_id: int, from_node_id: int | None = None, to_node_id: int | None = None): + if self.edge.df is not None: + if from_node_id is not None: + self.edge.df.loc[edge_id, ["from_node_id"]] = from_node_id + if to_node_id is not None: + self.edge.df.loc[edge_id, ["to_node_id"]] = to_node_id + + self.reset_edge_geometry(edge_ids=[edge_id]) + + def merge_basins(self, basin_id: int, to_basin_id: int, are_connected=True): + for node_id in (basin_id, to_basin_id): + if node_id not in self.basin.node.df.index: + raise ValueError(f"{node_id} is not a basin") + + if are_connected: + paths = [i for i in nx.all_shortest_paths(self.graph, basin_id, to_basin_id) if len(i) == 3] + + if len(paths) == 0: + raise ValueError(f"basin {basin_id} not a direct neighbor of basin {to_basin_id}") + + # remove flow-node and connected edges + for path in paths: + self.remove_node(path[1], remove_edges=True) + + # get a complete edge-list to modify + edge_ids = self.edge.df[self.edge.df.from_node_id == basin_id].index.to_list() + edge_ids += self.edge.df[self.edge.df.to_node_id == basin_id].index.to_list() + + # correct edge from and to attributes + self.edge.df.loc[self.edge.df.from_node_id == basin_id, "from_node_id"] = to_basin_id + self.edge.df.loc[self.edge.df.to_node_id == basin_id, "to_node_id"] = to_basin_id + + # reset edge geometries + self.reset_edge_geometry(edge_ids=edge_ids) + + # merge area if basin has any assigned to it + if basin_id in self.basin.area.df.node_id.to_numpy(): + poly = self.basin.area.df.set_index("node_id").at[basin_id, "geometry"] + + # if to_basin_id has area we union both areas + if to_basin_id in self.basin.area.df.node_id.to_numpy(): + poly = poly.union(self.basin.area.df.set_index("node_id").at[to_basin_id, "geometry"]) + if isinstance(poly, Polygon): + poly = MultiPolygon([poly]) + self.basin.area.df.loc[self.basin.area.df.node_id == to_basin_id, ["geometry"]] = poly + + # else we add a record to basin + else: + if isinstance(poly, Polygon): + poly = MultiPolygon([poly]) + self.basin.area.df.loc[self.basin.area.df.index.max() + 1] = {"node_id": to_basin_id, "geometry": poly} + + # finally we remove the basin + self.remove_node(basin_id) + + if self.basin.area.df.crs is None: + self.basin.area.df.crs = self.crs diff --git a/src/ribasim_nl/ribasim_nl/network_validator.py b/src/ribasim_nl/ribasim_nl/network_validator.py index 72472c5..13cb07d 100644 --- a/src/ribasim_nl/ribasim_nl/network_validator.py +++ b/src/ribasim_nl/ribasim_nl/network_validator.py @@ -1,5 +1,6 @@ from dataclasses import dataclass +import geopandas as gpd from ribasim import Model @@ -68,6 +69,31 @@ def node_internal_basin(self): mask = self.node_df.apply(lambda row: check_internal_basin(row, self.edge_df), axis=1) return self.node_df[mask] + def node_invalid_connectivity(self, tolerance: float = 1.0) -> gpd.GeoDataFrame: + """Check if node_from and node_to are correct on edge""" + node_df = self.node_df + invalid_edges_df = self.edge_incorrect_connectivity() + invalid_nodes = [] + for row in invalid_edges_df.itertuples(): + geoms = row.geometry.boundary.geoms + + for idx, attr in ((0, "from_node_id"), (1, "to_node_id")): + node_id = getattr(row, attr) + point = geoms[idx] + if node_id in node_df.index: + if point.distance(node_df.at[node_id, "geometry"]) > tolerance: + invalid_nodes += [{"node_id": node_id, "geometry": point}] + else: + invalid_nodes += [{"node_id": node_id, "geometry": point}] + + if invalid_nodes: + df = gpd.GeoDataFrame(invalid_nodes, crs=node_df.crs) + df.drop_duplicates(inplace=True) + else: + df = gpd.GeoDataFrame({"node_id": []}, geometry=gpd.GeoSeries(), crs=node_df.crs) + + return df + 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)]