Parameterize nzv

notebooks/noorderzijlvest/01_fix_model.py

@@ -9,6 +9,8 @@
 from shapely.ops import snap, split
 
 from ribasim_nl import CloudStorage, Model, Network, NetworkValidator
+from ribasim_nl.case_conversions import pascal_to_snake_case
+from ribasim_nl.gkw import get_data_from_gkw
 from ribasim_nl.reset_static_tables import reset_static_tables
 
 cloud = CloudStorage()
@@ -290,8 +292,44 @@ def get_network_node(point):
 
 model.remove_unassigned_basin_area()
 
+
+# %% TabulatedRatingCurve to Outlet
+
+# TabulatedRatingCurve to Outlet
+for row in model.node_table().df[model.node_table().df.node_type == "TabulatedRatingCurve"].itertuples():
+    node_id = row.Index
+    model.update_node(node_id=node_id, node_type="Outlet")
+
+# %% sanitize node-tables
+node_columns = model.basin.node.columns() + ["meta_code_waterbeheerder", "meta_categorie"]
+
+# name to code
+model.outlet.node.df.loc[:, "meta_code_waterbeheerder"] = model.outlet.node.df.name
+model.pump.node.df.loc[:, "meta_code_waterbeheerder"] = model.pump.node.df.name
+
+df = get_data_from_gkw(authority="Noorderzijlvest", layers=["gemaal", "stuw", "sluis"])
+df.set_index("code", inplace=True)
+names = df["naam"]
+names.loc["KSL011"] = "R.J. Cleveringensluizen"
+
+# remove names and clean columns
+for node_type in model.node_table().df.node_type.unique():
+    table = getattr(model, pascal_to_snake_case(node_type))
+
+    if "meta_code_waterbeheerder" in table.node.df.columns:
+        table.node.df.loc[:, "name"] = table.node.df["meta_code_waterbeheerder"].apply(
+            lambda x: names[x] if x in names.index.to_numpy() else ""
+        )
+    else:
+        table.node.df.name = ""
+    columns = [col for col in table.node.df.columns if col in node_columns]
+    table.node.df = table.node.df[columns]
+
+
 #  %% write model
 model.use_validation = True
 model.write(ribasim_toml)
 model.report_basin_area()
 model.report_internal_basins()
+
+# %%

notebooks/noorderzijlvest/02_parameterize_model.py

@@ -0,0 +1,39 @@
+# %%
+import time
+
+from ribasim_nl import CloudStorage, Model
+
+cloud = CloudStorage()
+authority = "Noorderzijlvest"
+short_name = "nzv"
+
+
+static_data_xlsx = cloud.joinpath(
+    "Noorderzijlvest",
+    "verwerkt",
+    "parameters",
+    "static_data.xlsx",
+)
+ribasim_dir = cloud.joinpath(authority, "modellen", f"{authority}_fix_model")
+ribasim_toml = ribasim_dir / f"{short_name}.toml"
+
+# you need the excel, but the model should be local-only by running 01_fix_model.py
+cloud.synchronize(filepaths=[static_data_xlsx])
+cloud.synchronize(filepaths=[ribasim_dir], check_on_remote=False)
+
+# %%
+
+# read
+model = Model.read(ribasim_toml)
+
+start_time = time.time()
+# %%
+# parameterize
+model.parameterize(static_data_xlsx=static_data_xlsx, precipitation_mm_per_day=10)
+print("Elapsed Time:", time.time() - start_time, "seconds")
+
+# %%
+
+# Write model
+ribasim_toml = cloud.joinpath(authority, "modellen", f"{authority}_parameterized_model", f"{short_name}.toml")
+model.write(ribasim_toml)

src/ribasim_nl/ribasim_nl/downstream.py

@@ -0,0 +1,43 @@
+from collections import deque
+
+import networkx as nx
+
+
+def downstream_nodes(graph: nx.DiGraph, node_id: int, stop_at_outlet: bool = False):
+    """Efficiently find all downstream nodes in a directed graph starting from a given node,
+    stopping traversal at nodes stopping at the next outlet.
+
+    Parameters
+    ----------
+    - graph (nx.DiGraph): The directed graph.
+    - start_node: The node to start the search from.
+    - stop_at_outlet (bool): To stop at the next inlet(s)
+
+    Returns
+    -------
+    - set: A set of all downstream nodes excluding the starting node.
+    """  # noqa: D205
+    visited = set()  # To keep track of visited nodes
+    node_ids = set({node_id})  # To store the result
+
+    # BFS using a deque
+    queue = deque([node_id])
+
+    while queue:
+        current_node = queue.popleft()
+
+        # Avoid re-visiting nodes
+        if current_node in visited:
+            continue
+        visited.add(current_node)
+
+        # Check successors (downstream neighbors)
+        for successor in graph.successors(current_node):
+            if successor not in visited:
+                node_ids.add(successor)
+
+                # Stop traversal if 'function' is 'outlet'
+                if (not stop_at_outlet) | (not graph.nodes[successor].get("function") == "outlet"):
+                    queue.append(successor)
+
+    return node_ids

src/ribasim_nl/ribasim_nl/model.py

@@ -16,8 +16,11 @@
 from shapely.geometry.base import BaseGeometry
 
 from ribasim_nl.case_conversions import pascal_to_snake_case
+from ribasim_nl.downstream import downstream_nodes
 from ribasim_nl.geometry import split_basin
+from ribasim_nl.parametrization.parameterize import Parameterize
 from ribasim_nl.run_model import run
+from ribasim_nl.upstream import upstream_nodes
 
 manning_data = manning_resistance.Static(length=[100], manning_n=[0.04], profile_width=[10], profile_slope=[1])
 level_data = level_boundary.Static(level=[0])
@@ -75,6 +78,14 @@ class Model(Model):
     _basin_results: Results | None = None
     _basin_outstate: Results | None = None
     _graph: nx.Graph | None = None
+    _parameterize: Parameterize | None = None
+
+    def __init__(self, **data):
+        super().__init__(**data)
+        self._parameterize = Parameterize(model=self)
+
+    def parameterize(self, **kwargs):
+        self._parameterize.run(**kwargs)
 
     @property
     def basin_results(self):
@@ -94,14 +105,31 @@ def basin_outstate(self):
     def graph(self):
         # create a DiGraph from edge-table
         if self._graph is None:
-            self._graph = nx.from_pandas_edgelist(
+            graph = nx.from_pandas_edgelist(
                 df=self.edge.df[["from_node_id", "to_node_id"]],
                 source="from_node_id",
                 target="to_node_id",
                 create_using=nx.DiGraph,
             )
+            if "meta_function" not in self.node_table().df.columns:
+                node_attributes = {node_id: {"function": ""} for node_id in self.node_table().df.index}
+            else:
+                node_attributes = (
+                    self.node_table()
+                    .df.rename(columns={"meta_function": "function"})[["function"]]
+                    .to_dict(orient="index")
+                )
+            nx.set_node_attributes(graph, node_attributes)
+
+            self._graph = graph
+
         return self._graph
 
+    @property
+    def reset_graph(self):
+        self._graph = None
+        return self.graph
+
     @property
     def next_node_id(self):
         return self.node_table().df.index.max() + 1
@@ -126,28 +154,35 @@ def update_state(self, time_stamp: pd.Timestamp | None = None):
         self.basin.state.df = df
 
     # methods relying on networkx. Discuss making this all in a subclass of Model
-    def _upstream_nodes(self, node_id):
+    def _upstream_nodes(self, node_id, **kwargs):
         # get upstream nodes
-        return list(nx.traversal.bfs_tree(self.graph, node_id, reverse=True))
+        #     return list(nx.traversal.bfs_tree(self.graph, node_id, reverse=True))
+        return upstream_nodes(graph=self.graph, node_id=node_id, **kwargs)
 
-    def _downstream_nodes(self, node_id):
+    def _downstream_nodes(self, node_id, **kwargs):
         # get downstream nodes
-        return list(nx.traversal.bfs_tree(self.graph, node_id))
+        return downstream_nodes(graph=self.graph, node_id=node_id, **kwargs)
+        # return list(nx.traversal.bfs_tree(self.graph, node_id))
 
-    def get_upstream_basins(self, node_id):
+    def get_upstream_basins(self, node_id, **kwargs):
         # get upstream basin area
-        upstream_node_ids = self._upstream_nodes(node_id)
+        upstream_node_ids = self._upstream_nodes(node_id, **kwargs)
         return self.basin.area.df[self.basin.area.df.node_id.isin(upstream_node_ids)]
 
-    def get_upstream_edges(self, node_id):
+    def get_downstream_basins(self, node_id, **kwargs):
+        # get upstream basin area
+        downstream_node_ids = self._downstream_nodes(node_id, **kwargs)
+        return self.basin.area.df[self.basin.area.df.node_id.isin(downstream_node_ids)]
+
+    def get_upstream_edges(self, node_id, **kwargs):
         # get upstream edges
-        upstream_node_ids = self._upstream_nodes(node_id)
+        upstream_node_ids = self._upstream_nodes(node_id, **kwargs)
         mask = self.edge.df.from_node_id.isin(upstream_node_ids[1:]) & self.edge.df.to_node_id.isin(upstream_node_ids)
         return self.edge.df[mask]
 
-    def get_downstream_edges(self, node_id):
+    def get_downstream_edges(self, node_id, **kwargs):
         # get upstream edges
-        downstream_node_ids = self._downstream_nodes(node_id)
+        downstream_node_ids = self._downstream_nodes(node_id, **kwargs)
         mask = self.edge.df.from_node_id.isin(downstream_node_ids) & self.edge.df.to_node_id.isin(
             downstream_node_ids[1:]
         )
@@ -501,7 +536,7 @@ def update_basin_area(self, node_id: int, geometry: Polygon | MultiPolygon, basi
 
         self.basin.area.df.loc[mask, ["node_id"]] = node_id
 
-    def add_basin_area(self, geometry: MultiPolygon, node_id: int | None = None):
+    def add_basin_area(self, geometry: MultiPolygon, node_id: int | None = None, meta_streefpeil: float | None = None):
         # if node_id is None, get an available node_id
         if pd.isna(node_id):
             basin_df = self.basin.node.df[self.basin.node.df.within(geometry)]
@@ -524,7 +559,10 @@ def add_basin_area(self, geometry: MultiPolygon, node_id: int | None = None):
                 raise ValueError(f"geometry-type {geometry.geom_type} is not valid. Provide (Multi)Polygon instead")
 
         # if all correct, assign
-        area_df = gpd.GeoDataFrame({"node_id": [node_id], "geometry": [geometry]}, crs=self.crs)
+        data = {"node_id": [node_id], "geometry": [geometry]}
+        if meta_streefpeil is not None:
+            data = {**data, "meta_streefpeil": [meta_streefpeil]}
+        area_df = gpd.GeoDataFrame(data, crs=self.crs)
         area_df.index.name = "fid"
         area_df.index += self.basin.area.df.index.max() + 1
         self.basin.area.df = pd.concat([self.basin.area.df, area_df])

src/ribasim_nl/ribasim_nl/parametrization/basin_tables.py

@@ -0,0 +1,76 @@
+import pandas as pd
+
+from ribasim_nl.model import Model
+from ribasim_nl.parametrization.conversions import round_to_precision
+from ribasim_nl.parametrization.empty_table import empty_table_df
+
+
+def update_basin_static(
+    model: Model, precipitation_mm_per_day: float | None = None, evaporation_mm_per_day: float | None = None
+):
+    """Add precipitation and/or evaporation to the model.basin.static table from basin.area
+
+    Args:
+        model (Model): Ribasim Model
+        precipitation_mm_per_day (float | None, optional): Precipitation in mm/day. Defaults to None.
+        evaporation_mm_per_day (float | None, optional): Evaporation in mm/day. Defaults to None.
+    """
+    static_df = empty_table_df(model, node_type="Basin", table_type="Static", fill_value=0)
+
+    area = model.basin.area.df.set_index("node_id").geometry.area
+    if precipitation_mm_per_day is not None:
+        precipitation = area * (precipitation_mm_per_day * 0.001 / 86400)  # m3/s
+        static_df.loc[:, "precipitation"] = precipitation[static_df.node_id].to_numpy()
+    if evaporation_mm_per_day is not None:
+        evaporation = area * (evaporation_mm_per_day * 0.001 / 86400)  # m3/s
+        static_df.loc[:, "evaporation"] = evaporation[static_df.node_id].to_numpy()
+
+
+def update_basin_profile(
+    model: Model,
+    percentages_map: dict = {"hoofdwater": 90, "doorgaand": 10, "bergend": 3},
+    default_percentage: int = 10,
+    profile_depth=2,
+):
+    # read profile from basin-table
+    profile = model.basin.area.df.copy()
+
+    # determine the profile area, which is also used for the profile
+    # profile["area"] = profile["geometry"].area * (default_percentage / 100)
+    profile = profile[["node_id", "meta_streefpeil", "geometry"]]
+    profile = profile.rename(columns={"meta_streefpeil": "level"})
+
+    # get open-water percentages per category
+    profile["percentage"] = default_percentage
+    for category, percentage in percentages_map.items():
+        node_ids = model.basin.node.df.loc[model.basin.node.df.meta_categorie == category].index.to_numpy()
+        profile.loc[profile.node_id.isin(node_ids), "percentage"] = percentage
+
+    # calculate area at invert from percentage
+    profile["area"] = profile["geometry"].area * profile["percentage"] / 100
+    profile.drop(columns=["geometry", "percentage"], inplace=True)
+
+    # define the profile-bottom
+    profile_bottom = profile.copy()
+    profile_bottom["area"] = 0.1
+    profile_bottom["level"] -= profile_depth
+
+    # define the profile slightly above the bottom of the bakje
+    profile_slightly_above_bottom = profile.copy()
+    profile_slightly_above_bottom["level"] -= profile_depth - 0.01  # remain one centimeter above the bottom
+
+    # combine all profiles by concatenating them, and sort on node_id, level and area.
+    profile_df = pd.concat([profile_bottom, profile_slightly_above_bottom, profile])
+    profile_df = profile_df.sort_values(by=["node_id", "level", "area"], ascending=True).reset_index(drop=True)
+    profile_df.loc[:, "area"] = profile_df["area"].apply(round_to_precision, args=(0.1,))
+
+    model.basin.profile.df = profile_df
+
+
+def update_basin_state(model: Model):
+    """Update basin state by max profile-level
+
+    Args:
+        model (Model): Ribasim Model
+    """
+    model.basin.state.df = model.basin.profile.df.groupby("node_id").max().reset_index()[["node_id", "level"]]