De Dommel met bergend gebied (#158)

Een model van De Dommel met bergingsknopen volgens methode Vd Gaast.
Resulteert in model DeDommel versie 2024.9.0, beschikbaar via:
https://deltares.thegood.cloud/f/119263

fixes/updates:
- #102
- #157
- #97

notebooks/de_dommel/01_fix_model_network.py

@@ -4,7 +4,7 @@
 from ribasim import Node
 from ribasim.nodes import basin, level_boundary, manning_resistance, outlet
 from ribasim_nl import CloudStorage, Model, NetworkValidator
-from shapely.geometry import Point
+from shapely.geometry import Point, Polygon
 
 cloud = CloudStorage()
 
@@ -117,13 +117,13 @@
 model.level_boundary.add(boundary_node, [level_data])
 
 model.edge.add(model.tabulated_rating_curve[614], model.level_boundary[28])
-
+# %%
 # see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292014475
 model.remove_node(node_id=1898, remove_edges=True)
 
 # see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292017813
-for node_id in [1891, 989, 1058]:
-    model.remove_node(node_id, remove_edges=True)
+model.update_node(989, "Outlet", [outlet.Static(flow_rate=[0])])
+model.update_node(1891, "LevelBoundary", [level_data])
 
 # see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291988317
 # for from_node_id, to_node_id in [799, 1580, 625, 1123, 597, 978]:
@@ -156,11 +156,49 @@
 if not network_validator.node_internal_basin().empty:
     raise Exception("nog steeds interne basins")
 
+
+# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2367724440
+gdf = gpd.read_file(
+    cloud.joinpath("DeDommel", "verwerkt", "4_ribasim", "hydamo.gpkg"),
+    layer="stuw",
+    fid_as_index=True,
+)
+kst = gdf.loc[35]
+geometry = Point(kst.geometry.x, kst.geometry.y)
+name = kst.CODE
+meta_object_type = "stuw"
+
+outlet_node_id = model.next_node_id
+
+kst_node = model.outlet.add(
+    Node(node_id=outlet_node_id, geometry=geometry, name=name, meta_object_type=meta_object_type),
+    [outlet_data],
+)
+
+
+gdf = gpd.read_file(
+    cloud.joinpath("DeDommel", "verwerkt", "4_ribasim", "hydamo.gpkg"),
+    layer="hydroobject",
+    engine="pyogrio",
+    fid_as_index=True,
+)
+geometry = gdf.at[2822, "geometry"].interpolate(0.5, normalized=True)
+basin_node_id = model.next_node_id
+basin_node = model.basin.add(
+    Node(node_id=basin_node_id, geometry=geometry, meta_krw_name="Witte Loop/Peelrijt", meta_krw_id="NL27_KD_3_2"),
+    basin_data,
+)
+
+model.remove_edge(from_node_id=664, to_node_id=8, remove_disconnected_nodes=False)
+model.edge.add(model.manning_resistance[664], basin_node)
+model.edge.add(basin_node, kst_node)
+model.edge.add(kst_node, model.level_boundary[8])
+
+# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2293486609
 df = network_validator.edge_incorrect_type_connectivity(
     from_node_type="ManningResistance", to_node_type="LevelBoundary"
 )
 
-# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2293486609
 for node_id in df.from_node_id:
     model.update_node(node_id, "Outlet", [outlet.Static(flow_rate=[100])])
 
@@ -174,9 +212,50 @@
 
 model.basin.area.df = model.basin.area.df[~model.basin.area.df.node_id.isin(model.unassigned_basin_area.node_id)]
 
-# # %% write model
+
+# fix basin area
+
+# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2370880128
+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]
+geoseries = gpd.GeoSeries(holes, crs=28992)
+
+drainage_areas_df = gpd.read_file(
+    cloud.joinpath("DeDommel", "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(geoseries):
+    # 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.buffer(-10).within(drainage_area)].set_index(
+            "node_id"
+        )
+        intersecting_basins_df = selected_basins_df.intersection(geometry.buffer(10))
+        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)
+        )
+        model.basin.area.df.loc[model.basin.area.df.node_id == assigned_basin_id, "geometry"] = geometry
+
+#  %% write model
 model.edge.df.reset_index(drop=True, inplace=True)
 model.edge.df.index.name = "edge_id"
-ribasim_toml = ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_model_network", "model.toml")
+ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_model_network", "model.toml")
 
 model.write(ribasim_toml)

notebooks/de_dommel/03_fix_basin_area.py

@@ -23,13 +23,14 @@
 dissolved_area_gdf.to_file(cloud.joinpath("DeDommel", "verwerkt", "water_area.gpkg"))
 
 # %%
-basin_df = model.basin.node.df
-basin_area_df = gpd.read_file(
-    cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg"), engine="pyogrio", fid_as_index=True
-)
+basin_area_gpkg = cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg")
+basin_area_df = model.basin.area.df
+basin_area_df.to_file(basin_area_gpkg)
 basin_area_df.set_index("node_id", inplace=True)
 
+basin_df = model.basin.node.df
 
+# %%
 data = []
 ignore_basins = [1278, 1228, 1877, 1030]
 row = next(i for i in basin_df.itertuples() if i.Index == 1230)
@@ -76,6 +77,8 @@
 area_df = gpd.GeoDataFrame(data, crs=model.basin.node.df.crs)
 area_df = area_df[~area_df.is_empty]
 area_df.index.name = "fid"
+mask = area_df.geometry.type == "Polygon"
+area_df.loc[mask, "geometry"] = area_df.geometry[mask].apply(lambda x: MultiPolygon([x]))
 model.basin.area.df = area_df
 # %%
 ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_areas", "model.toml")

notebooks/de_dommel/04_parameterize_model.py

@@ -357,6 +357,8 @@ def get_area_and_profile(node_id):
         flow_rate = round(
             basin_area_df.at[model.upstream_node_id(node_id), "geometry"].area * 0.015 / 86400, 2
         )  # 15mm/day of upstream areay
+    else:
+        flow_rate = flow_rate / 60  # m3/min -> m3/s
 
     data = pump.Static(flow_rate=[flow_rate], min_upstream_level=min_upstream_level)