From bc0aaea391f0ca4959eb2966e3d1bc780ecd44dc Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Fri, 20 Dec 2024 11:04:10 +0100 Subject: [PATCH] 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 a958908f..2ca69ea6 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 439e374d..4f7122c0 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 138df43c..90eb27b4 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)