Merge pull request #401 from DHI/netcdf-3d

Support for 3d via xarray (NetCDF etc)

modelskill/model/grid.py

@@ -148,7 +148,7 @@ def _extract_point(
         given a PointObservation. No time interpolation is done!"""
         method: str = spatial_method or "linear"
 
-        x, y = observation.x, observation.y
+        x, y, z = observation.x, observation.y, observation.z
         if (x is None) or (y is None):
             raise ValueError(
                 f"PointObservation '{observation.name}' cannot be used for extraction "
@@ -163,7 +163,12 @@ def _extract_point(
         assert isinstance(self.data, xr.Dataset)
 
         # TODO: avoid runtrip to pandas if possible (potential loss of metadata)
-        ds = self.data.interp(coords=dict(x=x, y=y), method=method)  # type: ignore
+        if "z" in self.data.dims and z is not None:
+            ds = self.data.interp(
+                coords=dict(x=x, y=y, z=z), method=method  # type: ignore
+            )
+        else:
+            ds = self.data.interp(coords=dict(x=x, y=y), method=method)  # type: ignore
         # TODO: exclude aux cols in dropna
         df = ds.to_dataframe().drop(columns=["x", "y"]).dropna()
         if len(df) == 0:

modelskill/utils.py

@@ -17,6 +17,8 @@
     "north": "y",
     "x": "x",
     "y": "y",
+    "z": "z",
+    "depth": "z",
 }
 TIME_COORDINATE_NAME_MAPPING = {
     "t": "time",

tests/integration/test_integration_grid.py

@@ -0,0 +1,30 @@
+import modelskill as ms
+import pytest
+import xarray as xr
+
+def test_extract_point_from_3d():
+    mr = ms.GridModelResult(
+        "tests/testdata/cmems_mod_med_phy-sal_anfc_4.2km-3D_PT1H-m_1705916517624.nc",
+        item="so",
+        name="MedSea",
+    )
+
+    point_ds = xr.open_dataset("tests/testdata/aegean_sea_salinity_ts.nc")
+
+    # TODO use x,y,z without explicitly setting them (NetCDF has them as coordinates)
+    obs = ms.PointObservation(
+        point_ds,
+        x=float(point_ds.longitude),
+        y=float(point_ds.latitude),
+        z=float(point_ds.depth),
+        item="so",
+    )
+
+    cmp = ms.match(obs=obs, mod=mr, spatial_method="nearest")
+    assert cmp.quantity.name == "Salinity"
+    assert cmp.quantity.unit == "0.001"
+
+    sc = cmp.score()
+
+    # "Observed" data is extracted from the 3D model result, so the score should be 0.0
+    assert sc["MedSea"] == pytest.approx(0.0)

tests/model/test_grid.py

@@ -218,3 +218,31 @@ def test_grid_with_directional_data_with_cf_metadata_is_directional_by_default()
         "tests/testdata/SW/CMEMS_DutchCoast_2017-10-28.nc", item="VMDR"
     )
     assert mr.quantity.is_directional
+
+
+def test_extract_point_from_3d():
+    mr = ms.GridModelResult(
+        "tests/testdata/cmems_mod_med_phy-sal_anfc_4.2km-3D_PT1H-m_1705916517624.nc",
+        item="so",
+        name="MedSea",
+    )
+
+    point_ds = xr.open_dataset("tests/testdata/aegean_sea_salinity_ts.nc")
+
+    # TODO use x,y,z without explicitly setting them (NetCDF has them as coordinates)
+    obs = ms.PointObservation(
+        point_ds,
+        x=float(point_ds.longitude),
+        y=float(point_ds.latitude),
+        z=float(point_ds.depth),
+        item="so",
+    )
+
+    cmp = ms.match(obs=obs, mod=mr, spatial_method="nearest")
+    assert cmp.quantity.name == "Salinity"
+    assert cmp.quantity.unit == "0.001"
+
+    sc = cmp.score()
+
+    # "Observed" data is extracted from the 3D model result, so the score should be 0.0
+    assert sc["MedSea"] == pytest.approx(0.0)