Fix dommel edges (#137)

commits for #136:
- `reset_edge_geometry()` for ribasim_nl.Model class to reset
edge_geometries to a straight line
- re-organized processing scripts:
   - `01_fix_model_network.py` contains all for #102 
   - `02_fix_edges.py`: fixes edges to follow HydroObject
- `03_fix_basin_area.py`: assigns water_areas as basin / area replacing
entire catchment area
  
Results in model-version 2024.8.3, available at:
https://deltares.thegood.cloud/f/118078

notebooks/de_dommel/02_fix_edges.py

@@ -0,0 +1,61 @@
+# %%
+import geopandas as gpd
+from ribasim_nl import CloudStorage, Model, Network
+
+cloud = CloudStorage()
+
+
+# %% load model
+ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel", "model.toml")
+model = Model.read(ribasim_toml)
+
+# %% network from HydroObjects
+network_gpkg = cloud.joinpath("DeDommel", "verwerkt", "network.gpkg")
+if network_gpkg.exists():
+    network = Network.from_network_gpkg(network_gpkg)
+else:
+    network = Network.from_lines_gpkg(
+        cloud.joinpath("DeDommel", "verwerkt", "2_voorbewerking", "hydamo.gpkg"), layer="hydroobject"
+    )
+
+    network.to_file(network_gpkg)
+# %% edges follow HydroObjects
+model.reset_edge_geometry()
+node_df = model.node_table().df.set_index("node_id")
+data = []
+for row in model.edge.df.itertuples():
+    try:
+        # get or add node_from
+        from_point = node_df.at[row.from_node_id, "geometry"]
+        distance = network.nodes.distance(from_point)
+        if distance.min() < 0.1:
+            node_from = distance.idxmin()
+        else:
+            node_from = network.add_node(from_point, max_distance=10, align_distance=1)
+
+        # get or add node_to
+        to_point = node_df.at[row.to_node_id, "geometry"]
+        distance = network.nodes.distance(to_point)
+        if distance.min() < 0.1:
+            node_to = distance.idxmin()
+        else:
+            node_to = network.add_node(to_point, max_distance=10, align_distance=1)
+
+        if (node_from is not None) and (node_to is not None):
+            # get line geometry
+            geometry = network.get_line(node_from, node_to)
+
+            # replace edge geometry
+            model.edge.df.loc[row.Index, ["geometry"]] = geometry
+        else:
+            print(f"edge not updated for {row.Index} as node_from and node_to cannot be found")
+            data += [row]
+    except:  # noqa: E722
+        print("edge not updated due to Exception")
+        data += [row]
+        continue
+
+gpd.GeoDataFrame(data, crs=28992).to_file("rare_edges.gpkg")
+
+# %% write model
+# model.write(ribasim_toml)

notebooks/de_dommel/03_fix_basin_area.py

@@ -0,0 +1,81 @@
+# %%
+import geopandas as gpd
+import pandas as pd
+from ribasim_nl import CloudStorage, Model
+from shapely.geometry import MultiPolygon
+
+cloud = CloudStorage()
+
+
+# %% load model
+ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel", "model.toml")
+model = Model.read(ribasim_toml)
+
+
+area_gdf = gpd.read_file(cloud.joinpath("DeDommel", "verwerkt", "watervlakken", "LWW_2023_A_water_vlak_V.shp"))
+
+
+named_area_gdf = area_gdf.dissolve("NAAM", dropna=True).reset_index()
+unnamed_area_poly = area_gdf[area_gdf.NAAM.isna()].buffer(0.1).union_all().buffer(-0.1)
+unnamed_area_gdf = gpd.GeoDataFrame(geometry=list(unnamed_area_poly.geoms), crs=area_gdf.crs)
+
+dissolved_area_gdf = pd.concat([unnamed_area_gdf, named_area_gdf])
+dissolved_area_gdf.to_file(cloud.joinpath("DeDommel", "verwerkt", "water_area.gpkg"))
+
+# %%
+basin_df = model.basin.node.df.set_index("node_id")
+basin_area_df = gpd.read_file(
+    cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg"), engine="pyogrio", fid_as_index=True
+)
+basin_area_df.set_index("node_id", inplace=True)
+
+
+data = []
+ignore_basins = [1278, 1228, 1877, 1030]
+row = next(i for i in basin_df.itertuples() if i.Index == 1230)
+for row in basin_df.itertuples():
+    if row.Index not in ignore_basins:
+        edges_mask = (model.edge.df.from_node_id == row.Index) | (model.edge.df.to_node_id == row.Index)
+        edges_geom = model.edge.df[edges_mask].union_all()
+
+        selected_areas = dissolved_area_gdf[dissolved_area_gdf.intersects(edges_geom)]
+
+        basin_geom = gpd.clip(selected_areas, basin_area_df.at[row.Index, "geometry"]).union_all()
+        if isinstance(basin_geom, MultiPolygon):
+            basin_geom = MultiPolygon(i for i in basin_geom.geoms if i.area > 100)
+
+        data += [{"node_id": row.Index, "geometry": basin_geom}]
+
+        # set name to basin if empty
+        name = ""
+        area_df = selected_areas[selected_areas.contains(row.geometry)]
+        if not area_df.empty:
+            name = area_df.iloc[0].NAAM
+        model.basin.node.df.loc[model.basin.node[row.Index].index[0], ["name"]] = name
+        # assign name to edges if defined
+        model.edge.df.loc[edges_mask, ["name"]] = name
+
+# Manuals
+for node_id, name in [(1228, "Beatrixkanaal"), (1278, "Eindhovens kanaal")]:
+    data += [{"node_id": node_id, "geometry": dissolved_area_gdf[dissolved_area_gdf.NAAM == name].iloc[0].geometry}]
+    edges_mask = (model.edge.df.from_node_id == node_id) | (model.edge.df.to_node_id == node_id)
+    model.edge.df.loc[edges_mask, "name"] = name
+
+node_id, name = (1030, "Reusel")
+data += [
+    {
+        "node_id": node_id,
+        "geometry": dissolved_area_gdf[dissolved_area_gdf.NAAM == name]
+        .clip(basin_area_df.loc[[1015, 1030]].union_all())
+        .union_all(),
+    }
+]
+edges_mask = (model.edge.df.from_node_id == node_id) | (model.edge.df.to_node_id == node_id)
+model.edge.df.loc[edges_mask, "name"] = name
+
+area_df = gpd.GeoDataFrame(data, crs=model.basin.node.df.crs)
+area_df = area_df[~area_df.is_empty]
+model.basin.area.df = area_df
+# %%
+
+model.write(ribasim_toml)

src/ribasim_nl/ribasim_nl/model.py

@@ -5,7 +5,7 @@
 from pydantic import BaseModel
 from ribasim import Model, Node
 from ribasim.geometry.edge import NodeData
-from shapely.geometry import Point
+from shapely.geometry import LineString, Point
 from shapely.geometry.base import BaseGeometry
 
 from ribasim_nl.case_conversions import pascal_to_snake_case
@@ -308,3 +308,14 @@ def fix_unassigned_basin_area(self, method: str = "within", distance: float = 10
                 raise ValueError("Assign Basin Area to your model first")
         else:
             raise ValueError("Assign a Basin Node to your model first")
+
+    def reset_edge_geometry(self, edge_ids: list | None = None):
+        node_df = self.node_table().df.set_index("node_id")
+        if edge_ids is not None:
+            df = self.edge.df[self.edge.df.index.isin(edge_ids)]
+        else:
+            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"]])
+            self.edge.df.loc[row.Index, ["geometry"]] = geometry

src/ribasim_nl/ribasim_nl/network.py

@@ -442,7 +442,7 @@ def move_node(
             )
             return None
 
-    def add_node(self, point: Point, max_distance: float):
+    def add_node(self, point: Point, max_distance: float, align_distance: float = 100):
         # set _graph undirected to None
         self._graph_undirected = None
 
@@ -469,7 +469,7 @@ def add_node(self, point: Point, max_distance: float):
             self.add_link(node_from, node_id, us_geometry)
             self.add_link(node_id, node_to, ds_geometry)
 
-            return self.move_node(point, max_distance=max_distance, align_distance=100)
+            return self.move_node(point, max_distance=max_distance, align_distance=align_distance)
         else:
             logger.warning(
                 f"No Node added. Closest edge: {edge_id}, distance > max_distance ({edge_distance} > {max_distance})"