Support specifying coordinate reference system (CRS) of geometries

core/src/config.jl

@@ -120,6 +120,7 @@ end
 @option @addnodetypes struct Toml <: TableOption
     starttime::DateTime
     endtime::DateTime
+    crs::String
     ribasim_version::String
     input_dir::String
     results_dir::String

core/test/data/config_test.toml

@@ -1,5 +1,6 @@
 starttime = 2019-01-01
 endtime = 2019-12-31
+crs = "EPSG:28992"
 input_dir = "../../generated_testmodels/lhm"
 results_dir = "../../generated_testmodels/lhm"
 ribasim_version = "2024.6.1"

core/test/data/logging_test_loglevel_debug.toml

@@ -1,5 +1,6 @@
 starttime = 2019-01-01
 endtime = 2019-12-31
+crs = "EPSG:28992"
 input_dir = "."
 results_dir = "results"
 ribasim_version = "2024.6.1"

core/test/data/logging_test_no_loglevel.toml

@@ -1,5 +1,6 @@
 starttime = 2019-01-01
 endtime = 2019-12-31
+crs = "EPSG:28992"
 input_dir = "."
 results_dir = "results"
 ribasim_version = "2024.6.1"

core/test/docs.toml

@@ -2,6 +2,7 @@
 # can also be set to a date-time like 1979-05-27T07:32:00
 starttime = 2019-01-01 # required
 endtime = 2021-01-01   # required
+crs = "EPSG:4326"      # required
 
 # input files
 input_dir = "."         # required

core/test/io_test.jl

@@ -8,6 +8,7 @@
         database = "path/to/file",
         input_dir = ".",
         results_dir = "results",
+        crs = "EPSG:28992",
         ribasim_version = string(Ribasim.pkgversion(Ribasim)),
     )
     config = Ribasim.Config(toml, "model")
@@ -21,6 +22,7 @@
         database = "path/to/file",
         input_dir = "input",
         results_dir = "results",
+        crs = "EPSG:28992",
         ribasim_version = string(Ribasim.pkgversion(Ribasim)),
     )
     config = Ribasim.Config(toml, "model")
@@ -34,6 +36,7 @@
         database = "/path/to/file",
         input_dir = ".",
         results_dir = "results",
+        crs = "EPSG:28992",
         ribasim_version = string(Ribasim.pkgversion(Ribasim)),
     )
     config = Ribasim.Config(toml)

docs/python/examples.ipynb

@@ -76,7 +76,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\")"
+    "model = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:4326\")"
    ]
   },
   {
@@ -457,7 +457,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\")"
+    "model = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:4326\")"
    ]
   },
   {
@@ -771,10 +771,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = Model(\n",
-    "    starttime=\"2020-01-01\",\n",
-    "    endtime=\"2020-12-01\",\n",
-    ")"
+    "model = Model(starttime=\"2020-01-01\", endtime=\"2020-12-01\", crs=\"EPSG:4326\")"
    ]
   },
   {
@@ -1037,10 +1034,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = Model(\n",
-    "    starttime=\"2020-01-01\",\n",
-    "    endtime=\"2020-01-20\",\n",
-    ")"
+    "model = Model(starttime=\"2020-01-01\", endtime=\"2020-01-20\", crs=\"EPSG:4326\")"
    ]
   },
   {
@@ -1430,10 +1424,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model = Model(\n",
-    "    starttime=\"2020-01-01\",\n",
-    "    endtime=\"2020-02-01\",\n",
-    ")"
+    "model = Model(starttime=\"2020-01-01\", endtime=\"2020-02-01\", crs=\"EPSG:4326\")"
    ]
   },
   {

python/ribasim/ribasim/geometry/edge.py

@@ -10,7 +10,6 @@
 from pandera.dtypes import Int32
 from pandera.typing import Series
 from pandera.typing.geopandas import GeoDataFrame, GeoSeries
-from pydantic import model_validator
 from shapely.geometry import LineString, MultiLineString, Point
 
 from ribasim.input_base import SpatialTableModel
@@ -45,13 +44,6 @@ class Config:
 class EdgeTable(SpatialTableModel[EdgeSchema]):
     """Defines the connections between nodes."""
 
-    @model_validator(mode="after")
-    def empty_table(self) -> "EdgeTable":
-        if self.df is None:
-            self.df = GeoDataFrame[EdgeSchema]()
-        self.df.set_geometry("geometry", inplace=True)
-        return self
-
     def add(
         self,
         from_node: NodeData,
@@ -68,6 +60,8 @@ def add(
         edge_type = (
             "control" if from_node.node_type in SPATIALCONTROLNODETYPES else "flow"
         )
+        assert self.df is not None
+
         table_to_append = GeoDataFrame[EdgeSchema](
             data={
                 "from_node_type": pd.Series([from_node.node_type], dtype=str),
@@ -79,12 +73,10 @@ def add(
                 "subnetwork_id": pd.Series([subnetwork_id], dtype=pd.Int32Dtype()),
             },
             geometry=geometry_to_append,
+            crs=self.df.crs,
         )
 
-        if self.df is None:
-            self.df = table_to_append
-        else:
-            self.df = GeoDataFrame[EdgeSchema](pd.concat([self.df, table_to_append]))
+        self.df = GeoDataFrame[EdgeSchema](pd.concat([self.df, table_to_append]))
 
     def get_where_edge_type(self, edge_type: str) -> NDArray[np.bool_]:
         assert self.df is not None

python/ribasim/ribasim/input_base.py

@@ -1,5 +1,4 @@
 import re
-import warnings
 from abc import ABC, abstractmethod
 from collections.abc import Callable, Generator
 from contextlib import closing
@@ -39,7 +38,6 @@
 delimiter = " / "
 
 gpd.options.io_engine = "pyogrio"
-warnings.filterwarnings("ignore", category=UserWarning, module="pyogrio")
 
 context_file_loading: ContextVar[dict[str, Any]] = ContextVar(
     "file_loading", default={}
@@ -134,15 +132,6 @@ def set_filepath(self, filepath: Path) -> None:
         self.filepath = filepath
         self.model_config["validate_assignment"] = True
 
-    @abstractmethod
-    def _save(self, directory: DirectoryPath, input_dir: DirectoryPath) -> None:
-        """Save this instance to disk.
-
-        This method needs to be implemented by any class deriving from
-        FileModel.
-        """
-        raise NotImplementedError()
-
     @classmethod
     @abstractmethod
     def _load(cls, filepath: Path | None) -> dict[str, Any]:
@@ -228,9 +217,7 @@ def _load(cls, filepath: Path | None) -> dict[str, Any]:
             return {}
 
     def _save(
-        self,
-        directory: DirectoryPath,
-        input_dir: DirectoryPath,
+        self, directory: DirectoryPath, input_dir: DirectoryPath, crs: str
     ) -> None:
         # TODO directory could be used to save an arrow file
         db_path = context_file_loading.get().get("database")
@@ -239,9 +226,9 @@ def _save(
             self._write_arrow(self.filepath, directory, input_dir)
         elif db_path is not None:
             self.sort()
-            self._write_table(db_path)
+            self._write_geopackage(db_path, crs)
 
-    def _write_table(self, temp_path: Path) -> None:
+    def _write_geopackage(self, temp_path: Path, crs: str) -> None:
         """
         Write the contents of the input to a database.
 
@@ -366,16 +353,16 @@ def _from_db(cls, path: Path, table: str):
 
             return df
 
-    def _write_table(self, path: Path) -> None:
+    def _write_geopackage(self, path: Path, crs: str) -> None:
         """
-        Write the contents of the input to a database.
+        Write the contents of the input to the GeoPackage.
 
         Parameters
         ----------
         path : Path
         """
         assert self.df is not None
-        self.df.to_file(path, layer=self.tablename(), driver="GPKG", mode="a")
+        self.df.set_crs(crs).to_file(path, layer=self.tablename(), driver="GPKG")
 
 
 class ChildModel(BaseModel):
@@ -436,9 +423,9 @@ def node_ids(self) -> set[int]:
             node_ids.update(table.node_ids())
         return node_ids
 
-    def _save(self, directory: DirectoryPath, input_dir: DirectoryPath, **kwargs):
+    def _save(self, directory: DirectoryPath, input_dir: DirectoryPath, crs: str):
         for table in self._tables():
-            table._save(directory, input_dir)
+            table._save(directory, input_dir, crs)
 
     def _repr_content(self) -> str:
         """Generate a succinct overview of the content.

python/ribasim/ribasim/model.py

@@ -9,6 +9,7 @@
 import tomli
 import tomli_w
 from matplotlib import pyplot as plt
+from pandera.typing.geopandas import GeoDataFrame
 from pydantic import (
     DirectoryPath,
     Field,
@@ -39,7 +40,7 @@
     Terminal,
     UserDemand,
 )
-from ribasim.geometry.edge import EdgeTable
+from ribasim.geometry.edge import EdgeSchema, EdgeTable
 from ribasim.geometry.node import NodeTable
 from ribasim.input_base import (
     ChildModel,
@@ -86,6 +87,7 @@ class Model(FileModel):
     user_demand: UserDemand = Field(default_factory=UserDemand)
 
     edge: EdgeTable = Field(default_factory=EdgeTable)
+    crs: str
 
     @model_validator(mode="after")
     def set_node_parent(self) -> "Model":
@@ -97,6 +99,13 @@ def set_node_parent(self) -> "Model":
             setattr(v, "_parent_field", k)
         return self
 
+    @model_validator(mode="after")
+    def ensure_edge_table_is_present(self) -> "Model":
+        if self.edge.df is None:
+            self.edge.df = GeoDataFrame[EdgeSchema]()
+        self.edge.df.set_geometry("geometry", inplace=True, crs=self.crs)
+        return self
+
     @field_serializer("input_dir", "results_dir")
     def serialize_path(self, path: Path) -> str:
         return str(path)
@@ -155,7 +164,7 @@ def _save(self, directory: DirectoryPath, input_dir: DirectoryPath):
         db_path.parent.mkdir(parents=True, exist_ok=True)
         db_path.unlink(missing_ok=True)
         context_file_loading.get()["database"] = db_path
-        self.edge._save(directory, input_dir)
+        self.edge._save(directory, input_dir, self.crs)
 
         node = self.node_table()
         # Temporarily require unique node_id for #1262
@@ -165,14 +174,14 @@ def _save(self, directory: DirectoryPath, input_dir: DirectoryPath):
         node.df.set_index("node_id", drop=False, inplace=True)
         node.df.sort_index(inplace=True)
         node.df.index.name = "fid"
-        node._save(directory, input_dir)
+        node._save(directory, input_dir, self.crs)
 
         for sub in self._nodes():
-            sub._save(directory, input_dir)
+            sub._save(directory, input_dir, self.crs)
 
     def node_table(self) -> NodeTable:
         """Compute the full NodeTable from all node types."""
-        df_chunks = [node.node.df for node in self._nodes()]
+        df_chunks = [node.node.df.set_crs(self.crs) for node in self._nodes()]
         df = pd.concat(df_chunks, ignore_index=True)
         node_table = NodeTable(df=df)
         node_table.sort()
@@ -392,12 +401,6 @@ def to_xugrid(self):
             name="node_index",
         )
 
-        if node_df.crs is None:
-            # TODO: can be removed when CRS is required, #1254
-            projected = False
-        else:
-            projected = node_df.crs.is_projected
-
         grid = xugrid.Ugrid1d(
             node_x=node_df.geometry.x,
             node_y=node_df.geometry.y,
@@ -409,7 +412,7 @@ def to_xugrid(self):
                 )
             ),
             name="ribasim",
-            projected=projected,
+            projected=node_df.crs.is_projected,
             crs=node_df.crs,
         )
 

python/ribasim/tests/test_io.py

@@ -117,7 +117,7 @@ def test_sort(level_setpoint_with_minmax, tmp_path):
     model.write(tmp_path / "basic/ribasim.toml")
     # write sorts the model in place
     assert table.df.iloc[0]["greater_than"] == 5.0
-    model_loaded = ribasim.Model(filepath=tmp_path / "basic/ribasim.toml")
+    model_loaded = ribasim.Model.read(filepath=tmp_path / "basic/ribasim.toml")
     table_loaded = model_loaded.discrete_control.condition
     edge_loaded = model_loaded.edge
     assert table_loaded.df.iloc[0]["greater_than"] == 5.0
@@ -153,7 +153,9 @@ def test_roundtrip(trivial, tmp_path):
 def test_datetime_timezone():
     # Due to a pydantic issue, a time zone was added.
     # https://github.com/Deltares/Ribasim/issues/1282
-    model = ribasim.Model(starttime="2000-01-01", endtime="2001-01-01 00:00:00")
+    model = ribasim.Model(
+        starttime="2000-01-01", endtime="2001-01-01 00:00:00", crs="EPSG:28992"
+    )
     assert isinstance(model.starttime, datetime)
     assert isinstance(model.endtime, datetime)
     assert model.starttime.tzinfo is None

python/ribasim/tests/test_model.py

@@ -6,6 +6,7 @@
 import pytest
 import xugrid
 from pydantic import ValidationError
+from pyproj import CRS
 from ribasim.config import Solver
 from ribasim.geometry.edge import NodeData
 from ribasim.input_base import esc_id
@@ -139,6 +140,7 @@ def test_tabulated_rating_curve_model(tabulated_rating_curve, tmp_path):
     basin_area = tabulated_rating_curve.basin.area.df
     assert basin_area is not None
     assert basin_area.geometry.geom_type.iloc[0] == "Polygon"
+    assert basin_area.crs == CRS.from_epsg(28992)
     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(
@@ -186,6 +188,16 @@ def test_node_table(basic):
     assert df.subnetwork_id.dtype == pd.Int32Dtype()
     assert df.node_type.iloc[0] == "Basin"
     assert df.node_type.iloc[-1] == "Terminal"
+    assert df.crs == CRS.from_epsg(28992)
+
+
+def test_edge_table(basic):
+    model = basic
+    df = model.edge.df
+    assert df.geometry.is_unique
+    assert df.from_node_id.dtype == np.int32
+    assert df.subnetwork_id.dtype == pd.Int32Dtype()
+    assert df.crs == CRS.from_epsg(28992)
 
 
 def test_indexing(basic):
@@ -224,7 +236,7 @@ def test_xugrid(basic, tmp_path):
     assert isinstance(uds, xugrid.UgridDataset)
     assert uds.grid.edge_dimension == "ribasim_nEdges"
     assert uds.grid.node_dimension == "ribasim_nNodes"
-    assert uds.grid.crs is None
+    assert uds.grid.crs == CRS.from_epsg(28992)
     assert uds.node_id.dtype == np.int32
     uds.ugrid.to_netcdf(tmp_path / "ribasim.nc")
     uds = xugrid.open_dataset(tmp_path / "ribasim.nc")