Migrate from_copernicusmarine to use from_sgrid_conventions (#2434)

Co-authored-by: Erik van Sebille <e.vansebille@uu.nl>

Author: VeckoTheGecko

src/parcels/_core/fieldset.py

@@ -216,47 +216,28 @@ def from_copernicusmarine(ds: xr.Dataset):
             )
 
         ds = _rename_coords_copernicusmarine(ds)
-        grid = XGrid(
-            xgcm.Grid(
-                ds,
-                coords={
-                    "X": {
-                        "left": "lon",
-                    },
-                    "Y": {
-                        "left": "lat",
-                    },
-                    "Z": {
-                        "left": "depth",
-                    },
-                    "T": {
-                        "center": "time",
-                    },
-                },
-                autoparse_metadata=False,
-                **_DEFAULT_XGCM_KWARGS,
-            ),
-            mesh="spherical",
-        )
-
-        fields = {}
-        if "U" in ds.data_vars and "V" in ds.data_vars:
-            fields["U"] = Field("U", ds["U"], grid, XLinear)
-            fields["V"] = Field("V", ds["V"], grid, XLinear)
+        if "W" in ds.data_vars:
+            # Negate W to convert from up positive to down positive (as that's the direction of positive z)
+            ds["W"] -= ds["W"]
 
-            if "W" in ds.data_vars:
-                ds["W"] -= ds[
-                    "W"
-                ]  # Negate W to convert from up positive to down positive (as that's the direction of positive z)
-                fields["W"] = Field("W", ds["W"], grid, XLinear)
-                fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"])
-            else:
-                fields["UV"] = VectorField("UV", fields["U"], fields["V"])
-
-        for varname in set(ds.data_vars) - set(fields.keys()):
-            fields[varname] = Field(varname, ds[varname], grid, XLinear)
-
-        return FieldSet(list(fields.values()))
+        if "grid" in ds.cf.cf_roles:
+            raise ValueError(
+                "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset."
+            )
+        ds["grid"] = xr.DataArray(
+            0,
+            attrs=sgrid.Grid2DMetadata(  # use dummy *_center dimensions - this is A grid data (all defined on nodes)
+                cf_role="grid_topology",
+                topology_dimension=2,
+                node_dimensions=("lon", "lat"),
+                face_dimensions=(
+                    sgrid.DimDimPadding("x_center", "lon", sgrid.Padding.LOW),
+                    sgrid.DimDimPadding("y_center", "lat", sgrid.Padding.LOW),
+                ),
+                vertical_dimensions=(sgrid.DimDimPadding("z_center", "depth", sgrid.Padding.LOW),),
+            ).to_attrs(),
+        )
+        return FieldSet.from_sgrid_conventions(ds, mesh="spherical")
 
     def from_fesom2(ds: ux.UxDataset):
         """Create a FieldSet from a FESOM2 uxarray.UxDataset.

src/parcels/_core/utils/sgrid.py

@@ -378,7 +378,13 @@ def parse_sgrid(ds: xr.Dataset):
     xgcm_coords = {}
     for dim_dim_padding, axis in zip(dimensions, "XYZ", strict=False):
         xgcm_position = SGRID_PADDING_TO_XGCM_POSITION[dim_dim_padding.padding]
-        xgcm_coords[axis] = {"center": dim_dim_padding.dim1, xgcm_position: dim_dim_padding.dim2}
+
+        coords = {}
+        for pos, dim in [("center", dim_dim_padding.dim1), (xgcm_position, dim_dim_padding.dim2)]:
+            # only include dimensions in dataset (ignore dimensions in metadata that may not exist - e.g., due to `.isel`)
+            if dim in ds.dims:
+                coords[pos] = dim
+        xgcm_coords[axis] = coords
 
     return (ds, {"coords": xgcm_coords})
 

tests/utils/test_sgrid.py

@@ -188,7 +188,7 @@ def test_parse_sgrid_2d(grid_metadata: sgrid.Grid2DMetadata):
     """Test the ingestion of datasets in XGCM to ensure that it matches the SGRID metadata provided"""
     ds = dummy_sgrid_2d_ds(grid_metadata)
 
-    ds, xgcm_kwargs = sgrid.parse_sgrid(ds)
+    _, xgcm_kwargs = sgrid.parse_sgrid(ds)
     grid = xgcm.Grid(ds, autoparse_metadata=False, **xgcm_kwargs)
 
     for ddp, axis in zip(grid_metadata.face_dimensions, ["X", "Y"], strict=True):