From 848eed8f221a0ccf826fd2e3be6a2fb3ab6451f9 Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Thu, 12 Dec 2024 10:27:26 +0100 Subject: [PATCH 1/5] fix_graph (#202) --- src/ribasim_nl/ribasim_nl/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index 5510fea..9c4a739 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -754,7 +754,7 @@ def merge_basins( ) 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] + paths = [i for i in nx.all_shortest_paths(nx.Graph(self.graph), node_id, to_node_id) if len(i) == 3] if len(paths) == 0: raise ValueError(f"basin {node_id} not a direct neighbor of basin {to_node_id}") From 1508c4b61e1239703a5c63aead8714ec6acbe44e Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Thu, 12 Dec 2024 12:10:40 +0100 Subject: [PATCH 2/5] merge_basins (#203) --- src/ribasim_nl/ribasim_nl/model.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index 9c4a739..fe47d8b 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -754,6 +754,7 @@ def merge_basins( ) if are_connected and (to_node_type != "FlowBoundary"): + self._graph = None # set self._graph to None, so it will regenerate on currend edge-table paths = [i for i in nx.all_shortest_paths(nx.Graph(self.graph), node_id, to_node_id) if len(i) == 3] if len(paths) == 0: From 2dd1288343215487d4ef9c0439c05591a584323e Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Thu, 12 Dec 2024 13:00:04 +0100 Subject: [PATCH 3/5] Wdod (#201) edit-acties voor WDOD --------- Co-authored-by: ngoorden <49673619+ngoorden@users.noreply.github.com> Co-authored-by: Martijn Visser --- .../drents_overijsselse_delta/00_get_model.py | 2 +- .../01_fix_model_network.py | 29 ++++++++++--------- src/ribasim_nl/ribasim_nl/model.py | 10 +++---- 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/notebooks/drents_overijsselse_delta/00_get_model.py b/notebooks/drents_overijsselse_delta/00_get_model.py index 2cf928b..912d599 100644 --- a/notebooks/drents_overijsselse_delta/00_get_model.py +++ b/notebooks/drents_overijsselse_delta/00_get_model.py @@ -12,4 +12,4 @@ 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")) + ribasim_toml.replace(ribasim_toml.with_name(f"{short_name}.toml")) diff --git a/notebooks/drents_overijsselse_delta/01_fix_model_network.py b/notebooks/drents_overijsselse_delta/01_fix_model_network.py index 5ec802b..645d001 100644 --- a/notebooks/drents_overijsselse_delta/01_fix_model_network.py +++ b/notebooks/drents_overijsselse_delta/01_fix_model_network.py @@ -25,20 +25,22 @@ layer="duikersifonhevel", ) -split_line_gdf = gpd.read_file( - cloud.joinpath(authority, "verwerkt", "fix_user_data.gpkg"), layer="split_basins", fid_as_index=True -) - # Load node edit data model_edits_url = cloud.joinurl(authority, "verwerkt", "model_edits.gpkg") model_edits_path = cloud.joinpath(authority, "verwerkt", "model_edits.gpkg") if not model_edits_path.exists(): cloud.download_file(model_edits_url) +# Load node edit data +fix_user_data_url = cloud.joinurl(authority, "verwerkt", "fix_user_data.gpkg") +fix_user_data_path = cloud.joinpath(authority, "verwerkt", "fix_user_data.gpkg") +if not fix_user_data_path.exists(): + cloud.download_file(fix_user_data_url) + +split_line_gdf = gpd.read_file( + cloud.joinpath(authority, "verwerkt", fix_user_data_path), 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) @@ -339,15 +341,16 @@ "remove_basin_area", "split_basin", "merge_basins", + "add_basin", "update_node", "add_basin_area", - "add_basin", "update_basin_area", "redirect_edge", "reverse_edge", "deactivate_node", "move_node", "remove_node", + "connect_basins", ] actions = [i for i in actions if i in gpd.list_layers(model_edits_path).name.to_list()] @@ -364,16 +367,16 @@ kwargs = {k: v for k, v in row._asdict().items() if k in keywords} method(**kwargs) +# remove unassigned basin area +model.fix_unassigned_basin_area() +model.remove_unassigned_basin_area() -# %% Reset static tables - -# Reset static tables model = reset_static_tables(model) - # %% write model model.use_validation = True model.write(ribasim_toml) model.invalid_topology_at_node().to_file(ribasim_toml.with_name("invalid_topology_at_connector_nodes.gpkg")) - +model.report_basin_area() +model.report_internal_basins() # %% diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index fe47d8b..2a90279 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -785,17 +785,17 @@ def merge_basins( 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 isinstance(poly, Polygon): + poly = MultiPolygon([poly]) + # if to_node_id has area we union both areas - if to_node_id in self.basin.area.df.node_id.to_numpy(): + if len(self.basin.area.df.loc[self.basin.area.df.node_id == to_node_id]) == 1: 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, From 030524bdc0f07362a564867b4ab7ad6ced7eb8b1 Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Fri, 20 Dec 2024 11:02:47 +0100 Subject: [PATCH 4/5] Hea fixes (#208) Co-authored-by: ngoorden <49673619+ngoorden@users.noreply.github.com> --- .../hunze_en_aas/01_fix_model_network.py | 54 ++++++++++++++++++- notebooks/upload_feedback_formulieren.py | 10 ---- src/ribasim_nl/ribasim_nl/model.py | 5 +- src/ribasim_nl/ribasim_nl/run_model.py | 47 ++++++++++++++++ 4 files changed, 103 insertions(+), 13 deletions(-) create mode 100644 src/ribasim_nl/ribasim_nl/run_model.py diff --git a/notebooks/hunze_en_aas/01_fix_model_network.py b/notebooks/hunze_en_aas/01_fix_model_network.py index 39b956b..f13df40 100644 --- a/notebooks/hunze_en_aas/01_fix_model_network.py +++ b/notebooks/hunze_en_aas/01_fix_model_network.py @@ -27,6 +27,10 @@ if not model_edits_path.exists(): cloud.download_file(model_edits_url) +# Load area file to fill basin area holes +ribasim_areas_path = cloud.joinpath(authority, "verwerkt", "4_ribasim", "areas.gpkg") +ribasim_areas_gdf = gpd.read_file(ribasim_areas_path, fid_as_index=True, layer="areas") + # %% some stuff we'll need again manning_data = manning_resistance.Static(length=[100], manning_n=[0.04], profile_width=[10], profile_slope=[1]) @@ -90,18 +94,26 @@ # Reset static tables model = reset_static_tables(model) +# fix unassigned basin area +model.fix_unassigned_basin_area() +model.explode_basin_area() +# fix unassigned basin area +model.fix_unassigned_basin_area() # %% + actions = [ + "remove_basin_area", "remove_node", "add_basin", + "update_node", "add_basin_area", "update_basin_area", - "merge_basins", "reverse_edge", + "redirect_edge", + "merge_basins", "move_node", "connect_basins", - "update_node", "deactivate_node", ] actions = [i for i in actions if i in gpd.list_layers(model_edits_path).name.to_list()] @@ -117,9 +129,47 @@ method(**kwargs) +# %% 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.remove_unassigned_basin_area() + # %% write model model.use_validation = True model.write(ribasim_toml) model.invalid_topology_at_node().to_file(ribasim_toml.with_name("invalid_topology_at_connector_nodes.gpkg")) +model.report_basin_area() +model.report_internal_basins() # %% diff --git a/notebooks/upload_feedback_formulieren.py b/notebooks/upload_feedback_formulieren.py index a25af36..2ffe0a9 100644 --- a/notebooks/upload_feedback_formulieren.py +++ b/notebooks/upload_feedback_formulieren.py @@ -8,17 +8,7 @@ cloud = CloudStorage() WATER_AUTHORITIES = [ - "AaenMaas", - "BrabantseDelta", - "DeDommel", "DrentsOverijsselseDelta", - "HunzeenAas", - "Limburg", - "Noorderzijlvest", - "RijnenIJssel", - "StichtseRijnlanden", - "ValleienVeluwe", - "Vechtstromen", ] FEEDBACK_XLS = cloud.joinpath("Basisgegevens", "feedbackformulier", "Feedback Formulier.xlsx") diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index 2a90279..0bb0159 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -713,7 +713,7 @@ def remove_unassigned_basin_area(self): 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 + self.basin.area.df = df def explode_basin_area(self, remove_z=True): df = self.basin.area.df.explode().reset_index(drop=True) @@ -891,3 +891,6 @@ def invalid_topology_at_node(self, edge_type: str = "flow") -> gpd.GeoDataFrame: return gpd.GeoDataFrame( [], columns=["node_id", "node_type", "exception"], geometry=gpd.GeoSeries(crs=self.crs) ).set_index("node_id") + + +# %% diff --git a/src/ribasim_nl/ribasim_nl/run_model.py b/src/ribasim_nl/ribasim_nl/run_model.py new file mode 100644 index 0000000..a750e71 --- /dev/null +++ b/src/ribasim_nl/ribasim_nl/run_model.py @@ -0,0 +1,47 @@ +import os +import subprocess +from pathlib import Path + +# TODO: add ribasim_exe so it can be used if ribasim is not part of env path +# TODO: check if ribasim is in path, stop if not and ribasim_exe is not provided +# TODO: raise FileNotFoundError if toml_path does not exist. User + + +def run( + toml_path: Path, + stream_output: bool = True, + returncode: bool = True, +): + """To run a Ribasim model + + Args: + toml_path (Path): path to your ribasim toml-file + stream_output (bool, optional): stream output in IDE. Defaults to False. + returncode (bool, optional): return return code after running model. Defaults to True. + + """ + env = os.environ.copy() + + input = "" + proc = subprocess.Popen( + ["ribasim", toml_path.as_posix()], + cwd=toml_path.parent.as_posix(), + env=env, + stdin=subprocess.PIPE, + stdout=subprocess.PIPE, + encoding="ascii", + ) + if stream_output: + with proc: + proc.stdin.write(input) + proc.stdin.close() + for line in proc.stdout: + print(line, end="") + outs = None + else: + outs, _ = proc.communicate(input) + + if returncode: + return proc.returncode + else: + return outs From bc0aaea391f0ca4959eb2966e3d1bc780ecd44dc Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Fri, 20 Dec 2024 11:04:10 +0100 Subject: [PATCH 5/5] fix_nzv (#209) Co-authored-by: Martijn Visser --- .../noorderzijlvest/01b_fix_basin_area.py | 62 +++++++++++++++---- src/ribasim_nl/ribasim_nl/geometry.py | 59 +++++++++++++++++- src/ribasim_nl/ribasim_nl/network.py | 11 ++-- 3 files changed, 113 insertions(+), 19 deletions(-) diff --git a/notebooks/noorderzijlvest/01b_fix_basin_area.py b/notebooks/noorderzijlvest/01b_fix_basin_area.py index a958908..2ca69ea 100644 --- a/notebooks/noorderzijlvest/01b_fix_basin_area.py +++ b/notebooks/noorderzijlvest/01b_fix_basin_area.py @@ -5,6 +5,7 @@ import pandas as pd from networkx import all_shortest_paths, shortest_path from shapely.geometry import MultiLineString +from shapely.ops import snap, split from ribasim_nl import CloudStorage, Model, Network @@ -32,9 +33,27 @@ cloud_storage.joinpath( authority_name, "verwerkt", "5_D_HYDRO_export", "hydroobjecten", "Noorderzijlvest_hydroobjecten.shp" ), + use_fid_as_indes=True, ) + +points = ( + model.node_table().df[model.node_table().df.node_type.isin(["TabulatedRatingCurve", "Outlet", "Pump"])].geometry +) + +for row in lines_gdf.itertuples(): + line = row.geometry + snap_points = points[points.distance(line) < 0.1] + snap_points = snap_points[snap_points.distance(line.boundary) > 0.1] + if not snap_points.empty: + snap_point = snap_points.union_all() + line = snap(line, snap_point, 1e-8) + split_lines = split(line, snap_point) + + lines_gdf.loc[row.Index, ["geometry"]] = split_lines + +lines_gdf = lines_gdf.explode(index_parts=False, ignore_index=True) lines_gdf.crs = 28992 -network = Network(lines_gdf) +network = Network(lines_gdf.copy()) network.to_file(cloud_storage.joinpath(authority_name, "verwerkt", "network.gpkg")) # %% add snap_point to he_df @@ -64,6 +83,9 @@ he_outlet_df.set_index("HEIDENT", inplace=True) he_df.set_index("HEIDENT", inplace=True) +# niet altijd ligt de coordinaat goed +he_outlet_df.loc["GPGKST0470", ["geometry"]] = model.manning_resistance[892].geometry + he_outlet_df.to_file(cloud_storage.joinpath(authority_name, "verwerkt", "HydrologischeEenheden_v45_outlets.gpkg")) # %% Edit network @@ -78,7 +100,16 @@ cloud_storage.download_file(model_edits_url) -for action in gpd.list_layers(model_edits_path).name.to_list(): +for action in [ + "merge_basins", + "remove_node", + "update_node", + "reverse_edge", + "connect_basins", + "move_node", + "add_basin", + "remove_edge", +]: print(action) # get method and args method = getattr(model, action) @@ -118,9 +149,7 @@ def find_basin_id(kwk_code): def get_network_node(point): node = network.move_node(point, max_distance=1, align_distance=10) if node is None: - node = network.add_node(point, max_distance=1) - if node is None: - node = network.nodes.distance(point).idxmin() + node = network.add_node(point, max_distance=1, align_distance=10) return node @@ -184,7 +213,7 @@ def get_network_node(point): ) model.edge.df.loc[mask, ["geometry"]] = edge - mask = he_df.node_id.isna() & (he_outlet_df.distance(MultiLineString(data)) < 1) + mask = he_df.node_id.isna() & (he_outlet_df.distance(MultiLineString(data)) < 0.75) he_df.loc[mask, ["node_id"]] = node_id # %% add last missings @@ -201,23 +230,34 @@ def get_network_node(point): he_df.loc[row.Index, ["node_id"]] = basin_node_id -# %% renew model Basin Area - -# based on actions above we update the Basin / Area table: -# basins we found on KWKuit defined in the he-table - data = [] for node_id, df in he_df[he_df["node_id"].notna()].groupby("node_id"): geometry = df.union_all() streefpeil = df["OPVAFWZP"].min() data += [{"node_id": node_id, "meta_streefpeil": streefpeil, "geometry": geometry}] + df = gpd.GeoDataFrame(data, crs=model.crs) df.loc[:, "geometry"] = df.buffer(0.1).buffer(-0.1) df.index.name = "fid" model.basin.area.df = df +for action in ["remove_basin_area", "add_basin_area"]: + 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) + +model.remove_unassigned_basin_area() + # %% + 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/src/ribasim_nl/ribasim_nl/geometry.py b/src/ribasim_nl/ribasim_nl/geometry.py index 439e374..4f7122c 100644 --- a/src/ribasim_nl/ribasim_nl/geometry.py +++ b/src/ribasim_nl/ribasim_nl/geometry.py @@ -1,9 +1,10 @@ +# %% """Functions to apply on a shapely.geometry""" from typing import get_type_hints -from shapely.geometry import LineString, MultiPolygon, Point, Polygon -from shapely.ops import polygonize, polylabel +from shapely.geometry import LineString, MultiLineString, MultiPolygon, Point, Polygon +from shapely.ops import polygonize, polylabel, snap, split from ribasim_nl.generic import _validate_inputs @@ -193,3 +194,57 @@ def edge(point_from: Point, point_to: Point) -> LineString: if point_to.has_z: point_to = Point(point_to.x, point_to.y) return LineString((point_from, point_to)) + + +def project_point(line: LineString, point: Point, tolerance: float = 1) -> Point: + """Projects a point to a LineString + + Args: + line (LineString): LineString to project point on + point (Point): Point to project on LineString + tolerance (float): Tolerance of point to line + + Returns + ------- + Point: Projected Point + """ + if point.distance(line) > tolerance: + raise ValueError(f"point > {tolerance} from line ({point.distance(line)})") + return line.interpolate(line.project(point)) + + +def split_line(line: LineString, point: Point, tolerance: float = 0.1) -> MultiLineString | LineString: + """Split a line into a 2 geometry multiline + + Args: + line (LineString): input line + point (Point): point to split on + tolerance (float, optional): tolerance of point to line. Defaults to 0.1. + + Returns + ------- + MultiLineString | LineString: in case LineString is split, return MultiLineString with 2 LineStrings + """ + # if point is within tolerance of line.boundary we don't split + if line.boundary.distance(point) < tolerance: + return line + + # try if a normal split works + result = split(line, point) + if len(list(result.geoms)) == 2: + return MultiLineString(result) + + # if not, we try it again + else: + # project the point on the line, checking if it's not too far of first + point = project_point(line, point, tolerance) + + # we snap the line to the point so it will have the point coordinate + line = snap(line, point, 1e-8) + + # now we should be able to split + result = split(line, point) + if len(list(result.geoms)) == 2: + return MultiLineString(result) + else: + return line diff --git a/src/ribasim_nl/ribasim_nl/network.py b/src/ribasim_nl/ribasim_nl/network.py index 138df43..90eb27b 100644 --- a/src/ribasim_nl/ribasim_nl/network.py +++ b/src/ribasim_nl/ribasim_nl/network.py @@ -9,10 +9,10 @@ import pandas as pd from geopandas import GeoDataFrame, GeoSeries from networkx import DiGraph, Graph, NetworkXNoPath, shortest_path, traversal -from shapely import force_2d from shapely.geometry import LineString, Point, box from shapely.ops import snap, split +from ribasim_nl.geometry import drop_z, split_line from ribasim_nl.styles import add_styles_to_geopackage logger = logging.getLogger(__name__) @@ -93,9 +93,8 @@ def validate_inputs(self): self.lines_gdf = self.lines_gdf.explode(index_parts=False) # remove z-coordinates - self.lines_gdf.loc[:, "geometry"] = gpd.GeoSeries( - force_2d(self.lines_gdf.geometry.array), crs=self.lines_gdf.crs - ) + if self.lines_gdf.has_z.any(): + self.lines_gdf.loc[:, "geometry"] = self.lines_gdf.geometry.apply(lambda x: drop_z(x) if x.has_z else x) @classmethod def from_lines_gpkg(cls, gpkg_file: str | Path, layer: str | None = None, **kwargs): @@ -471,7 +470,7 @@ def add_node(self, point: Point, max_distance: float, align_distance: float = 10 # add edges self.graph.remove_edge(node_from, node_to) - us_geometry, ds_geometry = split(snap(edge_geometry, node_geometry, 0.01), node_geometry).geoms + us_geometry, ds_geometry = split_line(edge_geometry, node_geometry).geoms self.add_link(node_from, node_id, us_geometry) self.add_link(node_id, node_to, ds_geometry) @@ -650,7 +649,7 @@ def to_file(self, path: str | Path): raise ValueError(f"{path} is not a GeoPackage, please provide a file with extention 'gpkg'") # write nodes and links - self.nodes.to_file(path, layer="nodes", engine="pyogrio") + self.nodes[["geometry", "type"]].to_file(path, layer="nodes", engine="pyogrio") self.links.to_file(path, layer="links", engine="pyogrio") # add styles add_styles_to_geopackage(path)