Remove model.to_fews, add CF attrs to netcdf input (#2938)

Fixes #2935

I noticed we didn't add any CF metadata to the NetCDF input files, this
fixes that and tests it.

Author: visr

python/ribasim/ribasim/input_base.py

@@ -630,11 +630,25 @@ def _write_netcdf(self, filepath: Path) -> None:
             raise err
         elif "time" in cols:
             ds = self.df.set_index(["time", "node_id"]).to_xarray()
+            ds["time"].attrs.update(
+                {"standard_name": "time", "axis": "T", "long_name": "time"}
+            )
         elif "node_id" in cols:
             ds = self.df.set_index(["node_id"]).to_xarray()
+            ds["node_id"].attrs.update(
+                {"cf_role": "timeseries_id", "long_name": "station identification code"}
+            )
         else:
             raise err
 
+        ds.attrs.update(
+            {
+                "Conventions": "CF-1.12",
+                "title": "Ribasim model input",
+                "references": "https://ribasim.org",
+            }
+        )
+
         # Write to NetCDF file
         ds.to_netcdf(filepath)
 

python/ribasim/ribasim/model.py

@@ -63,7 +63,6 @@
 )
 from ribasim.utils import (
     MissingOptionalModule,
-    _add_cf_attributes,
     _concat,
     _link_lookup,
     _node_lookup,
@@ -895,120 +894,3 @@ def _add_allocation(self, uds):
         )
 
         return uds
-
-    def to_fews(
-        self,
-        region_home: str | PathLike[str],
-        add_network: bool = True,
-        add_results: bool = True,
-    ) -> None:
-        """
-        Write the model network and results into files used by Delft-FEWS.
-
-        ** Warning: This method is experimental and is likely to change. **
-
-        To run this method, the model needs to be written to disk, and have results.
-        The Node, Link and Basin / area tables are written to shapefiles in the REGION_HOME/Config directory.
-        The results are written to NetCDF files in the REGION_HOME/Modules directory.
-        The netCDF files are NetCDF4 with CF-conventions.
-
-        Parameters
-        ----------
-        region_home: str | PathLike[str]
-            Path to the Delft-FEWS REGION_HOME directory.
-        add_network: bool, optional
-            Write shapefiles representing the network, enabled by default.
-        add_results: bool, optional
-            Write the results to NetCDF files, enabled by default.
-        """
-        region_home = DirectoryPath(region_home)
-        if add_network:
-            self._network_to_fews(region_home)
-        if add_results:
-            self._results_to_fews(region_home)
-
-    def _network_to_fews(self, region_home: DirectoryPath) -> None:
-        """Write the Node and Link tables to shapefiles for use in Delft-FEWS."""
-        df_link = self.link.df
-        df_node = self.node.df
-        assert df_link is not None
-        assert df_node is not None
-
-        df_basin_area = self.basin.area.df
-        if df_basin_area is None:
-            # Fall back to the Basin points if the area polygons are not set
-            df_basin_area = df_node[df_node["node_type"] == "Basin"]
-
-        network_dir = region_home / "Config/MapLayerFiles/{ModelId}"
-        network_dir.mkdir(parents=True, exist_ok=True)
-        link_path = network_dir / "{ModelId}Links.shp"
-        node_path = network_dir / "{ModelId}Nodes.shp"
-        basin_area_path = network_dir / "{ModelId}Areas.shp"
-
-        with warnings.catch_warnings():
-            warnings.filterwarnings(
-                "ignore", "Normalized/laundered field name", RuntimeWarning
-            )
-            warnings.filterwarnings(
-                "ignore",
-                "Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.",
-                UserWarning,
-            )
-            df_link.to_file(link_path)
-            df_node.to_file(node_path)
-            df_basin_area.to_file(basin_area_path)
-
-    def _results_to_fews(self, region_home: DirectoryPath) -> None:
-        """Convert the model results to NetCDF with CF-conventions for importing into Delft-FEWS."""
-        # Delft-FEWS doesn't support our UGRID from `model.to_xugrid` yet,
-        # so we convert Arrow to regular CF-NetCDF4.
-
-        ext = self.results_extension
-        basin_path = self.results_path / f"basin{ext}"
-        flow_path = self.results_path / f"flow{ext}"
-        concentration_path = self.results_path / f"concentration{ext}"
-
-        if ext == ".arrow":
-            basin_df = pd.read_feather(basin_path)
-            flow_df = pd.read_feather(flow_path)
-        else:
-            basin_df = xr.open_dataset(basin_path).to_dataframe().reset_index()
-            flow_df = xr.open_dataset(flow_path).to_dataframe().reset_index()
-
-        ds_basin = basin_df.set_index(["time", "node_id"]).to_xarray()
-        _add_cf_attributes(ds_basin, timeseries_id="node_id")
-        ds_basin["level"].attrs.update({"units": "m"})
-        ds_basin["storage"].attrs.update({"units": "m3"})
-        ds_basin["relative_error"].attrs.update({"units": "1"})
-
-        flow_rate_variables = [
-            "inflow_rate",
-            "outflow_rate",
-            "storage_rate",
-            "precipitation",
-            "evaporation",
-            "drainage",
-            "infiltration",
-            "balance_error",
-        ]
-        for var in flow_rate_variables:
-            ds_basin[var].attrs.update({"units": "m3 s-1"})
-
-        ds_flow = flow_df.set_index(["time", "link_id"]).to_xarray()
-        _add_cf_attributes(ds_flow, timeseries_id="link_id")
-        ds_flow["flow_rate"].attrs.update({"units": "m3 s-1"})
-
-        results_dir = region_home / "Modules/ribasim/{ModelId}/work/results"
-        results_dir.mkdir(parents=True, exist_ok=True)
-        ds_basin.to_netcdf(results_dir / "basin.nc")
-        ds_flow.to_netcdf(results_dir / "flow.nc")
-
-        if concentration_path.is_file():
-            if ext == ".arrow":
-                df = pd.read_feather(concentration_path)
-            else:
-                df = xr.open_dataset(concentration_path).to_dataframe().reset_index()
-            ds = df.set_index(["time", "node_id", "substance"]).to_xarray()
-            _add_cf_attributes(ds, timeseries_id="node_id", realization="substance")
-            ds["concentration"].attrs.update({"units": "g m-3"})
-            ds.to_netcdf(results_dir / "concentration.nc")

python/ribasim/ribasim/utils.py

@@ -90,48 +90,6 @@ def _concat(dfs, **kwargs):
         return pd.concat(dfs, **kwargs)
 
 
-def _add_cf_attributes(ds, timeseries_id: str, realization: str | None = None) -> None:
-    """
-    Add CF attributes to an xarray.Dataset.
-
-    Parameters
-    ----------
-    ds : xarray.Dataset
-        The dataset to which CF attributes will be added.
-    timeseries_id : str
-        The name of the variable that identifies the timeseries.
-    realization : str | None, optional
-        The name of the variable representing realizations (e.g., "substance"), if applicable.
-
-    Returns
-    -------
-    None
-    """
-    ds.attrs.update(
-        {
-            "Conventions": "CF-1.12",
-            "title": "Ribasim model results",
-            "references": "https://ribasim.org",
-        }
-    )
-    ds["time"].attrs.update({"standard_name": "time", "axis": "T", "long_name": "time"})
-    ds[timeseries_id].attrs.update(
-        {"cf_role": "timeseries_id", "long_name": "station identification code"}
-    )
-    if realization:
-        # Use realization as the standard name as recommended by ECMWF.
-        # axis = "E" is not currently enabled since it seemed to confuse Delft-FEWS.
-        # https://confluence.ecmwf.int/display/COPSRV/Metadata+recommendations+for+encoding+NetCDF+products+based+on+CF+convention#MetadatarecommendationsforencodingNetCDFproductsbasedonCFconvention-3.4Realizationdiscretecoordinates
-        ds[realization].attrs.update(
-            {
-                "standard_name": "realization",
-                "units": "1",
-                "long_name": "substance name",
-            }
-        )
-    return ds
-
-
 class UsedIDs(BaseModel):
     """A helper class to manage globally unique node IDs.
 

python/ribasim/tests/test_io.py

@@ -8,6 +8,7 @@
 import pytest
 import ribasim
 import tomli
+import xarray as xr
 from pandas import DataFrame
 from pandas.testing import assert_frame_equal
 from pydantic import ValidationError
@@ -387,7 +388,17 @@ def test_roundtrip_netcdf(tabulated_rating_curve_control, tmp_path):
     model2 = Model.read(model1dir / "ribasim.toml")
 
     assert (model1dir / "input/database.gpkg").is_file()
-    assert (model1dir / "input/basin-state.nc").is_file()
+    basin_state_path = model1dir / "input/basin-state.nc"
+    assert basin_state_path.is_file()
+
+    # check CF attributes in the NetCDF file
+    with xr.open_dataset(basin_state_path) as ds:
+        assert ds.attrs["Conventions"] == "CF-1.12"
+        assert ds.attrs["title"] == "Ribasim model input"
+        assert ds.attrs["references"] == "https://ribasim.org"
+        node_attrs = ds["node_id"].attrs
+        assert node_attrs["cf_role"] == "timeseries_id"
+        assert node_attrs["long_name"] == "station identification code"
 
     df1 = model1.basin.state.df
     df2 = model2.basin.state.df

python/ribasim/tests/test_model.py

@@ -213,29 +213,6 @@ def test_to_xugrid(model, tmp_path):
         model.to_xugrid(add_flow=True, add_allocation=True)
 
 
-@pytest.mark.parametrize(
-    "model",
-    [basic_model(), outlet_model(), pid_control_equation_model(), trivial_model()],
-)
-def test_to_fews(model, tmp_path):
-    region_home = tmp_path
-    network_dir = region_home / "Config/MapLayerFiles/{ModelId}"
-
-    with pytest.raises(FileNotFoundError, match="Model must be written to disk"):
-        model.to_fews(region_home)
-
-    model.write(tmp_path / "model/ribasim.toml")
-    model.to_fews(region_home, add_results=False)
-    assert (network_dir / "{ModelId}Links.dbf").is_file()
-    assert (network_dir / "{ModelId}Links.shp").is_file()
-    assert (network_dir / "{ModelId}Nodes.dbf").is_file()
-    assert (network_dir / "{ModelId}Nodes.shp").is_file()
-
-    # Cannot test results=True without results
-    with pytest.raises(FileNotFoundError, match=r"Cannot find basin_state.nc"):
-        model.to_fews(region_home, add_results=True)
-
-
 def test_to_crs(bucket: Model):
     model = bucket