diff --git a/notebooks/noorderzijlvest/01b_fix_basin_area.py b/notebooks/noorderzijlvest/01b_fix_basin_area.py index a958908..4abf93c 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,38 @@ 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 +# # split lines at points +# for point in points: +# line_idx = lines_gdf.distance(point).idxmin() +# line = lines_gdf.at[line_idx, "geometry"] +# if line.has_z: +# line = drop_z(line) +# line = split_line(line=line, point=point, tolerance=0.1) +# lines_gdf.loc[line_idx, ["geometry"]] = line + +lines_gdf = lines_gdf.explode(index_parts=False, ignore_index=True) lines_gdf.crs = 28992 -network = Network(lines_gdf) +lines_gdf.to_file("lines_out.gpkg", index=True, fid="fid") + +network = Network(lines_gdf.copy()) +network.lines_gdf.to_file("lines_out2.gpkg") network.to_file(cloud_storage.joinpath(authority_name, "verwerkt", "network.gpkg")) # %% add snap_point to he_df @@ -64,6 +94,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 @@ -118,9 +151,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 +215,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,18 +232,15 @@ 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 = pd.concat( + [gpd.GeoDataFrame(data, crs=model.crs), gpd.read_file(model_edits_path, layer="add_basin_area")], ignore_index=True +) df.loc[:, "geometry"] = df.buffer(0.1).buffer(-0.1) df.index.name = "fid" model.basin.area.df = df @@ -220,4 +248,5 @@ def get_network_node(point): # %% 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)