Add optional 'Basin / area' spatial table (#1200)

Fixes #960
It looks like this when loaded in QGIS:


![image](https://github.com/Deltares/Ribasim/assets/4471859/d0c4ad97-9a91-483d-9568-2c3dbb8b3d5b)

The color is still random as I didn't implement a `renderer` yet, though
that can be a follow up.
I did have to take special care to make sure these polygons were not
plotted on top, hiding the model.

MultiPolygons and missing polygons also work.

docs/core/usage.qmd

@@ -310,6 +310,16 @@ Internally this get converted to two functions, $A(S)$ and $h(S)$, by integratin
 The minimum area cannot be zero to avoid numerical issues.
 The maximum area is used to convert the precipitation flux into an inflow.
 
+## Basin / area
+
+The optional area table is not used during computation, but provides a place to associate areas in the form of polygons to Basins.
+Using this makes it easier to recognize which water or land surfaces are represented by Basins.
+
+column    | type                    | restriction
+--------- | ----------------------- | -----------
+node_id   | Int                     | sorted
+geom      | Polygon or MultiPolygon | (optional)
+
 ## Basin / subgrid
 
 The subgrid_level table defines a piecewise linear interpolation from a basin water level to a subgrid element water level.

python/ribasim/ribasim/config.py

@@ -2,7 +2,8 @@
 
 from pydantic import Field
 
-from ribasim.input_base import ChildModel, NodeModel, TableModel
+from ribasim.geometry import BasinAreaSchema
+from ribasim.input_base import ChildModel, NodeModel, SpatialTableModel, TableModel
 
 # These schemas are autogenerated
 from ribasim.schemas import (
@@ -174,6 +175,10 @@ class Basin(NodeModel):
         default_factory=TableModel[BasinSubgridSchema],
         json_schema_extra={"sort_keys": ["subgrid_id", "basin_level"]},
     )
+    area: SpatialTableModel[BasinAreaSchema] = Field(
+        default_factory=SpatialTableModel[BasinAreaSchema],
+        json_schema_extra={"sort_keys": ["node_id"]},
+    )
 
 
 class ManningResistance(NodeModel):

python/ribasim/ribasim/geometry/__init__.py

@@ -1,4 +1,5 @@
+from ribasim.geometry.area import BasinAreaSchema
 from ribasim.geometry.edge import Edge
 from ribasim.geometry.node import Node
 
-__all__ = ["Edge", "Node"]
+__all__ = ["BasinAreaSchema", "Edge", "Node"]

python/ribasim/ribasim/geometry/area.py

@@ -0,0 +1,12 @@
+from typing import Any
+
+import pandera as pa
+from pandera.typing import Series
+from pandera.typing.geopandas import GeoSeries
+
+from ribasim.schemas import _BaseSchema
+
+
+class BasinAreaSchema(_BaseSchema):
+    node_id: Series[int]
+    geometry: GeoSeries[Any] = pa.Field(default=None, nullable=True)

python/ribasim/tests/test_model.py

@@ -129,6 +129,9 @@ def test_node_ids_unsequential(basic):
 
 def test_tabulated_rating_curve_model(tabulated_rating_curve, tmp_path):
     model_orig = tabulated_rating_curve
+    basin_area = tabulated_rating_curve.basin.area.df
+    assert basin_area is not None
+    assert basin_area.geometry.geom_type[0] == "Polygon"
     model_orig.write(tmp_path / "tabulated_rating_curve/ribasim.toml")
     model_new = Model.read(tmp_path / "tabulated_rating_curve/ribasim.toml")
     pd.testing.assert_series_equal(

python/ribasim_testmodels/ribasim_testmodels/basic.py

@@ -287,34 +287,9 @@ def tabulated_rating_curve_model() -> ribasim.Model:
     Only the upstream Basin receives a (constant) precipitation.
     """
 
-    # Setup the basins:
-    profile = pd.DataFrame(
-        data={
-            "node_id": [1, 1, 4, 4],
-            "area": [0.01, 1000.0] * 2,
-            "level": [0.0, 1.0] * 2,
-        }
-    )
-
-    # Convert steady forcing to m/s
-    # 2 mm/d precipitation
-    seconds_in_day = 24 * 3600
-    precipitation = 0.002 / seconds_in_day
-    # only the upstream basin gets precipitation
-    static = pd.DataFrame(
-        data={
-            "node_id": [1, 4],
-            "precipitation": [precipitation, 0.0],
-        }
-    )
-    state = pd.DataFrame(
-        data={"node_id": static["node_id"], "level": 0.04471158417652035}
-    )
-
-    basin = ribasim.Basin(profile=profile, static=static, state=state)
-
     # Set up a rating curve node:
     # Discharge: lose 1% of storage volume per day at storage = 1000.0.
+    seconds_in_day = 24 * 3600
     q1000 = 1000.0 * 0.01 / seconds_in_day
 
     rating_curve = ribasim.TabulatedRatingCurve(
@@ -354,18 +329,60 @@ def tabulated_rating_curve_model() -> ribasim.Model:
     )
     node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])
 
-    node_id, node_type = ribasim.Node.node_ids_and_types(basin, rating_curve)
-
     # Make sure the feature id starts at 1: explicitly give an index.
     node = ribasim.Node(
         df=gpd.GeoDataFrame(
-            data={"node_type": node_type},
-            index=pd.Index(node_id, name="fid"),
+            data={
+                "node_type": [
+                    "Basin",
+                    "TabulatedRatingCurve",
+                    "TabulatedRatingCurve",
+                    "Basin",
+                ]
+            },
+            index=pd.Index([1, 2, 3, 4], name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
         )
     )
 
+    # Setup the basins:
+    profile = pd.DataFrame(
+        data={
+            "node_id": [1, 1, 4, 4],
+            "area": [0.01, 1000.0] * 2,
+            "level": [0.0, 1.0] * 2,
+        }
+    )
+
+    # Convert steady forcing to m/s
+    # 2 mm/d precipitation
+    precipitation = 0.002 / seconds_in_day
+    # only the upstream basin gets precipitation
+    static = pd.DataFrame(
+        data={
+            "node_id": [1, 4],
+            "precipitation": [precipitation, 0.0],
+        }
+    )
+    state = pd.DataFrame(
+        data={"node_id": static["node_id"], "level": 0.04471158417652035}
+    )
+    # Turn the basin node point geometries into polygons describing the area
+    assert node.df is not None
+    basin_area = (
+        node.df[node.df.node_type == "Basin"]
+        .geometry.buffer(1.0)
+        .reset_index(drop=True)
+    )
+    area = gpd.GeoDataFrame(
+        data={"node_id": static.node_id},
+        geometry=basin_area,
+        crs="EPSG:28992",
+    )
+
+    basin = ribasim.Basin(profile=profile, static=static, state=state, area=area)
+
     # Setup the edges:
     from_id = np.array([1, 1, 2, 3], dtype=np.int64)
     to_id = np.array([2, 3, 4, 4], dtype=np.int64)

ribasim_qgis/core/nodes.py

@@ -4,7 +4,7 @@
 The classes specify:
 
 * The (unabbreviated) name
-* The type of geometry (No geometry, point, linestring, polygon)
+* The type of geometry ("No Geometry", "Point", "LineString", "Polygon")
 * The required attributes of the attribute table
 
 Each node layer is (optionally) represented in multiple places:
@@ -263,7 +263,7 @@ def attributes(cls) -> list[QgsField]:
 
     @classmethod
     def geometry_type(cls) -> str:
-        return "Linestring"
+        return "LineString"
 
     @classmethod
     def input_type(cls) -> str:
@@ -366,7 +366,7 @@ def input_type(cls) -> str:
 
     @classmethod
     def geometry_type(cls) -> str:
-        return "No geometry"
+        return "No Geometry"
 
     @classmethod
     def attributes(cls) -> list[QgsField]:
@@ -421,6 +421,24 @@ def attributes(cls) -> list[QgsField]:
         ]
 
 
+class BasinArea(Input):
+    @classmethod
+    def input_type(cls) -> str:
+        return "Basin / area"
+
+    @classmethod
+    def is_spatial(cls):
+        return True
+
+    @classmethod
+    def geometry_type(cls) -> str:
+        return "Polygon"
+
+    @classmethod
+    def attributes(cls) -> list[QgsField]:
+        return [QgsField("node_id", QVariant.Int)]
+
+
 class BasinState(Input):
     @classmethod
     def input_type(cls) -> str:
@@ -692,7 +710,7 @@ def input_type(cls) -> str:
 
     @classmethod
     def geometry_type(cls) -> str:
-        return "LineString"
+        return "No Geometry"
 
     @classmethod
     def attributes(cls) -> list[QgsField]:
@@ -710,7 +728,7 @@ def input_type(cls) -> str:
 
     @classmethod
     def geometry_type(cls) -> str:
-        return "LineString"
+        return "No Geometry"
 
     @classmethod
     def attributes(cls) -> list[QgsField]:
@@ -732,7 +750,7 @@ def input_type(cls) -> str:
 
     @classmethod
     def geometry_type(cls) -> str:
-        return "LineString"
+        return "No Geometry"
 
     @classmethod
     def attributes(cls) -> list[QgsField]:

ribasim_qgis/widgets/dataset_widget.py

@@ -62,16 +62,17 @@ def items(self) -> list[QTreeWidgetItem]:
         root = self.invisibleRootItem()
         return [root.child(i) for i in range(root.childCount())]
 
-    def add_item(self, name: str):
+    def add_item(self, name: str) -> QTreeWidgetItem:
         item = QTreeWidgetItem()
         self.addTopLevelItem(item)
         item.setText(0, name)
         return item
 
-    def add_node_layer(self, element: Input) -> None:
+    def add_node_layer(self, element: Input) -> QTreeWidgetItem:
         # These are mandatory elements, cannot be unticked
         item = self.add_item(name=element.input_type())
-        item.element = element
+        item.element = element  # type: ignore[attr-defined]
+        return item
 
     def remove_geopackage_layers(self) -> None:
         """
@@ -260,16 +261,32 @@ def load_geopackage(self) -> None:
         self.dataset_tree.clear()
         geo_path = get_database_path_from_model_file(self.path)
         nodes = load_nodes_from_geopackage(geo_path)
-        for node_layer in nodes.values():
-            self.dataset_tree.add_node_layer(node_layer)
+
         name = self.path.stem
         self.ribasim_widget.create_groups(name)
-        for item in self.dataset_tree.items():
+
+        # Make sure "Node", "Edge", "Basin / area" are the top three layers
+        node = nodes.pop("Node")
+        item = self.dataset_tree.add_node_layer(node)
+        self.add_item_to_qgis(item)
+
+        edge = nodes.pop("Edge")
+        item = self.dataset_tree.add_node_layer(edge)
+        self.add_item_to_qgis(item)
+
+        basin_area_layer = nodes.pop("Basin / area", None)
+        if basin_area_layer is not None:
+            item = self.dataset_tree.add_node_layer(basin_area_layer)
+            self.add_item_to_qgis(item)
+
+        # Add the remaining layers
+        for node_layer in nodes.values():
+            item = self.dataset_tree.add_node_layer(node_layer)
             self.add_item_to_qgis(item)
 
         # Connect node and edge layer to derive connectivities.
-        self.node_layer = nodes["Node"].layer
-        self.edge_layer = nodes["Edge"].layer
+        self.node_layer = node.layer
+        self.edge_layer = edge.layer
         self.edge_layer.editingStopped.connect(self.explode_and_connect)
         return
 

utils/templates/model.py.jinja

@@ -4,11 +4,13 @@ import pandera as pa
 from pandera.dtypes import Timestamp
 from pandera.typing import Series
 
+
 class _BaseSchema(pa.DataFrameModel):
     class Config:
         add_missing_columns = True
         coerce = True
 
+
 {% for m in models %}
 class {{m[:name]}}Schema(_BaseSchema):
 {% for f in m[:fields] %}