direct basins to ribasim_nl

notebooks/rijkswaterstaat/netwerk.py

@@ -1,7 +1,7 @@
 # %%
 import geopandas as gpd
 from ribasim_nl import CloudStorage
-from ribasim_nl.geodataframe import join_by_poly_overlay, split_basins
+from ribasim_nl.geodataframe import direct_basins, join_by_poly_overlay, split_basins
 
 cloud = CloudStorage()
 
@@ -53,4 +53,32 @@
 
 rws_krw_line_gdf.to_file(rws_krw_lines)
 
-# %%
+# %% direct basins
+
+basin_ident = "owmident"
+link_ident = "Name"
+
+basins_gdf = gpd.read_file(
+    cloud.joinpath("Rijkswaterstaat", "verwerkt", "krw_basins_vlakken.gpkg")
+)
+
+network_gdf = gpd.read_file(
+    cloud.joinpath("Basisgegevens", "lsm3-j18_5v6", "shapes", "network_Branches.shp")
+)
+network_gdf.set_crs(28992, inplace=True)
+drop_duplicates = True
+
+poly_directions_gdf = direct_basins(basins_gdf, network_gdf, basin_ident, link_ident)
+
+
+poly_directions_gdf.to_file(
+    cloud.joinpath("Rijkswaterstaat", "verwerkt", "krw_basins_verbindingen.gpkg")
+)
+
+# %% snap nodes
+
+# %% build graph
+
+# %% build_network
+
+# %% A(h)

src/ribasim_nl/ribasim_nl/geodataframe.py

@@ -1,3 +1,4 @@
+"""All functions resulting in a geopandas.GeoDataFrame"""
 from typing import Literal, Union
 
 import geopandas as gpd
@@ -114,3 +115,129 @@ def split_basins(basins_gdf: GeoDataFrame, lines_gdf: GeoDataFrame) -> GeoDataFr
             ignore_index=True,
         )
     return basins_gdf
+
+
+def direct_basins(
+    basins_gdf: GeoDataFrame,
+    network_gdf: GeoDataFrame,
+    basin_ident: str,
+    link_ident: str,
+    drop_duplicates: bool = True,
+) -> GeoDataFrame:
+    """Find directions of basins 'basins_gdf' (polygons) by 'network_gdf'(linestrings)
+
+    Parameters
+    ----------
+    basins_gdf : GeoDataFrame
+        GeoDataFrame with basins
+    network_gdf : GeoDataFrame
+        GeoDataFrame with network trough basins
+    basin_ident : str
+        column in basins_gdf to with unique values ending up in the result as 'us_basin' and 'ds_basin'
+    link_ident : str
+        column in network_gdf with unique values ending up in the result as 'link_ident'
+    drop_duplicates : bool, optional
+        Duplicate directions, multiple links between two basins will be removed, by default True
+
+    Returns
+    -------
+    GeoDataFrame
+        Selection of network_gdf representing links between basins in basins_gdf
+    """
+
+    def find_intersecting_basins(line):
+        """Find intersecting basins on a LineString"""
+        intersecting_basins_gdf = basins_gdf.iloc[
+            basins_gdf.sindex.intersection(line.bounds)
+        ]
+        intersecting_basins_gdf = intersecting_basins_gdf[
+            intersecting_basins_gdf["geometry"].intersects(line)
+        ]
+        return intersecting_basins_gdf[basin_ident]
+
+    def find_ident(point):
+        """Find basin containing a Point"""
+        # filter by spatial index
+        containing_basins_gdf = basins_gdf.iloc[
+            basins_gdf.sindex.intersection(point.bounds)
+        ]
+
+        # filter by contain and take first row
+        containing_basins_gdf = containing_basins_gdf[
+            containing_basins_gdf["geometry"].contains(point)
+        ]
+
+        if containing_basins_gdf.empty:
+            return None
+        elif len(containing_basins_gdf) == 1:
+            return getattr(containing_basins_gdf.iloc[0], basin_ident)
+        else:
+            print(containing_basins_gdf)
+
+    # list of dicts to be filled with 'us_basin', 'ds_basin' and 'geometry'
+    data = []
+
+    # we iterate per basin to determine intersections
+    for poly_row in basins_gdf.itertuples():
+        # identification for the basin
+        poly_ident = getattr(poly_row, basin_ident)
+
+        ## select intersecting lines, use spatial index (and then exact intersection)
+        intersecting_network_gdf = network_gdf.iloc[
+            network_gdf.sindex.intersection(poly_row.geometry.bounds)
+        ]
+        intersecting_network_gdf = intersecting_network_gdf[
+            intersecting_network_gdf["geometry"].intersects(poly_row.geometry.boundary)
+        ]
+
+        # iterate per link, to find basin directions
+        for line_row in intersecting_network_gdf.itertuples():
+            # find upstream and downstream points
+            us_point, ds_point = line_row.geometry.boundary.geoms
+
+            # find upstream and downstream basins
+            us_ident = find_ident(us_point)
+            ds_ident = find_ident(ds_point)
+
+            # determine if/how a link should be added to data
+            if ((us_ident is not None) and (ds_ident is not None)) and (
+                us_ident != ds_ident
+            ):
+                link_id = getattr(line_row, link_ident)
+
+                if poly_ident not in [us_ident, ds_ident]:
+                    data += [
+                        {
+                            "link_ident": link_id,
+                            "us_basin": us_ident,
+                            "ds_basin": poly_ident,
+                            "geometry": line_row.geometry,
+                        },
+                        {
+                            "link_ident": link_id,
+                            "us_basin": poly_ident,
+                            "ds_basin": ds_ident,
+                            "geometry": line_row.geometry,
+                        },
+                    ]
+                elif len(find_intersecting_basins(line_row.geometry)) == 2:
+                    data += [
+                        {
+                            "link_ident": link_id,
+                            "us_basin": us_ident,
+                            "ds_basin": ds_ident,
+                            "geometry": line_row.geometry,
+                        }
+                    ]
+                else:
+                    print(f"we skip a case for link: '{link_id}'")
+
+    poly_directions_gdf = gpd.GeoDataFrame(data, crs=basins_gdf.crs)
+
+    # if we found duplicates, we may want to drop them
+    if drop_duplicates:
+        poly_directions_gdf.drop_duplicates(
+            ["us_basin", "ds_basin"], keep="first", inplace=True
+        )
+
+    return poly_directions_gdf

src/ribasim_nl/ribasim_nl/geometry.py

@@ -1,4 +1,3 @@
-# %%
 """Functions to apply on a shapely.geometry"""
 from typing import List, Union, get_type_hints