From d1f8ea135f44e36fd5b752d0a1a7cefa457c683f Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Tue, 20 Aug 2024 20:22:17 +0200 Subject: [PATCH 1/2] model with edges --- .../{get_model.py => 00_get_model.py} | 0 .../{get_verwerkt.py => 00_get_verwerkt.py} | 0 ...odify_model.py => 01_fix_model_network.py} | 0 notebooks/de_dommel/02_parameterize_model.py | 61 +++++++++++++++++++ src/ribasim_nl/ribasim_nl/model.py | 13 +++- src/ribasim_nl/ribasim_nl/network.py | 4 +- 6 files changed, 75 insertions(+), 3 deletions(-) rename notebooks/de_dommel/{get_model.py => 00_get_model.py} (100%) rename notebooks/de_dommel/{get_verwerkt.py => 00_get_verwerkt.py} (100%) rename notebooks/de_dommel/{modify_model.py => 01_fix_model_network.py} (100%) create mode 100644 notebooks/de_dommel/02_parameterize_model.py diff --git a/notebooks/de_dommel/get_model.py b/notebooks/de_dommel/00_get_model.py similarity index 100% rename from notebooks/de_dommel/get_model.py rename to notebooks/de_dommel/00_get_model.py diff --git a/notebooks/de_dommel/get_verwerkt.py b/notebooks/de_dommel/00_get_verwerkt.py similarity index 100% rename from notebooks/de_dommel/get_verwerkt.py rename to notebooks/de_dommel/00_get_verwerkt.py diff --git a/notebooks/de_dommel/modify_model.py b/notebooks/de_dommel/01_fix_model_network.py similarity index 100% rename from notebooks/de_dommel/modify_model.py rename to notebooks/de_dommel/01_fix_model_network.py diff --git a/notebooks/de_dommel/02_parameterize_model.py b/notebooks/de_dommel/02_parameterize_model.py new file mode 100644 index 0000000..fc37cf7 --- /dev/null +++ b/notebooks/de_dommel/02_parameterize_model.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) diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index b1d45a5..cd9cfc7 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/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 diff --git a/src/ribasim_nl/ribasim_nl/network.py b/src/ribasim_nl/ribasim_nl/network.py index 0e991d9..6e62809 100644 --- a/src/ribasim_nl/ribasim_nl/network.py +++ b/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})" From 410895f674b71c055646442b7f1df4c661117181 Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Wed, 21 Aug 2024 17:11:30 +0200 Subject: [PATCH 2/2] fix edges and basins --- ..._parameterize_model.py => 02_fix_edges.py} | 0 notebooks/de_dommel/03_fix_basin_area.py | 81 +++++++++++++++++++ 2 files changed, 81 insertions(+) rename notebooks/de_dommel/{02_parameterize_model.py => 02_fix_edges.py} (100%) create mode 100644 notebooks/de_dommel/03_fix_basin_area.py diff --git a/notebooks/de_dommel/02_parameterize_model.py b/notebooks/de_dommel/02_fix_edges.py similarity index 100% rename from notebooks/de_dommel/02_parameterize_model.py rename to notebooks/de_dommel/02_fix_edges.py diff --git a/notebooks/de_dommel/03_fix_basin_area.py b/notebooks/de_dommel/03_fix_basin_area.py new file mode 100644 index 0000000..504362f --- /dev/null +++ b/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)