IJsselmeer model

environment.yml

@@ -14,7 +14,7 @@ dependencies:
   - matplotlib
   - mypy
   - openpyxl
-  - pandas!=2.1.0
+  - pandas<2.1.0
   - pandas-stubs
   - pandera
   - pip
@@ -27,6 +27,7 @@ dependencies:
   - pytest
   - pytest-cov
   - python>=3.9
+  - rastersats
   - ruff
   - shapely>=2.0
   - tomli

notebooks/ijsselmeermodel/basins.py

@@ -0,0 +1,140 @@
+# %%
+import os
+from pathlib import Path
+
+import geopandas as gpd
+import pandas as pd
+import ribasim
+from hydamo import code_utils
+from ribasim_nl.utils.geometry import cut_basin, drop_z
+from ribasim_nl.utils.geoseries import basins_to_points
+
+DATA_DIR = Path(os.getenv("RIBASIM_NL_DATA_DIR"))
+MODEL_DIR = Path(os.getenv("RIBASIM_NL_MODEL_DIR")) / "ijsselmeer"
+
+DEFAULT_AREA = [0.01, 1000.0]
+DEFAULT_LEVEL = [0.0, 1.0]
+DEFAULT_EVAPORATION = 0
+DEFAULT_PRECIPITATION = 0.002 / 86400
+START_NODE_ID = 1
+basins = []
+
+
+def add_basin(**kwargs):
+    global basins
+    kwargs["geometry"] = drop_z(kwargs["geometry"])
+    basins += [kwargs]
+
+
+krw_ids = [
+    "NL92_IJSSELMEER",
+    "NL92_MARKERMEER",
+    "NL92_RANDMEREN_ZUID",
+    "NL92_RANDMEREN_OOST",
+    "NL92_KETELMEER_VOSSEMEER",
+    "NL92_ZWARTEMEER",
+]
+
+rws_krw_gpkg = DATA_DIR / r"KRW/krw-oppervlaktewaterlichamen-nederland-vlakken.gpkg"
+rws_krw_gdf = gpd.read_file(rws_krw_gpkg).set_index("owmident")
+
+
+start_node_id = 1
+
+# rws_krw_gdf.loc[krw_ids].explore()
+
+krw_cutlines_gdf = gpd.read_file(MODEL_DIR / "model_data.gpkg", layer="krw_cutlines")
+
+
+def user_id(code, wbh_code, code_postfix=None):
+    code = code.split("_", 1)[-1]
+    if code_postfix:
+        code = f"{code}_{code_postfix }"
+    return code_utils.generate_model_id(code, "basin", wbh_code=wbh_code)
+
+
+for row in rws_krw_gdf.loc[krw_ids].itertuples():
+    # row = next(rws_krw_gdf.loc[krw_ids].itertuples())
+    code = row.Index
+    basin_polygon = row.geometry
+    if code in krw_cutlines_gdf.owmident.to_numpy():
+        if code in krw_cutlines_gdf.owmident.to_numpy():
+            if basin_polygon.geom_type == "Polygon":
+                for cut_line in (
+                    krw_cutlines_gdf[krw_cutlines_gdf.owmident == row.Index]
+                    .sort_values("cut_order")
+                    .itertuples()
+                ):
+                    # cut_line = krw_cutlines_gdf[krw_cutlines_gdf.owmident == row.Index].sort_values("cut_order").geometry[0]
+                    basin_multi_polygon = cut_basin(basin_polygon, cut_line.geometry)
+                    geometry = basin_multi_polygon.geoms[0]
+                    add_basin(
+                        user_id=user_id(code, "80", cut_line.cut_order),
+                        geometry=geometry,
+                        krw=True,
+                    )
+                    basin_polygon = basin_multi_polygon.geoms[1]
+                add_basin(
+                    user_id=user_id(code, "80", cut_line.cut_order + 1),
+                    geometry=basin_polygon,
+                    krw=True,
+                )
+            else:
+                raise TypeError(
+                    f"basin_polygon not of correct type {basin_polygon.geom_type}"
+                )
+    else:
+        add_basin(user_id=user_id(code, "80"), geometry=basin_polygon, lsw=True)
+
+gdf = gpd.read_file(DATA_DIR / r"Zuiderzeeland/Oplevering LHM/peilgebieden.gpkg")
+
+# Define your selection criteria
+mask = (
+    (gdf["GPGIDENT"] == "LVA.01")
+    | (gdf["GPGIDENT"] == "3.01")
+    | (gdf["GPGIDENT"] == "LAGE AFDELING")
+    | (gdf["GPGIDENT"] == "HOGE AFDELING")
+)
+
+# Process the selected polygons and add centroids and attributes
+for index, row in gdf[mask].iterrows():
+    add_basin(user_id=user_id(row.GPGIDENT, "37"), geometry=row.geometry)
+
+# Also, save the selected polygons to the new GeoPackage
+# "sel_peilgebieden_gdf.to_file(output_geopackage, driver="GPKG")
+
+
+basins_gdf = gpd.GeoDataFrame(basins, crs=28992)
+basins_gdf["node_id"] = basins_gdf.index + START_NODE_ID
+basins_gdf.to_file(MODEL_DIR / "basins.gpkg", layer="basins_areas")
+basins_gdf.loc[:, "geometry"] = basins_to_points(basins_gdf["geometry"])
+basins_gdf.to_file(MODEL_DIR / "basins.gpkg", layer="basins")
+
+
+## %% generate profiles
+data = []
+for row in basins_gdf.itertuples():
+    data += list(
+        zip(
+            [row.node_id] * len(DEFAULT_AREA),
+            [row.user_id] * len(DEFAULT_AREA),
+            DEFAULT_AREA,
+            DEFAULT_LEVEL,
+        )
+    )
+    profile_df = pd.DataFrame(data, columns=["node_id", "user_id", "area", "level"])
+
+## %% generate static
+static_df = basins_gdf[["node_id", "user_id"]].copy()
+static_df["drainage"] = 0
+static_df["potential_evaporation"] = DEFAULT_EVAPORATION
+static_df["infiltration"] = 0
+static_df["precipitation"] = DEFAULT_PRECIPITATION
+static_df["urban_runoff"] = 0
+
+## %%
+profile_df["remarks"] = profile_df["user_id"]
+static_df["remarks"] = static_df["user_id"]
+
+## %%
+basin = ribasim.Basin(profile=profile_df, static=static_df)

notebooks/ijsselmeermodel/brongegevens.ipynb

@@ -0,0 +1,44 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geopandas as gpd\n",
+    "from hydamo import HyDAMO\n",
+    "from pathlib import Path\n",
+    "import os"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "DATA_DIR = os.getenv(\"RIBASIM_NL_DATA_DIR\")\n",
+    "\n",
+    "# file-paths\n",
+    "kunstwerken_gpkg = Path(DATA_DIR) / \"nl_kunstwerken.gpkg\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kunstwerken_gdf = gpd.read_file(kunstwerken_gpkg)"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

notebooks/ijsselmeermodel/ijsselmeer_model.py

@@ -0,0 +1,84 @@
+import os
+from pathlib import Path
+
+import fiona
+import geopandas as gpd
+import pandas as pd
+
+DATA_DIR = Path(os.getenv("RIBASIM_NL_DATA_DIR"))
+MODEL_DIR = Path(os.getenv("RIBASIM_NL_MODEL_DIR")) / "ijsselmeer"
+
+
+# %% kunstwerken Zuiderzeeland
+
+
+# Inlezen waterschapsgrenzen
+shapefile = DATA_DIR / r"nederland/Waterschapsgrenzen.shp"
+waterschapsgrenzen_gdf = gpd.read_file(shapefile)
+
+# Filter the GeoDataFrame to select the feature with value 'Waterschap Zuiderzeeland'
+selected_waterschap_gdf = waterschapsgrenzen_gdf[
+    waterschapsgrenzen_gdf["waterschap"] == "Waterschap Zuiderzeeland"
+]
+
+output_gpkg = MODEL_DIR / "ZZL_grens.gpkg"
+selected_waterschap_gdf.to_file(output_gpkg, driver="GPKG")
+
+# Path to the GeoPackage file
+gpkg_path = DATA_DIR / r"uitlaten_inlaten.gpkg"
+
+# List available layers in the GeoPackage
+layers = fiona.listlayers(gpkg_path)
+print(layers)
+
+# Select the desired layers
+desired_layers = ["gemaal", "stuw", "sluis"]
+
+# Read the selected layers into GeoDataFrames
+uitlaten_inlaten_gdf = {}
+with fiona.open(gpkg_path, "r") as gpkg:
+    for layer_name in desired_layers:
+        if layer_name in layers:
+            gdf = gpd.read_file(gpkg_path, layer=layer_name)
+            uitlaten_inlaten_gdf[layer_name] = gdf
+
+# Spatial operations
+selected_waterschap_gdf = selected_waterschap_gdf.to_crs(
+    uitlaten_inlaten_gdf[desired_layers[0]].crs
+)
+
+# Perform the spatial join for the first layer
+points_within_waterschap_gdf = gpd.sjoin(
+    uitlaten_inlaten_gdf[desired_layers[0]], selected_waterschap_gdf, op="within"
+)
+
+# Combine the results from different layers (if needed)
+dfs_to_concat = []
+for layer_name in desired_layers[1:]:
+    result = gpd.sjoin(
+        uitlaten_inlaten_gdf[layer_name], selected_waterschap_gdf, op="within"
+    )
+    dfs_to_concat.append(result)
+
+# Concatenate the DataFrames
+points_within_waterschap_gdf = pd.concat(
+    [points_within_waterschap_gdf] + dfs_to_concat, ignore_index=True
+)
+
+# Drop the 'OBJECTID' column if it exists
+if "OBJECTID" in points_within_waterschap_gdf:
+    points_within_waterschap_gdf = points_within_waterschap_gdf.drop(
+        columns=["OBJECTID"]
+    )
+
+# Output to a GeoPackage
+output_gpkg = MODEL_DIR / "inlaten_uitlaten_ZZL.gpkg"
+points_within_waterschap_gdf.to_file(output_gpkg, driver="GPKG")
+
+# %% waterlopen Zuiderzeeland
+
+# Handmatig bewerkte waterlopen
+gpkg_path = MODEL_DIR / "primaire_waterlopen.gpkg"
+
+
+# %%

notebooks/ijsselmeermodel/kunstwerken.py

@@ -0,0 +1,76 @@
+# %%
+import os
+from pathlib import Path
+
+import fiona
+import geopandas as gpd
+import pandas as pd
+
+DATA_DIR = Path(os.getenv("RIBASIM_NL_DATA_DIR"))
+MODEL_DIR = Path(os.getenv("RIBASIM_NL_MODEL_DIR")) / "ijsselmeer"
+
+# %% kunstwerken Zuiderzeeland
+
+
+# Inlezen waterschapsgrenzen
+shapefile = DATA_DIR / r"nederland/Waterschapsgrenzen.shp"
+waterschapsgrenzen_gdf = gpd.read_file(shapefile)
+
+# Filter the GeoDataFrame to select the feature with value 'Waterschap Zuiderzeeland'
+selected_waterschap_gdf = waterschapsgrenzen_gdf[
+    waterschapsgrenzen_gdf["waterschap"] == "Waterschap Zuiderzeeland"
+]
+
+output_gpkg = MODEL_DIR / "ZZL_grens.gpkg"
+selected_waterschap_gdf.to_file(output_gpkg, driver="GPKG")
+
+# Path to the GeoPackage file
+gpkg_path = DATA_DIR / r"uitlaten_inlaten.gpkg"
+
+# List available layers in the GeoPackage
+layers = fiona.listlayers(gpkg_path)
+print(layers)
+
+# Select the desired layers
+desired_layers = ["gemaal", "stuw", "sluis"]
+
+# Read the selected layers into GeoDataFrames
+uitlaten_inlaten_gdf = {}
+with fiona.open(gpkg_path, "r") as gpkg:
+    for layer_name in desired_layers:
+        if layer_name in layers:
+            gdf = gpd.read_file(gpkg_path, layer=layer_name)
+            uitlaten_inlaten_gdf[layer_name] = gdf
+
+# Spatial operations
+selected_waterschap_gdf = selected_waterschap_gdf.to_crs(
+    uitlaten_inlaten_gdf[desired_layers[0]].crs
+)
+
+# Perform the spatial join for the first layer
+points_within_waterschap_gdf = gpd.sjoin(
+    uitlaten_inlaten_gdf[desired_layers[0]], selected_waterschap_gdf, op="within"
+)
+
+# Combine the results from different layers (if needed)
+dfs_to_concat = []
+for layer_name in desired_layers[1:]:
+    result = gpd.sjoin(
+        uitlaten_inlaten_gdf[layer_name], selected_waterschap_gdf, op="within"
+    )
+    dfs_to_concat.append(result)
+
+# Concatenate the DataFrames
+points_within_waterschap_gdf = pd.concat(
+    [points_within_waterschap_gdf] + dfs_to_concat, ignore_index=True
+)
+
+# Drop the 'OBJECTID' column if it exists
+if "OBJECTID" in points_within_waterschap_gdf:
+    points_within_waterschap_gdf = points_within_waterschap_gdf.drop(
+        columns=["OBJECTID"]
+    )
+
+# Output to a GeoPackage
+output_gpkg = MODEL_DIR / "inlaten_uitlaten_ZZL.gpkg"
+points_within_waterschap_gdf.to_file(output_gpkg, driver="GPKG")