From 31b72df73730360645ee6a1b846913b8fc97e399 Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Mon, 2 Dec 2024 16:54:55 +0100 Subject: [PATCH] Fix vlakken (#194) updates maken alle waterschapsmodellen rekenend. Deze branch nog niet weggooien; Basin / Area moet nog op orde gemaakt worden voor NZV, WDOD en HeA: image --------- Co-authored-by: ngoorden <49673619+ngoorden@users.noreply.github.com> --- notebooks/aa_en_maas/01_fix_model_network.py | 12 +- notebooks/aa_en_maas/01b_fix_basin_area.py | 243 ++++++++++++++++ .../amstel_gooi_en_vecht/00_get_model.py | 19 ++ notebooks/basin_area_nodes.py | 51 +++- .../brabantse_delta/01b_get_basin_holes.py | 24 ++ .../brabantse_delta/01c_clip_basin_holes.py | 46 +++ .../brabantse_delta/01d_fix_basin_area.py | 99 +++++++ notebooks/de_dommel/01_fix_model_network.py | 8 +- notebooks/de_dommel/01b_fix_basin_area.py | 53 ++++ notebooks/de_dommel/rare_edges.gpkg | Bin 106496 -> 0 bytes .../01_fix_model_network.py | 65 +++-- .../hunze_en_aas/01_fix_model_network.py | 204 ++++++++++++++ notebooks/koppelen/01_bergend_gebied.py | 1 + notebooks/limburg/01_fix_model_network.py | 4 +- notebooks/limburg/01b_fix_basin_area.py | 121 ++++++++ .../noorderzijlvest/01_fix_model_network.py | 3 +- .../noorderzijlvest/01b_fix_basin_area.py | 71 +++++ .../rijn_en_ijssel/01_fix_model_network.py | 3 + .../rijn_en_ijssel/01b_fix_basin_area.py | 102 +++++++ .../rijn_en_ijssel/01c_check_basin_holes.py | 26 ++ .../01_fix_model_network.py | 2 + .../stichtse_rijnlanden/01b_fix_basin_area.py | 58 ++++ .../vallei_en_veluwe/01_fix_model_network.py | 2 + .../vallei_en_veluwe/01b_fix_basin_area.py | 71 +++++ .../vechtstromen/01_fix_model_network.py | 3 +- notebooks/vechtstromen/01b_fix_basin_area.py | 58 ++++ src/ribasim_nl/ribasim_nl/model.py | 264 ++++++++++++++++-- .../ribasim_nl/network_validator.py | 10 +- 28 files changed, 1554 insertions(+), 69 deletions(-) create mode 100644 notebooks/aa_en_maas/01b_fix_basin_area.py create mode 100644 notebooks/amstel_gooi_en_vecht/00_get_model.py create mode 100644 notebooks/brabantse_delta/01b_get_basin_holes.py create mode 100644 notebooks/brabantse_delta/01c_clip_basin_holes.py create mode 100644 notebooks/brabantse_delta/01d_fix_basin_area.py create mode 100644 notebooks/de_dommel/01b_fix_basin_area.py delete mode 100644 notebooks/de_dommel/rare_edges.gpkg create mode 100644 notebooks/hunze_en_aas/01_fix_model_network.py create mode 100644 notebooks/limburg/01b_fix_basin_area.py create mode 100644 notebooks/noorderzijlvest/01b_fix_basin_area.py create mode 100644 notebooks/rijn_en_ijssel/01b_fix_basin_area.py create mode 100644 notebooks/rijn_en_ijssel/01c_check_basin_holes.py create mode 100644 notebooks/stichtse_rijnlanden/01b_fix_basin_area.py create mode 100644 notebooks/vallei_en_veluwe/01b_fix_basin_area.py create mode 100644 notebooks/vechtstromen/01b_fix_basin_area.py diff --git a/notebooks/aa_en_maas/01_fix_model_network.py b/notebooks/aa_en_maas/01_fix_model_network.py index 5d3307c..83db91e 100644 --- a/notebooks/aa_en_maas/01_fix_model_network.py +++ b/notebooks/aa_en_maas/01_fix_model_network.py @@ -114,10 +114,12 @@ model.remove_node(node_id, remove_edges=True) +# %% see: https://github.com/Deltares/Ribasim-NL/issues/149#issuecomment-2431933060 # Omkeren edgerichting -for edge_id in [131, 398, 406, 495, 513, 515, 894]: +for edge_id in [131, 398, 407, 495, 513, 515, 894]: model.reverse_edge(edge_id=edge_id) + # %% see: https://github.com/Deltares/Ribasim-NL/issues/149#issuecomment-2422164355 # Corrigeren netwerk bij Spuisluis Crèvecoeur @@ -234,6 +236,13 @@ model.tabulated_rating_curve.static.df = df +# %% see: https://github.com/Deltares/Ribasim-NL/issues/149#issuecomment-2431933060 +node_ids = [280, 335, 373, 879] +model.tabulated_rating_curve.static.df.loc[ + model.tabulated_rating_curve.static.df.node_id.isin([280, 335, 373, 879]), "active" +] = False + + # %% # level_boundaries updaten @@ -277,5 +286,6 @@ # %% write model model.use_validation = True model.write(ribasim_toml) +model.report_basin_area() # %% diff --git a/notebooks/aa_en_maas/01b_fix_basin_area.py b/notebooks/aa_en_maas/01b_fix_basin_area.py new file mode 100644 index 0000000..f27f28c --- /dev/null +++ b/notebooks/aa_en_maas/01b_fix_basin_area.py @@ -0,0 +1,243 @@ +# %% Import Libraries and Initialize Variables +import geopandas as gpd +import numpy as np +import pandas as pd +from ribasim.nodes import level_boundary, outlet + +from ribasim_nl import CloudStorage, Model, NetworkValidator + +# Initialize cloud storage and set authority/model parameters +cloud_storage = CloudStorage() +authority_name = "AaenMaas" +model_short_name = "aam" + +# Define the path to the Ribasim model configuration file +ribasim_model_dir = cloud_storage.joinpath(authority_name, "modellen", f"{authority_name}_fix_model_network") +ribasim_model_path = ribasim_model_dir / f"{model_short_name}.toml" +model = Model.read(ribasim_model_path) +network_validator = NetworkValidator(model) +# %% Load Input Geospatial Files +drainage_units_path = cloud_storage.joinpath( + authority_name, + "verwerkt", + "1_ontvangen_data", + "Na_levering_202404", + "afwateringseenheden_WAM", + "Afwateringseenheden.shp", +) +drainage_units_gdf = gpd.read_file(drainage_units_path, fid_as_index=True) + +# Load alternative drainage data +drainage_units_path = cloud_storage.joinpath( + authority_name, "aangeleverd", "Eerste_levering", "AfvoergebiedAanvoergebied.shp" +) +drainage_units_johnny_gdf = gpd.read_file(drainage_units_path, fid_as_index=True) + +# Load Ribasim model basin areas +ribasim_areas_path = cloud_storage.joinpath(authority_name, "verwerkt", "4_ribasim", "areas.gpkg") +ribasim_areas_gdf = gpd.read_file(ribasim_areas_path, fid_as_index=True, layer="areas") + +# Load node edit data +basin_node_edits_path = cloud_storage.joinpath(authority_name, "verwerkt", "model_edits.gpkg") +basin_node_edits_gdf = gpd.read_file(basin_node_edits_path, fid_as_index=True, layer="unassigned_basin_node") +basin_area_edits_gdf = gpd.read_file(basin_node_edits_path, fid_as_index=True, layer="unassigned_basin_area") +internal_basin_edits_gdf = gpd.read_file(basin_node_edits_path, fid_as_index=True, layer="internal_basins") + + +# replace unassigned basin area with baisn_area_edits +# 770 - 43 = 727 +model.basin.area.df = model.basin.area.df[~model.basin.area.df.index.isin(model.unassigned_basin_area.index)] +# 727 + 28 = 755 +df = basin_area_edits_gdf[basin_area_edits_gdf["to_node_id"].notna()] +df.loc[:, ["node_id"]] = df["to_node_id"].astype("int32") +model.basin.area.df = pd.concat([model.basin.area.df, df[["node_id", "geometry"]]]) + +# %% Assign Ribasim model ID's (dissolved areas) to the model basin areas (original areas with code) by overlapping the Ribasim area file baed on largest overlap +# then assign Ribasim node-ID's to areas with the same area code. Many nodata areas disappear by this method +combined_basin_areas_gdf = gpd.overlay(ribasim_areas_gdf, model.basin.area.df, how="union").explode() +combined_basin_areas_gdf["geometry"] = combined_basin_areas_gdf["geometry"].apply(lambda x: x if x.has_z else x) +combined_basin_areas_gdf["area"] = combined_basin_areas_gdf.geometry.area +non_null_basin_areas_gdf = combined_basin_areas_gdf[combined_basin_areas_gdf["node_id"].notna()] + +largest_area_node_ids = non_null_basin_areas_gdf.loc[ + non_null_basin_areas_gdf.groupby("code")["area"].idxmax(), ["code", "node_id"] +].reset_index(drop=True) + +combined_basin_areas_gdf = combined_basin_areas_gdf.merge(largest_area_node_ids, on="code", how="left") +combined_basin_areas_gdf["node_id"] = combined_basin_areas_gdf["node_id_y"] +combined_basin_areas_gdf.drop(columns=["node_id_x", "node_id_y"], inplace=True) +combined_basin_areas_gdf = combined_basin_areas_gdf.drop_duplicates(keep="first") +combined_basin_areas_gdf = combined_basin_areas_gdf.dissolve(by="code").reset_index() + +# %% The Ribasim model basins that have still nodata are being checked if they have overlap with aftwareringseenheden.shp. +# If overlap, they get ther Ribasim node-id where they have the largest overlap with +filtered_drainage_units_gdf = drainage_units_johnny_gdf[drainage_units_johnny_gdf["SOORTAFVOE"] != "Deelstroomgebied"] +filtered_drainage_units_gdf["geometry"] = filtered_drainage_units_gdf["geometry"].apply(lambda x: x if x.has_z else x) +filtered_drainage_units_gdf = filtered_drainage_units_gdf.to_crs(combined_basin_areas_gdf.crs) + +# Overlay filtered drainage units and updated basin areas +combined_basin_areas_johnny_gdf = gpd.overlay( + filtered_drainage_units_gdf, combined_basin_areas_gdf, how="union" +).explode() +combined_basin_areas_johnny_gdf = combined_basin_areas_johnny_gdf.dissolve(by="CODE") + +# Step 1: Separate unassigned from assigned units +unassigned_units_gdf = combined_basin_areas_gdf[combined_basin_areas_gdf["node_id"].isnull()] +unassigned_units_gdf["geometry"] = unassigned_units_gdf["geometry"].apply(lambda x: x if x.has_z else x) +assigned_units_gdf = combined_basin_areas_gdf[combined_basin_areas_gdf["node_id"].notna()] +assigned_units_gdf["geometry"] = assigned_units_gdf["geometry"].apply(lambda x: x if x.has_z else x) + +# Step 2: Calculate intersection areas between unassigned units and Johnny basins +overlap_gdf = gpd.overlay(combined_basin_areas_johnny_gdf, unassigned_units_gdf, how="union") + +# Step 3: Add overlap area for each polygon +overlap_gdf["overlap_area"] = overlap_gdf.geometry.area + +# Step 4: Select the largest overlap per code to assign `node_id` +largest_area_node_ids = overlap_gdf.loc[ + overlap_gdf.groupby("OBJECTID_2")["overlap_area"].idxmax(), ["node_id_1", "OBJECTID_2"] +].reset_index(drop=True) + +# Step 5: Merge largest node_id from overlaps back to unassigned units +unassigned_units_gdf = unassigned_units_gdf.merge( + largest_area_node_ids, left_on=["OBJECTID"], right_on=["OBJECTID_2"], how="outer" +) +unassigned_units_gdf["node_id"] = unassigned_units_gdf["node_id_1"] +unassigned_units_gdf.drop(columns=["node_id_1"], inplace=True) + +# Step 6: Merge unassigned node_ids back into the main dataset +basin_area_update = combined_basin_areas_gdf.merge( + unassigned_units_gdf[["OBJECTID", "node_id"]], + on="OBJECTID", + how="left", + suffixes=("", "_unassigned"), +) +# Step 7: Finalize missing `node_id` values from unassigned units +basin_area_update.loc[:, ["node_id"]] = basin_area_update["node_id"].fillna(basin_area_update["node_id_unassigned"]) +basin_area_update.drop(columns=["node_id_unassigned"], inplace=True) + +# %% If there are still nodata basins they are removed by assigning Nearest Basin ID +null_node_rows = basin_area_update[basin_area_update["node_id"].isnull()] + +if not null_node_rows.empty: + null_node_rows = null_node_rows.set_geometry(null_node_rows.geometry.centroid) + basin_area_update_centroid = basin_area_update.set_geometry(basin_area_update.geometry.centroid) + + nearest_basin = gpd.sjoin_nearest( + null_node_rows, + basin_area_update_centroid[basin_area_update_centroid["node_id"].notna()][["geometry", "node_id"]], + how="left", + distance_col="distance", + ) + basin_area_update.loc[basin_area_update["node_id"].isnull(), "node_id"] = nearest_basin["node_id_right"] +# basin_area_update.to_file("basin_area_update.gpkg", layer="basin_area_update") +# %% Based on basin_node_edits.gpkg, areas are assigned the Ribasim node_id that is in the file +basin_node_edits_notnull_gdf = basin_node_edits_gdf[basin_node_edits_gdf["to_area_code"].notna()] +merged_gdf = basin_area_update.merge( + basin_node_edits_notnull_gdf[["to_area_code", "node_id"]], + how="left", + left_on="code", + right_on="to_area_code", +) +merged_gdf["node_id"] = merged_gdf["node_id_y"].combine_first(merged_gdf["node_id_x"]) +merged_gdf.drop(columns=["node_id_x", "node_id_y"], inplace=True) + +# Dissolve geometries by `node_id` and save final GeoDataFrame +final_basins_gdf = merged_gdf.set_index("node_id").dissolve(by="node_id").reset_index() +# final_basins_gdf.to_file("basins_noholes.gpkg", layer="basins_noholes") + +final_basins_gdf.index.name = "fid" +model.basin.area.df = final_basins_gdf[["node_id", "geometry"]] + +# %% Check Differences in Node_ID Between Initial and Final Models +final_node_ids = final_basins_gdf["node_id"] +model_node_ids = model.basin.area.df["node_id"] +missing_in_model = final_basins_gdf[~final_basins_gdf["node_id"].isin(model_node_ids)] +missing_in_final = model.basin.area.df[~model.basin.area.df["node_id"].isin(final_node_ids)] +missing_gdf = pd.concat([missing_in_model, missing_in_final]) + +if "fid" in missing_gdf.columns: + missing_gdf = missing_gdf.rename(columns={"fid": "new_fid_name"}) + +# %% merge_basins +for row in basin_node_edits_gdf[basin_node_edits_gdf["to_node_id"].notna()].itertuples(): + if pd.isna(row.connected): + are_connected = True + else: + are_connected = row.connected + model.merge_basins(basin_id=row.node_id, to_basin_id=row.to_node_id, are_connected=are_connected) + +mask = internal_basin_edits_gdf["to_node_id"].notna() & internal_basin_edits_gdf["add_object"].isna() +for row in internal_basin_edits_gdf[mask].itertuples(): + if pd.isna(row.connected): + are_connected = True + else: + are_connected = row.connected + model.merge_basins(basin_id=row.node_id, to_basin_id=row.to_node_id, are_connected=are_connected) + +# %% add and connect nodes +for row in internal_basin_edits_gdf[internal_basin_edits_gdf.add_object.notna()].itertuples(): + from_basin_id = row.node_id + to_basin_id = row.to_node_id + if row.add_object == "stuw": + node_type = "TabulatedRatingCurve" + model.add_and_connect_node( + from_basin_id, int(to_basin_id), geometry=row.geometry, node_type=node_type, name=row.add_object_name + ) + +# %% reverse direction at node +for row in internal_basin_edits_gdf[internal_basin_edits_gdf["reverse_direction"]].itertuples(): + model.reverse_direction_at_node(node_id=row.node_id) + +# %% change node_type +for row in basin_node_edits_gdf[basin_node_edits_gdf["change_to_node_type"].notna()].itertuples(): + if row.change_to_node_type: + model.update_node(row.node_id, row.change_to_node_type, data=[level_boundary.Static(level=[0])]) + + +# %% corrigeren knoop-topologie +outlet_data = outlet.Static(flow_rate=[100]) +# 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]) + + +# %% +# basin-profielen/state 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 + + +model.write(ribasim_model_dir.with_stem("AaenMaas_fix_model_area") / "aam.toml") +model.report_basin_area() +model.report_internal_basins() +# %% diff --git a/notebooks/amstel_gooi_en_vecht/00_get_model.py b/notebooks/amstel_gooi_en_vecht/00_get_model.py new file mode 100644 index 0000000..36aa8fa --- /dev/null +++ b/notebooks/amstel_gooi_en_vecht/00_get_model.py @@ -0,0 +1,19 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +authority = "AmstelGooienVecht" +short_name = "agv" + +cloud = CloudStorage() + +model_url = cloud.joinurl(authority, "modellen", f"{authority}_parametrized_2024_11_20") +ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_parametrized_2024_11_20", "ribasim.toml") +if not ribasim_toml.exists(): + cloud.download_content(model_url) + +if ribasim_toml.exists(): # get a short_name version to differentiate QGIS layergroup + ribasim_toml.with_name(f"{short_name}.toml").write_text(ribasim_toml.read_text()) + +# %% diff --git a/notebooks/basin_area_nodes.py b/notebooks/basin_area_nodes.py index 403b0d0..e96b983 100644 --- a/notebooks/basin_area_nodes.py +++ b/notebooks/basin_area_nodes.py @@ -10,20 +10,45 @@ # %% 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]), - } - ] + ribasim_dir = cloud.joinpath(authority, "modellen", f"{authority}_fix_model_network") + if ribasim_dir.exists(): + ribasim_toml = next(ribasim_dir.glob("*.toml"), None) + if ribasim_toml is not None: + print(authority) + model = Model.read(ribasim_toml) + try: + valid = model.invalid_topology_at_node().empty + except KeyError: + valid = False + + unassigned_basin_area = model.unassigned_basin_area + if not unassigned_basin_area.empty: + unassigned_basin_area.to_file( + ribasim_dir.joinpath("basin_node_area_fouten.gpkg"), layer="area_niet_toegekend" + ) + area_mismatch = len(model.basin.node.df) - len(model.basin.area.df) + + unassigned_basin_node = model.basin.node.df[~model.basin.node.df.index.isin(model.basin.area.df.node_id)] + if not unassigned_basin_node.empty: + unassigned_basin_node.to_file( + ribasim_dir.joinpath("basin_node_area_fouten.gpkg"), layer="node_niet_toegekend" + ) + + data += [ + { + "waterschap": authority, + "model_valide": valid, + "basin_niet_toegekend": len(unassigned_basin_area), + "basin_knopen": len(model.basin.node.df), + "basin_vlakken": len(model.basin.area.df), + "basin_verschil": area_mismatch, + "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) +df.to_excel(cloud.joinpath("stand_modellen.xlsx"), index=False) + +# %% diff --git a/notebooks/brabantse_delta/01b_get_basin_holes.py b/notebooks/brabantse_delta/01b_get_basin_holes.py new file mode 100644 index 0000000..4dd4538 --- /dev/null +++ b/notebooks/brabantse_delta/01b_get_basin_holes.py @@ -0,0 +1,24 @@ +# %% +import geopandas as gpd +from shapely.geometry import MultiPolygon, Polygon + +from ribasim_nl import CloudStorage, Model + +cloud = CloudStorage() + +authority = "BrabantseDelta" +short_name = "wbd" + +ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_fix_model_network", f"{short_name}.toml") +model = Model.read(ribasim_toml) + +afwateringseenheden_df = gpd.read_file( + cloud.joinpath(authority, "verwerkt", "4_ribasim", "hydamo.gpkg"), layer="afwateringseenheden" +) +afwateringseenheden_poly = afwateringseenheden_df.buffer(0.01).buffer(-0.01).union_all() +waterschap_poly = MultiPolygon([Polygon(i.exterior) for i in afwateringseenheden_poly.geoms]) + +basin_polygon = model.basin.area.df.union_all() +holes_poly = waterschap_poly.difference(basin_polygon) +holes_df = gpd.GeoSeries(holes_poly.geoms, crs=model.basin.area.df.crs) +holes_df.to_file(cloud.joinpath(authority, "verwerkt", "basin_gaten.gpkg")) diff --git a/notebooks/brabantse_delta/01c_clip_basin_holes.py b/notebooks/brabantse_delta/01c_clip_basin_holes.py new file mode 100644 index 0000000..be6747b --- /dev/null +++ b/notebooks/brabantse_delta/01c_clip_basin_holes.py @@ -0,0 +1,46 @@ +# %% +import geopandas as gpd +import pandas as pd +from shapely.geometry import LineString, MultiPolygon, Polygon + +from ribasim_nl import CloudStorage +from ribasim_nl.geometry import split_basin, split_basin_multi_polygon + +authority = "BrabantseDelta" +cloud = CloudStorage() + +area_df = gpd.read_file(cloud.joinpath(authority, "verwerkt", "basin_gaten.gpkg"), fid_as_index=True) +cut_lines_df = gpd.read_file(cloud.joinpath(authority, "verwerkt", "knip_basin_gaten.gpkg"), fid_as_index=True) + + +new_area_df = gpd.GeoSeries([], crs=area_df.crs) + +for area_fid, cut_lines_select_df in cut_lines_df.groupby("poly_fid"): + # iterate trough cut_lines per area creating a series with geometries + series = gpd.GeoSeries([area_df.at[area_fid, "geometry"]], crs=area_df.crs) + + for row in cut_lines_select_df.itertuples(): + if not len(series[series.intersects(row.geometry)]) == 1: + raise ValueError(f"line with fid {row.Index} intersects basin with fid {area_fid} more than once") + + # split polygon into two polygons + fid = series[series.intersects(row.geometry)].index[0] + area_poly = series.loc[fid] + if isinstance(area_poly, MultiPolygon): + try: + result = split_basin_multi_polygon(series.loc[fid], line=row.geometry) + except ValueError: + result = split_basin_multi_polygon(series.loc[fid], line=LineString(row.geometry.boundary.geoms)) + elif isinstance(area_poly, Polygon): + result = split_basin(area_poly, line=row.geometry) + result = list(result.geoms) + + # original geometry gets first polygon + series.loc[fid] = result[0] + + # we add the second polygon to the series + series = pd.concat([series, gpd.GeoSeries([result[1]], crs=series.crs)], ignore_index=True) + + new_basin_area_df = pd.concat([new_area_df, series], ignore_index=True) + +new_basin_area_df.to_file(cloud.joinpath(authority, "verwerkt", "basin_gaten_geknipt.gpkg")) diff --git a/notebooks/brabantse_delta/01d_fix_basin_area.py b/notebooks/brabantse_delta/01d_fix_basin_area.py new file mode 100644 index 0000000..d56eef5 --- /dev/null +++ b/notebooks/brabantse_delta/01d_fix_basin_area.py @@ -0,0 +1,99 @@ +# %% Import Libraries and Initialize Variables + +import inspect + +import geopandas as gpd + +from ribasim_nl import CloudStorage, Model, NetworkValidator + +# Initialize cloud storage and set authority/model parameters +cloud_storage = CloudStorage() +authority_name = "BrabantseDelta" +model_short_name = "wbd" + +# Define the path to the Ribasim model configuration file +ribasim_model_dir = cloud_storage.joinpath(authority_name, "modellen", f"{authority_name}_fix_model_network") +ribasim_model_path = ribasim_model_dir / f"{model_short_name}.toml" + +ribasim_areas_path = cloud_storage.joinpath(authority_name, "verwerkt", "4_ribasim", "areas.gpkg") +ribasim_areas_gdf = gpd.read_file(ribasim_areas_path, fid_as_index=True, layer="areas") + +model = Model.read(ribasim_model_path) +network_validator = NetworkValidator(model) + +# Load node edit data +model_edits_url = cloud_storage.joinurl(authority_name, "verwerkt", "model_edits.gpkg") +model_edits_path = cloud_storage.joinpath(authority_name, "verwerkt", "model_edits.gpkg") +if not model_edits_path.exists(): + cloud_storage.download_file(model_edits_url) + +# %% +for action in gpd.list_layers(model_edits_path).name: + print(action) + # get method and args + method = getattr(model, action) + keywords = inspect.getfullargspec(method).args + df = gpd.read_file(model_edits_path, layer=action, fid_as_index=True) + for row in df.itertuples(): + # filter kwargs by keywords + kwargs = {k: v for k, v in row._asdict().items() if k in keywords} + method(**kwargs) + +# %% + +# remove unassigned basin area +model.remove_unassigned_basin_area() + +# %% 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") + +# 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") + +# %% Assign Ribasim model ID's (dissolved areas) to the model basin areas (original areas with code) by overlapping the Ribasim area file baed on largest overlap +# then assign Ribasim node-ID's to areas with the same area code. Many nodata areas disappear by this method +# Create the overlay of areas +combined_basin_areas_gdf = gpd.overlay(ribasim_areas_gdf, model.basin.area.df, how="union").explode() +combined_basin_areas_gdf["geometry"] = combined_basin_areas_gdf["geometry"].apply(lambda x: x if x.has_z else x) + +# Calculate area for each geometry +combined_basin_areas_gdf["area"] = combined_basin_areas_gdf.geometry.area + +# Separate rows with and without node_id +non_null_basin_areas_gdf = combined_basin_areas_gdf[combined_basin_areas_gdf["node_id"].notna()] + +# Find largest area node_ids for each code +largest_area_node_ids = non_null_basin_areas_gdf.loc[ + non_null_basin_areas_gdf.groupby("code")["area"].idxmax(), ["code", "node_id"] +] + +# Merge largest area node_ids back into the combined DataFrame +combined_basin_areas_gdf = combined_basin_areas_gdf.merge( + largest_area_node_ids, on="code", how="left", suffixes=("", "_largest") +) + +# Fill missing node_id with the largest_area node_id +combined_basin_areas_gdf["node_id"] = combined_basin_areas_gdf["node_id"].fillna( + combined_basin_areas_gdf["node_id_largest"] +) +combined_basin_areas_gdf.drop(columns=["node_id_largest"], inplace=True) +combined_basin_areas_gdf = combined_basin_areas_gdf.drop_duplicates() +combined_basin_areas_gdf = combined_basin_areas_gdf.dissolve(by="node_id").reset_index() +combined_basin_areas_gdf = combined_basin_areas_gdf[["node_id", "geometry"]] +combined_basin_areas_gdf.index.name = "fid" + +# %% +model.basin.area.df = combined_basin_areas_gdf + +# %% +model.use_validation = True +model.write(ribasim_model_dir.with_stem(f"{authority_name}_fix_model_area") / f"{model_short_name}.toml") +model.report_basin_area() +model.report_internal_basins() + +# %% diff --git a/notebooks/de_dommel/01_fix_model_network.py b/notebooks/de_dommel/01_fix_model_network.py index ad31018..436dca5 100644 --- a/notebooks/de_dommel/01_fix_model_network.py +++ b/notebooks/de_dommel/01_fix_model_network.py @@ -9,8 +9,10 @@ cloud = CloudStorage() -ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_2024_6_3", "model.toml") -database_gpkg = ribasim_toml.with_name("database.gpkg") +authority = "BrabantseDelta" +short_name = "wbd" + +ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_2024_6_3", f"{short_name}.toml") basin_data = [ @@ -260,3 +262,5 @@ ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_model_network", "model.toml") model.write(ribasim_toml) +model.report_basin_area() +model.report_internal_basins() diff --git a/notebooks/de_dommel/01b_fix_basin_area.py b/notebooks/de_dommel/01b_fix_basin_area.py new file mode 100644 index 0000000..d0122d9 --- /dev/null +++ b/notebooks/de_dommel/01b_fix_basin_area.py @@ -0,0 +1,53 @@ +# %% Import Libraries and Initialize Variables + +import inspect + +import geopandas as gpd + +from ribasim_nl import CloudStorage, Model, NetworkValidator + +# Initialize cloud storage and set authority/model parameters +cloud_storage = CloudStorage() +authority_name = "DeDommel" +model_short_name = "dommel" + +# Define the path to the Ribasim model configuration file +ribasim_model_dir = cloud_storage.joinpath(authority_name, "modellen", f"{authority_name}_fix_model_network") +ribasim_model_path = ribasim_model_dir / "model.toml" +model = Model.read(ribasim_model_path) +network_validator = NetworkValidator(model) + +# Load node edit data +model_edits_url = cloud_storage.joinurl(authority_name, "verwerkt", "model_edits.gpkg") +model_edits_path = cloud_storage.joinpath(authority_name, "verwerkt", "model_edits.gpkg") +if not model_edits_path.exists(): + cloud_storage.download_file(model_edits_url) + +# %% +for action in gpd.list_layers(model_edits_path).name: + # get method and args + method = getattr(model, action) + keywords = inspect.getfullargspec(method).args + df = gpd.read_file(model_edits_path, layer=action) + for row in df.itertuples(): + # filter kwargs by keywords + kwargs = {k: v for k, v in row._asdict().items() if k in keywords} + method(**kwargs) + +# %% 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") + +# 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") + +# %% +model.write(ribasim_model_dir.with_stem(f"{authority_name}_fix_model_area") / f"{model_short_name}.toml") +model.report_basin_area() +model.report_internal_basins() + +# %% diff --git a/notebooks/de_dommel/rare_edges.gpkg b/notebooks/de_dommel/rare_edges.gpkg deleted file mode 100644 index 141490d65b79acdb953d6866ff644372633e3683..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 106496 zcmeI53w#^Zndmio`aOBeiSsyCFp?TM@gt7ykPs&J#3QifM3xgw8iZnv?IC(N(iq#} zu{chkk1h8WE+x=i3VUG}?iSkZ50a2e0^#!9?UR;Ww%shHyDhXpFNM7!5GeP|XOp>?7%%@0{|!^mq)FI0vXZDcofhglR1^|0Xx>JFLnc@dUCEYrrJ zaQ8Nxp(7YVJ!~i27C~Vq*u(mmjosl$0C|;qp>zs~v7HGXc};cjdeJpPIuR53Eb33E zvSKDCq{OvIZKd2i)E`Tw(jwZ)qXeJLB5^2|LiJuh+QzZCroIrnZDpSPI#ln!2(h8| zy_eFisx25}^;8)5q|?c0g1?4OL`R0>xEe-Ce-fgp1Px*D=5~pHdI05Q&Jn0m$Fl*d zt3n~1b@@1T((6#23UghB0@E=Am}lvcIlZb+YcM}$f6P^^pf}C{-(*{v%UxSb-P9`g z1*4xC-WiO$e4kjoq@Q@OUpy6Fn zQB0!E;b0dN?n9TceR>7~T^I8QblC=Qwq!D@qa}@Q!H$rz1V0M1oA7WFYGZp8+;}(= z`IMN|q4xVXlsnx$wUk9j#rd)9?gXAmqp=Y&Ex(T%%@nOS##Br$EpfSderMS{*-9eJjt=^zq7X|&GyFg_JDx>8RnXs2>U3{!nz4m3d9#%zdz9_7 zxLoe$X6oi;xehYTY(zj?!5gcUDQ6)XJ>h`}EZg;Q3x+AXg->&yL{S{dzCy4F$LK;)y=j zi}U1hAD84)`G(gC&!}q~y{SKX?!FW(cDi>oPwPvXQB;qJHWc)Ul3>x97p06oMh(Ri z3mh)@s#Vm)_MBs;2Q%1O%77^m)7`sjT3mX)YQNJ*xh5A|U2dBbJVP60^+~wnq%(0mRTZ4R zDYbA=5Y6=Z9UsP1q3{S4z{}^ zKhYDw`?cOKy}T#Tvl$;ubO+nFdp6?p3O<4ATUOV50&5zZn;To3SJpQ*1sYmg@pD5< zeM3ut>5XvR;b5e1yGN>nC*WzWt6S~av7RVe^H#Qc?+FIYxP`D1OwK{ZPiqCw4 z-J$KC?wtu?_Xr>D5%~;{cjU3*A)()c>&b+fE;ho3w|np*k06f3`Dl6|n#QMELdp|p zYQPlM)UT;)k{@gGGW6qpbOv9d!=@7ALQE20+1S#&rm;nRwC0It`(yZcY#`PzrZZ9- zTk)f=IbUWVmdNr^9&fjW)SxGTEm2ojmsg4+k;0BKB#W$G(gH!$?YGLhwCECd}`UQoR+p&2A{=Y z+C<66TH)FYloOgp$>h9y2Z~qV?zT^uShB672equ0su$-)pf00e*l5C8%|00;m9An<7rh*720ORCG) zyYR!BiQ)TY`1l~7Rd1fj=3*M^8dlfVwbs@(Md};Z)-|kcY+Z??u3#)DS}wfvX_{Hk zQ@QZHy6-afJf00cfa1j=pI#kng7u>b$u7#11=1b_e#00KY&2mk>f00e*l5C8%|Kt~{d|F7Kt z|F(sGTPFf9fdCKy0zd!=00AHX1b_e#00KY&2mpcenm`3*tEO_t|9&|CKd;RRl?MVq z00;m9AOHk_01yBIKmZ5;0U$61fkNm1CoJ@dDa0Wf2mk>f00e*l5C8%|00;m9AOHk_ z01)_W6R4(K)$^;ZMMb#-f9dx>JA)w>|B`Yrgn!{gIsbplLcjIdZVJ>62mk>f00e*l z5C8%|00;m9AOHk_01%iVfjRc-`MU4^ldc88{r@v$2Vy_~2mk>f00e*l5C8%|00;m9 zAOHl;D+1>G|1Go5s|G-2fdCKy0zd!=00AHX1b_e#00KY&2mpaIBp{vt!~Xvap@BFM z00KY&2mk>f00e*l5C8%|00;nq^NfJ>`+upP-|I4I{5Hy-*omw2|G^9cEGQ=5{G_ z_a>w&lr^F5khvT$!Z{Gjv~eihy$u)85e%Uowv%m(pfD5cVSUWT?ri9@jzs`vWQHjc$L z_l4MPEAuM74%PcFLTspgZ~H6->+Bg}6;fF~BO0-t8CInX4Qnve!-ga0h*2r?Q^a^F z?W)>>y|a2MjC<1QWHiBF!za$xH%%pI2vyHa51^dPIRZUXX%kSz0uE>Kjtci&Jy2-a;Lkema+(`I6s!%oe)Gm8XFPQ^82Wm*PsNwWcJYZ%w{M^=%Gh zGufyRHybA@q@rUcxKix6IW{&{NEXMMW%E%EiplIJ#?oLWxoq4#X_C(wNNKOxIN9Nt zD8v%c3_lRfj%Sfi74&zMI^ElwX53FPd2?Tt>`}JQ;&Qp0o2i?Z<@&4i{mZeyN#;`lIMNMI`8*n=X$P@mL^v|g7qN9in$`%Z=dN{dzCy4F$LK;`LIl6Xwa|E-J~V@(r&Oo>9*~!yFp4OK%qo^JcZ7ApyCBdRIFG`t1j2en37C2n)Rja6p?K#Ix4`#5nl<`Wg zfIh-$c1V4)&MwmoJOg}89Lexl%Sb9WIkdBznBGoh=J7L5r@MF6w7B$o)qbasa!oF_ zy4+qbwSP3%w!C(cH#T#dNRr;srfCq~$;6^kSIU)8BEhcQ42OuI4{tICL?OxhyfuBbHObnVc%+6~ zThq0+X5b2Mz?({sdj0-@ek2bd{RkdN>IQC{q#Lje^MbiEWTfsO%^YVlx@v=KtKI4D z@=j}EVgQS#h7>CZWbVIGiu-wqbPxMR0 zv)%3$D8snn!k)aGF4p#b!37H2 zsu{68%iW(|_h)SXr>`11N0xo_XXK~Ke&^i$G_Q;){Z!pqeOi~y6t-$5&Wx>^N$02Y zX5^#AK<9o0qq z$<$SNt)X05tW-TD9ib)N7{CIq`BXdE7K!g`46G)&i0j_ z0U_=|>ydN@puA>6?aCCaX2nyobUhA9@zV6%6Nv`1Y;sr}mnUjP%I{Z?G03jxNsR;HJ63mNKuNrSd zt}3>=XkU-&=>;pK`eu&Q%;+oKg1{#T3L*EU_|cU)u=178^PKEN=~9=QVW^4fX;$LW z%D6CbDW#HNmah>P_LZn4IPDr-+6en6=p`<9TN^b|ldBCJ4#MH!)MuewG z!4Qj2C4!-j0P+N-;~)>&+u{i6DuPGe!{avt8GIF0>VSCDuNM;2g<_@J%U7hd{MnJ6 zDP9~+XLiX2NY|bG4!p6a=+n~neja_gVr%&uhmMoh)`3KNbX|zKSh&Ef>MO6j_|EJxig|U(>7S!X>s*qXL!5#Di7cH_f6Gwt&r;wFtedPI*wd*RG*oXhT zXg_||T15m&f-!^rXO7d)p(?v+qCp zXi(QGB1i;C1|n$xf&CR~Qd70mtm^$|lB?fEcIKvMFFRD#WM_{4+q@k=9(6NEntpJ) z^~*X{MgfwM2-=@@lu)xZRZGmOd`H_RznHKyqsP9r=$W*gdH$PE-(%n4X8!mg6?nEo zr^+ZmG7>?%=vZcXPE)nmtcuEf<(9-+8}nk|hYt;}w=sJUeC^4le_X&EdueN;o6@N= z3XqIM&>nYAP&aC-s?4h1UG`%7!UC#}hQD!d?{un2Mj~k6iB;X8saj-Kb@Pt~UMeo2 z>h!G-TBmIl$w&n4F|6u(P1QoPs*ykY=q`5wRn3bwEelWADw2^1+OKj>SPpBdD$S~{ z-_UgBsn;#cA9frGKJZ%$v$5#gi@$cin>qgX_dfH&gF0Io1xQ9BXur}v*K$x(afAef_PE{{eU9fBxZn9Um@WEPd|LKey^yWfUM8iJAh^nFa zswyfsUzOeV;NIs8sA@m{zQt6%nTRc2NH?aROM{$M&) zBqI^Dy=`CKM>Gku&5e3-Ys$CrE*nFbTd=sWPU1)9Y)RYZ^okPO6fwi@exS3*-on^k?Z^Y>TGd6#0y zv}av$^VU086qpx@AQ2!L$UJHN$hy{XT2obFRz((ABqv#DkrHyPQkQN-fMg(o)(@?% zI#uOnRUeFtd*9n@W4`{@bx#cMw=wNsO}tkAr3K7C+~|AzR{X4~G76B4M9}J?*16JJ ztIEu(c4J#D?Xok*CCPp~I=y^_&ex3sBqI^D&bEHjF;`PnYF0(o2mj$e_WL*PvuO*8 z*M6{mmXvIY(|~R%BKs%|Z5N8)IB2Jof$#zkANE zYn4%eWYh%bQa9RjYZ%wms=n+`{r3o-iXE$GuaRbXtpu_SF$$Pe&7uBGr^-3CD&vyG zq{=8@QZ<{pL8rAn~f>Gz<0wj{Ji0+UvGBcg~jvh$18W@XUzkNAQ2!L zh@jSspj70D{>qGW8HBFRR z)4N|{uIgGuF=R@;?vB*0$MA$mrdF*{M1Ulaj3mR=mP0Nt)`Yjrb8DpU|D)))E%*-} zKmZ5;0U!VbfB+Bx0zd!=00AHX1c1Qjgg^yltEO`IH29_ce+@pgr@ui5DvrbXKb-&P zZV-U;f8`Doa;#6T!jt2Day4GN@_+A7!|Ces} zBggp0I|<0KzBbqz1xQ9B2hI!@7BGyu(D?Xv;V-p@88*{`v4CSAPFQR$>6%r+Hb#` z+$b;j%&LAlt7Z7(H>}Lxe*f6~qif00e*l5C8%| z00;nqPnf_}l(Tw%?e6N_4g1m;$;v-VpnQ*~{CiPDqBy*E)vD3a(UmxPIM%-_Hps8U zKNqknJIwdjH`J|K&X3_AAIjq26o|_CMl<*ar)txw#CQp=F-8B21^>YV2mk>f00e*l z5C8%|00;m9AOHk_01!A22`qQuiwDKiUptV_|NV5Mg+5K6r2jy_Nxw$_iatTVO#hVr zA^igVH2nm9gnpF%F8#OkU(@%}e?i|t-%5Xt{t`V&-$0MkBAuqMrg{1*dON+9-a>cM z9rQ+eJ>7zfga;4+0zd!=00AHX1b_e#00KY&2mpcenSjl2vlQEWQs9+>3#Fh&3NDa> z3KmL1r4+cOV1X3Omx6gxFjorZNWp9=m?Z_Y z6jVq-xfGO1L8%m!NI|g_6iIz>&Ud(_P8)Z4 z-cIiDJxkSi8@Ko`Z^8fn!5u&4=Z>GAt;X9q*NP3CtK)xhM?)d*XzEop-oY(A?BN#v zxSKm%GR7UgkWu5EO8ifnaQs8Mc$X6Y;$Iu&iG`XLvcHPB29{eHqc$8A(OSvVM+qos)7 z*9dob(eJp!n|7)3bFr=dPcB}w*F0SI-Q4mYe}v0^ggf}mQB{6EwqZNB^eT!w^s9@x z!}DEg`~nKeX)p$3qONLt#d=%Fu!98W$sK!@f+djp)y4tv-fkABJ&1(EYY}?1V zg>T%<9eNx4P|4TT_(eFrOBcUKjj!Ss{-jYC|F9arm|OVMxGw&BHGT=U?Z0!&Z}4-E zovOh$Ca1DXv271<)ernD_dkzKU>pBfm0!j!KJ@q4#&>X!y*U@#n3RvOO-^vLT~Bex zUrp)kQjKjsjD6??_MJJ{_G{I0Je=pxpXNM&bs=|jD8n7yRi(x+=bXPi&N)x7;*PDX z;f{5_uf|`%l|Fe3SNejBd$zuZdv@(=HNJ+U%MN2-Im|uz9)o>_^rH*8#ndmbZyn_x zX?sQIS6-!my*Nklhf^!ma(tXS`T^(Wf6IOEg{9p0zW*&X-p{#x*K_WxOS$hIcml^0 z4r=n}|4G{aTV{R!+X*cM0zd!=00AHX1b_e#00KY&2mk>f@c%ObY5z}Khb;8N6^~Sm zl)qlSwe0z_l_ftb?k#F|-t36mZ?cWzWOx99b0V;}%LVk+I z`2;WW(E&atj%4^~e|jV(MsuH8$-Qa|vrL3Vk#Ml1gAF6k3=8leW>W;qx3irr&fXnD zo?Mw8w2|!yhEVV3cByRni#5`N8C;7L{=N(`i+b1y%Bvb$Z+Kse>d`ii4YMeb9_2GW zF}5?oM^mvR@5i!AoHwV@>%R!Gq4vFH7Z%1Lq_UXJIpFvc1rQv;P!Ai9eDXL}oE46V z?KQ3{wz_EVcGarNx0~d!MWydOQGF#BYjla}a*>8l;Ro1+A(T$US9&B(9u(@1pipmT zr|KYTu1(?YF6=*IhUaxYDp$^N9PgnvJwib!DTU%k9&ug{s%15W7vOBVt1B$S1#5+8*qYTizW;J&|ZI%O;1#ab3FK ze>PBXFY{^o13B$Vf0(RYE+g81{R|K_o*kmbPY@Ay51EmKPI8q`5LBM`kZu}AVRmyT z)5fA;DAGNp3vOjPd)XevCwF#-6^G`VUf4XnynM@kpud7zPxoQS<>d%f&{r?L(XXe(QS?An_x`NU@)_Za?8(=5n(vHBm1wa|%v`(cxGI zubrk#TSz8uW@9Hi#bBluviJ zVQfqah19tGdfdOh%IV&+g0jfpUYFgSz+Xlejg5$D`F&Jgg==3}myOmNA{I=HE^@h7 zte_?px`t%K+xdAE4T(s%Tt^Z1@<^aRoft``g-)sAeX pd.DataFrame: df = pd.read_feather(filepath) @@ -147,7 +177,7 @@ def remove_node(self, node_id: int, remove_edges: bool = False): for attr in table.model_fields.keys(): df = getattr(table, attr).df if df is not None: - if node_id in df.columns: + if "node_id" in df.columns: getattr(table, attr).df = df[df.node_id != node_id] else: getattr(table, attr).df = df[df.index != node_id] @@ -166,7 +196,7 @@ def remove_node(self, node_id: int, remove_edges: bool = False): if node_id in table._parent._used_node_ids: table._parent._used_node_ids.node_ids.remove(node_id) - def update_node(self, node_id, node_type, data, node_properties: dict = {}): + def update_node(self, node_id, node_type, data: list | None = None, node_properties: dict = {}): existing_node_type = self.node_table().df.at[node_id, "node_type"] # read existing table @@ -192,6 +222,8 @@ def update_node(self, node_id, node_type, data, node_properties: dict = {}): # add to table table = getattr(self, pascal_to_snake_case(node_type)) + if data is None: + data = getattr(DEFAULT_TABLES, pascal_to_snake_case(node_type)) table.add(Node(**node_dict), data) # sanitize node_dict @@ -296,6 +328,122 @@ 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 add_basin(self, node_id, geometry, tables=None, **kwargs): + # define node properties + if "name" in kwargs.keys(): + name = kwargs["name"] + kwargs.pop("name") + else: + name = "" + + node_properties = {k if k.startswith("meta_") else f"meta_{k}": v for k, v in kwargs.items()} + + # define tables, defaults if None + if tables is None: + tables = DEFAULT_TABLES.basin + + self.basin.add(Node(node_id=node_id, geometry=geometry, name=name, **node_properties), tables=tables) + + def connect_basins(self, from_basin_id, to_basin_id, node_type, geometry, tables=None, **kwargs): + self.add_and_connect_node( + from_basin_id=from_basin_id, + to_basin_id=to_basin_id, + node_type=node_type, + geometry=geometry, + tables=tables, + **kwargs, + ) + + def add_and_connect_node(self, from_basin_id, to_basin_id, geometry, node_type, tables=None, **kwargs): + # define node properties + if "name" in kwargs.keys(): + name = kwargs["name"] + kwargs.pop("name") + else: + name = "" + node_properties = {k if k.startswith("meta_") else f"meta_{k}": v for k, v in kwargs.items()} + + # define tables, defaults if None + if tables is None: + tables = getattr(DEFAULT_TABLES, pascal_to_snake_case(node_type)) + + # add node + node = getattr(self, pascal_to_snake_case(node_type)).add( + Node(geometry=geometry, name=name, **node_properties), tables=tables + ) + + # add edges from and to node + self.edge.add(self.get_node(from_basin_id), node) + self.edge.add(node, self.get_node(to_basin_id)) + + def add_basin_outlet(self, basin_id, geometry, node_type="Outlet", tables=None, **kwargs): + # define node properties + if "name" in kwargs.keys(): + name = kwargs["name"] + kwargs.pop("name") + else: + name = "" + + node_properties = {k if k.startswith("meta_") else f"meta_{k}": v for k, v in kwargs.items()} + + # define tables, defaults if None + if tables is None: + tables = getattr(DEFAULT_TABLES, pascal_to_snake_case(node_type)) + + # add outlet + node = getattr(self, pascal_to_snake_case(node_type)).add( + Node(geometry=geometry, name=name, **node_properties), tables=tables + ) + + # add edges from and to node + self.edge.add(self.basin[basin_id], node) + edge_geometry = self.edge.df.set_index(["from_node_id", "to_node_id"]).at[(basin_id, node.node_id), "geometry"] + + # add boundary + geometry = shapely.affinity.scale(edge_geometry, xfact=1.05, yfact=1.05, origin="center").boundary.geoms[1] + # geometry = edge_geometry.interpolate(1.05, normalized=True) + boundary_node = self.level_boundary.add(Node(geometry=geometry), tables=DEFAULT_TABLES.level_boundary) + self.edge.add(node, boundary_node) + + def reverse_direction_at_node(self, node_id): + for edge_id in self.edge.df[ + (self.edge.df.from_node_id == node_id) | (self.edge.df.to_node_id == node_id) + ].index: + self.reverse_edge(edge_id=edge_id) + + def update_basin_area(self, node_id: int, basin_area_fid: int | None = None): + if basin_area_fid is None: + basin_area_fid = self.basin.area.df.index.max() + 1 + self.basin.area.df.loc[basin_area_fid, ["node_id"]] = node_id + + def add_basin_area(self, geometry: MultiPolygon, node_id: int | None = None): + # if node_id is None, get an available node_id + if pd.isna(node_id): + basin_df = self.basin.node.df[self.basin.node.df.within(geometry)] + if basin_df.empty: + raise ValueError("No basin-node within basin area, specify node_id explicitly") + elif len(basin_df) > 1: + raise ValueError( + f"Multiple basin-nodes within area: {basin_df.index.to_numpy()}. Specify node_id explicitly" + ) + else: + node_id = basin_df.index[0] + elif node_id not in self.basin.node.df.index: + raise ValueError(f"Node_id {node_id} is not a basin") + + # check geometry and promote to mulitpolygon + if not geometry.geom_type == "MultiPolygon": + if geometry.geom_type == "Polygon": + geometry = MultiPolygon([geometry]) + else: + raise ValueError(f"geometry-type {geometry.geom_type} is not valid. Provide (Multi)Polygon instead") + + # if all correct, assign + area_df = gpd.GeoDataFrame({"node_id": [node_id], "geometry": [geometry]}, crs=self.crs) + area_df.index.name = "fid" + area_df.index += self.basin.area.df.index.max() + 1 + self.basin.area.df = pd.concat([self.basin.area.df, area_df]) + def move_node(self, node_id: int, geometry: Point): node_type = self.node_table().df.at[node_id, "node_type"] @@ -311,6 +459,18 @@ def move_node(self, node_id: int, geometry: Point): ].index.to_list() self.reset_edge_geometry(edge_ids=edge_ids) + def report_basin_area(self): + gpkg = self.filepath.with_name("basin_node_area_errors.gpkg") + self.unassigned_basin_area.to_file(gpkg, layer="unassigned_basin_area") + + unassigned_basin_node = self.basin.node.df[~self.basin.node.df.index.isin(self.basin.area.df.node_id)] + unassigned_basin_node.to_file(gpkg, layer="unassigned_basin_node") + + def report_internal_basins(self): + gpkg = self.filepath.with_name("internal_basins.gpkg") + df = self.basin.node.df[~self.basin.node.df.index.isin(self.edge.df.from_node_id)] + df.to_file(gpkg) + def find_closest_basin(self, geometry: BaseGeometry, max_distance: float | None) -> NodeData: """Find the closest basin_node.""" # only works when basin area are defined @@ -459,54 +619,98 @@ def redirect_edge(self, edge_id: int, from_node_id: int | None = None, to_node_i 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") + def deactivate_node(self, node_id: int): + node_type = self.get_node_type(node_id) + df = getattr(self, pascal_to_snake_case(node_type)).static.df + df.loc[df.node_id == node_id, ["active"]] = False - if are_connected: - paths = [i for i in nx.all_shortest_paths(self.graph, basin_id, to_basin_id) if len(i) == 3] + def remove_unassigned_basin_area(self): + df = self.basin.area.df[~self.basin.area.df.index.isin(self.unassigned_basin_area.index)] + if self.basin.area.df.node_id.duplicated().any(): + df = df.dissolve(by="node_id").reset_index() + df.index.name = "fid" + self.basin.area.df = df + + def merge_basins( + self, + basin_id: int | None = None, + node_id: int | None = None, + to_node_id: int | None = None, + to_basin_id: int | None = None, + are_connected=True, + ): + if basin_id is not None: + warnings.warn("basin_id is deprecated, use node_id instead", DeprecationWarning) + node_id = basin_id + + if to_basin_id is not None: + warnings.warn("to_basin_id is deprecated, use to_node_id instead", DeprecationWarning) + to_node_id = to_basin_id + + if node_id not in self.basin.node.df.index: + raise ValueError(f"{node_id} is not a basin") + to_node_type = self.node_table().df.at[to_node_id, "node_type"] + if to_node_type not in ["Basin", "FlowBoundary", "LevelBoundary"]: + raise ValueError( + f'{to_node_id} not of valid type: {to_node_type} not in ["Basin", "FlowBoundary", "LevelBoundary"]' + ) + + if are_connected and (to_node_type != "FlowBoundary"): + paths = [i for i in nx.all_shortest_paths(self.graph, node_id, to_node_id) if len(i) == 3] if len(paths) == 0: - raise ValueError(f"basin {basin_id} not a direct neighbor of basin {to_basin_id}") + raise ValueError(f"basin {node_id} not a direct neighbor of basin {to_node_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() + edge_ids = self.edge.df[self.edge.df.from_node_id == node_id].index.to_list() + edge_ids += self.edge.df[self.edge.df.to_node_id == node_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 + self.edge.df.loc[self.edge.df.from_node_id == node_id, "from_node_id"] = to_node_id + self.edge.df.loc[self.edge.df.to_node_id == node_id, "to_node_id"] = to_node_id + + # remove self-connecting edge in case we merge to flow-boundary + if to_node_type == "FlowBoundary": + mask = (self.edge.df.from_node_id == to_node_id) & (self.edge.df.to_node_id == to_node_id) + self.edge.df = self.edge.df[~mask] # 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_node_type == "Basin": + if node_id in self.basin.area.df.node_id.to_numpy(): + poly = self.basin.area.df.set_index("node_id").at[node_id, "geometry"] + + # if to_node_id has area we union both areas + if to_node_id in self.basin.area.df.node_id.to_numpy(): + poly = poly.union(self.basin.area.df.set_index("node_id").at[to_node_id, "geometry"]) + if isinstance(poly, Polygon): + poly = MultiPolygon([poly]) + self.basin.area.df.loc[self.basin.area.df.node_id == to_node_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_node_id, + "geometry": poly, + } - # 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 + if self.basin.area.df.crs is None: + self.basin.area.df.crs = self.crs - # 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} + # if node type is flow_boundary, we change type to LevelBoundary + if to_node_type == "FlowBoundary": + self.update_node(to_node_id, "LevelBoundary", data=[level_boundary.Static(level=[0.0])]) # 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 + self.remove_node(node_id) def invalid_topology_at_node(self, edge_type: str = "flow") -> gpd.GeoDataFrame: df_graph = self.edge.df diff --git a/src/ribasim_nl/ribasim_nl/network_validator.py b/src/ribasim_nl/ribasim_nl/network_validator.py index 057fee3..6f03683 100644 --- a/src/ribasim_nl/ribasim_nl/network_validator.py +++ b/src/ribasim_nl/ribasim_nl/network_validator.py @@ -12,14 +12,20 @@ def within_distance(row, gdf, tolerance=1.0) -> bool: def check_node_connectivity(row, node_df, tolerance=1.0) -> bool: invalid = True + # handle zero-length edges + if row.geometry.length == 0: + point_from = point_to = row.geometry.centroid + else: + point_from, point_to = row.geometry.boundary.geoms + # 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"]) + distance = point_from.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"]) + distance = point_to.distance(node_df.at[row.to_node_id, "geometry"]) invalid = distance > tolerance return invalid