IJsselmeer model (#26)

Notebooks and updates for generating the ijsselmeermodel: - included `layer` in generate_model_id to ensure uniqueness - output type annotaties in code_utils as discussed: #25 (comment) linked issues: #3: organize data #28: build model Discussion: https://github.com/Deltares/Ribasim-NL/discussions/27 --------- Co-authored-by: ngoorden <[email protected]> Co-authored-by: Martijn Visser <[email protected]>
Deltares · Nov 7, 2023 · ff0e5fe · ff0e5fe
1 parent b343c71
commit ff0e5fe
Show file tree

Hide file tree

Showing 26 changed files with 10,341 additions and 7,847 deletions.
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -17,12 +17,10 @@ repos:
     rev: v0.1.3
     hooks:
       - id: ruff
+        types_or: [python, pyi, jupyter]
         args: [--fix, --exit-non-zero-on-fix]
-  - repo: https://github.com/psf/black
-    rev: 23.10.1
-    hooks:
-      - id: black
-      - id: black-jupyter
+      - id: ruff-format
+        types_or: [python, pyi, jupyter]
   - repo: https://github.com/kynan/nbstripout
     rev: 0.6.1
     hooks:

diff --git a/.vscode/extensions.json b/.vscode/extensions.json
@@ -1,7 +1,6 @@
 {
     "recommendations": [
         "ms-python.python",
-        "ms-python.black-formatter",
         "ms-python.mypy-type-checker",
         "charliermarsh.ruff",
         "njpwerner.autodocstring"

diff --git a/.vscode/settings_template.json b/.vscode/settings_template.json
@@ -1,7 +1,7 @@
 {
     "notebook.formatOnSave.enabled": true,
     "[python]": {
-        "editor.defaultFormatter": "ms-python.black-formatter",
+        "editor.defaultFormatter": "charliermarsh.ruff",
         "editor.formatOnSave": true,
         "editor.codeActionsOnSave": {
             "source.fixAll": true

diff --git a/environment.yml b/environment.yml
@@ -14,7 +14,7 @@ dependencies:
   - matplotlib
   - mypy
   - openpyxl
-  - pandas!=2.1.0
+  - pandas
   - pandas-stubs
   - pandera
   - pip
@@ -27,6 +27,7 @@ dependencies:
   - pytest
   - pytest-cov
   - python>=3.9
+  - rasterstats
   - ruff
   - shapely>=2.0
   - tomli

diff --git a/notebooks/ijsselmeermodel/basins.py b/notebooks/ijsselmeermodel/basins.py
@@ -0,0 +1,142 @@
+# %%
+import os
+from pathlib import Path
+
+import geopandas as gpd
+import pandas as pd
+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
+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")
+
+# rws_krw_gdf.loc[krw_ids].explore()
+
+krw_cutlines_gdf = gpd.read_file(MODEL_DIR / "model_data.gpkg", layer="krw_cutlines")
+
+
+def strip_code(code):
+    return code.split("_", 1)[-1]
+
+
+def user_id(code, wbh_code, code_postfix=None):
+    code = strip_code(code)
+    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,
+                        rijkswater=strip_code(code),
+                    )
+                    basin_polygon = basin_multi_polygon.geoms[1]
+                add_basin(
+                    user_id=user_id(code, "80", cut_line.cut_order + 1),
+                    geometry=basin_polygon,
+                    rijkswater=strip_code(code),
+                )
+            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,
+            rijkswater=strip_code(code),
+        )
+
+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"),
+        peilvak=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.to_file(MODEL_DIR / "model_data.gpkg", layer="basin_area")
+basins_gdf.loc[:, "geometry"] = basins_to_points(basins_gdf["geometry"])
+basins_gdf.to_file(MODEL_DIR / "model_data.gpkg", layer="basin")
+
+## %% 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"]
diff --git a/notebooks/ijsselmeermodel/brongegevens.ipynb b/notebooks/ijsselmeermodel/brongegevens.ipynb
@@ -0,0 +1,44 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import geopandas as gpd"
+   ]
+  },
+  {
+   "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
+}
diff --git a/notebooks/ijsselmeermodel/ijsselmeer_model.py b/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"
+
+
+# %%
diff --git a/notebooks/ijsselmeermodel/kunstwerken.py b/notebooks/ijsselmeermodel/kunstwerken.py
@@ -0,0 +1,75 @@
+# %%
+import os
+from pathlib import Path
+
+import geopandas as gpd
+import pandas as pd
+from hydamo import code_utils
+
+DATA_DIR = Path(os.getenv("RIBASIM_NL_DATA_DIR"))
+MODEL_DIR = Path(os.getenv("RIBASIM_NL_MODEL_DIR")) / "ijsselmeer"
+EXCEL_FILE = "uitlaten_inlaten.xlsx"
+BGT_CODES = ["W0650", "P0024"]
+
+KUNSTWERKEN_XLSX = Path(DATA_DIR) / EXCEL_FILE
+MODEL_DATA_GPKG = Path(MODEL_DIR) / "model_data.gpkg"
+
+basin_gdf = gpd.read_file(MODEL_DATA_GPKG, layer="basin")
+
+kunstwerken_df = pd.read_excel(KUNSTWERKEN_XLSX)
+kunstwerken_df = kunstwerken_df.loc[kunstwerken_df.bgt_code.isin(BGT_CODES)]
+kunstwerken_gdf = gpd.GeoDataFrame(
+    kunstwerken_df,
+    geometry=gpd.points_from_xy(x=kunstwerken_df.x, y=kunstwerken_df.y),
+    crs=28992,
+)
+
+# %%
+pump_gdf = kunstwerken_gdf[kunstwerken_gdf.dm_type == "uitlaat"][
+    ["dm_capaciteit", "user_id", "peilvak", "rijkswater", "geometry"]
+].copy()
+pump_gdf.rename(
+    columns={"dm_capaciteit": "flow_rate", "peilvak": "id_from", "rijkswater": "id_to"},
+    inplace=True,
+)
+pump_gdf["id_from"] = pump_gdf["id_from"].apply(
+    lambda x: code_utils.generate_model_id(code=x, layer="basin", bgt_code="W0650")
+)
+pump_gdf[pump_gdf.flow_rate.isna()]["flow_rate"] = 0
+pump_gdf.to_file(MODEL_DIR / "model_data.gpkg", layer="pump")
+
+# %%
+outlet_gdf = (
+    kunstwerken_gdf[kunstwerken_gdf.dm_type == "inlaat"][
+        ["dm_capaciteit", "user_id", "peilvak", "rijkswater", "geometry"]
+    ]
+    .copy()
+    .rename(columns={"dm_capaciteit": "flow_rate"})
+)
+outlet_gdf.rename(
+    columns={"dm_capaciteit": "flow_rate", "peilvak": "id_to", "rijkswater": "id_from"},
+    inplace=True,
+)
+outlet_gdf["id_to"] = outlet_gdf["id_to"].apply(
+    lambda x: code_utils.generate_model_id(code=x, layer="basin", bgt_code="W0650")
+)
+outlet_gdf[outlet_gdf.flow_rate.isna()]["flow_rate"] = 0
+outlet_gdf.to_file(MODEL_DIR / "model_data.gpkg", layer="outlet")
+
+# %%
+
+resistance_gdf = gpd.read_file(
+    DATA_DIR.joinpath("ijsselmeergebied", "hydamo.gpkg"), layer="sluis"
+)
+
+resistance_gdf.rename(
+    columns={"rijkswater_naar": "id_to", "rijkswater_van": "id_from"}, inplace=True
+)
+
+resistance_gdf["user_id"] = resistance_gdf["code"].apply(
+    lambda x: code_utils.generate_model_id(code=x, layer="sluis", bgt_code="L0002")
+)
+resistance_gdf["resistance"] = 1000
+
+resistance_gdf.to_file(MODEL_DIR / "model_data.gpkg", layer="resistance")
+# %%