Fix dommel basin areas (#133)

All code to fix: #132 

Some properties/methods in ribasim_nl.Model:
- `unassigned_basin_area`: basin area not assigned to a valid basin
node_id
- `basin_node_without_area`: basins without an area
- `fix_unassigned_basin_area()`: Assign a Basin node_id to a Basin /
Area if the Area doesn't contain a basin node_id

Results in model 2024.8.2: https://deltares.thegood.cloud/f/118021

---------

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

notebooks/de_dommel/modify_model.py

@@ -189,18 +189,16 @@
 for node_id in df.from_node_id:
     model.update_node(node_id, "Outlet", [outlet.Static(flow_rate=[100])])
 
-# for row in df.itertuples():
-#     area_select_df = model.basin.area.df[model.basin.area.df.geometry.contains(row.geometry)]
-#     if len(area_select_df) == 1:
-#         area_id = area_select_df.iloc[0]["node_id"]
-#         if row.node_id != area_id:
-#             print(f"{row.node_id} == {area_id}")
-#     elif area_select_df.empty:
-#         raise Exception(f"Basin {row.node_id} not within Area")
-#     else:
-#         raise Exception(f"Basin {row.node_id} contained by multiple areas")
-
-# %% write model
+# # see: https://github.com/Deltares/Ribasim-NL/issues/132
+model.basin.area.df.loc[model.basin.area.df.duplicated("node_id"), ["node_id"]] = -1
+model.basin.area.df.reset_index(drop=True, inplace=True)
+model.fix_unassigned_basin_area()
+model.fix_unassigned_basin_area(method="closest", distance=100)
+model.fix_unassigned_basin_area()
+
+model.basin.area.df = model.basin.area.df[~model.basin.area.df.node_id.isin(model.unassigned_basin_area.node_id)]
+
+# # %% write model
 ribasim_toml = ribasim_toml.parents[1].joinpath("DeDommel", ribasim_toml.name)
 model.write(ribasim_toml)
 

src/ribasim_nl/ribasim_nl/model.py

@@ -84,6 +84,16 @@ def find_node_id(self, ds_node_id=None, us_node_id=None, **kwargs) -> int:
         else:
             return node_ids[0]
 
+    @property
+    def unassigned_basin_area(self):
+        """Get unassigned basin area"""
+        return self.basin.area.df[~self.basin.area.df.node_id.isin(self.basin.node.df.node_id)]
+
+    @property
+    def basin_node_without_area(self):
+        """Get basin node without area"""
+        return self.basin.node.df[~self.basin.node.df.node_id.isin(self.basin.area.df.node_id)]
+
     def get_node_type(self, node_id: int):
         return self.node_table().df.set_index("node_id").at[node_id, "node_type"]
 
@@ -262,3 +272,39 @@ def find_closest_basin(self, geometry: BaseGeometry, max_distance: float | None)
                 )
 
         return self.basin[basin_node_id]
+
+    def fix_unassigned_basin_area(self, method: str = "within", distance: float = 100):
+        """Assign a Basin node_id to a Basin / Area if the Area doesn't contain a basin node_id.
+
+        Args:
+            method (str): method to find basin node_id; `within` or `closest`. First start with `within`. Default is `within`
+            distance (float, optional): for method closest, the distance to find an unassigned basin node_id. Defaults to 100.
+        """
+        if self.basin.node.df is not None:
+            if self.basin.area.df is not None:
+                basin_area_df = self.basin.area.df[~self.basin.area.df.node_id.isin(self.basin.node.df.node_id)]
+
+                for row in basin_area_df.itertuples():
+                    if method == "within":
+                        # check if area contains basin-nodes
+                        basin_df = self.basin.node.df[self.basin.node.df.within(row.geometry)]
+
+                    elif method == "closest":
+                        basin_df = self.basin.node.df[self.basin.node.df.within(row.geometry)]
+                        # if method is `distance` and basin_df is emtpy we create a new basin_df
+                        if basin_df.empty:
+                            basin_df = self.basin.node.df[self.basin.node.df.distance(row.geometry) < distance]
+
+                    else:
+                        ValueError(f"Supported methods are 'within' or 'closest', got '{method}'.")
+
+                    # check if basin_nodes within area are not yet assigned an area
+                    basin_df = basin_df[~basin_df.node_id.isin(self.basin.area.df.node_id)]
+
+                    # if we have one node left we are done
+                    if len(basin_df) == 1:
+                        self.basin.area.df.loc[row.Index, ["node_id"]] = basin_df.iloc[0].node_id
+            else:
+                raise ValueError("Assign Basin Area to your model first")
+        else:
+            raise ValueError("Assign a Basin Node to your model first")