fix_nzv (#209)

Co-authored-by: Martijn Visser <[email protected]>

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()
 # %%

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

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)