FieldSet construction fails for 2D Copernicus dataset due to empty Z axis from SGRID parsing

[Dobbermann2]

Parcels version

main

Description

FieldSet.from_sgrid_conventions() fails when constructing a FieldSet from a 2D (lat,lon) Copernicus dataset.

Problem

In FieldSet.from_sgrid_conventions() the SGRID metadata is parsed using sgrid.parse_sgrid().

However, parse_sgrid returns a dictionary containing the key Z, even in absence of a vertical coordinate in the original dataset:

_, xgcm_kwargs = sgrid.parse_sgrid(ds)
print(xgcm_kwargs) # {'coords': {'X': {'right': 'lon'}, 'Y': {'right': 'lat'}, 'Z': {}}}

This coordinate dictionary is then passed to the xgcm.Grid constructor, returning a grid that now contains an incorrect Z axis.

xgcm_grid = xgcm.Grid(ds, autoparse_metadata=False, **xgcm_kwargs, **_DEFAULT_XGCM_KWARGS)
print(xgcm_grid.axes) 
'''
OrderedDict({'X': <xgcm.Axis 'X' (not periodic, boundary='fill')>
Axis Coordinates:
  * right    lon, 'Y': <xgcm.Axis 'Y' (not periodic, boundary='fill')>
Axis Coordinates:
  * right    lat, 'Z': <xgcm.Axis 'Z' (not periodic, boundary='fill')>
Axis Coordinates:, 'T': <xgcm.Axis 'T' (not periodic, boundary='fill')>
Axis Coordinates:
  * center   time})
'''

Then, a XGrid is constructed using the xgcm_grid. Because xgcm_grid.axes contains a Z axis, the XGrid constructor attempts to assert a valid depth. Since the dataset has no depth this raises a KeyError.

 if "Z" in grid.axes:
            assert_valid_depth(ds["depth"])   #KeyError: " No variable named 'depth'. Variables on the dataset include ['time', 'lat', 'lon', 'U', 'V', 'grid']"

Code sample

import xarray as xr
import glob
import parcels


ds_fields = xr.open_mfdataset(
    glob.glob(f'./data/metoffice_foam1_amm15_NWS_SSC_b202407*.nc'),
    data_vars="minimal",
    coords="minimal",
    compat="override",
)
    
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=dict(U=ds_fields["uo"], V=ds_fields["vo"]))
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)

Traceback (most recent call last):
  File "C:\Users\luukv\Documents\BONZ\development\Sandbox\main.py", line 23, in <module>  
    fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
  File "C:\Users\luukv\Documents\BONZ\development\Parcels\src\parcels\_core\fieldset.py", 
    grid = XGrid(xgcm_grid, mesh=mesh)
  File "C:\Users\luukv\Documents\BONZ\development\Parcels\src\parcels\_core\xgrid.py", line 122, in __init__
    assert_valid_depth(ds["depth"])
                       ~~^^^^^^^^^
  File "C:\Users\luukv\Documents\BONZ\development\Sandbox\.pixi\envs\default\Lib\site-packages\xarray\core\dataset.py", line 1357, in __getitem__
    raise KeyError(message) from e
KeyError: "No variable named 'depth'. Variables on the dataset include ['time', 'lat', 'lon', 'U', 'V', 'grid']"

metadata_cmems_mod_nws_phycur_anfc_1.5km2D_PT15M-i.json

[Dobbermann2]

A simple conditional check around the coord assignment should fix the issue

if coords:
    xgcm_coords[axis] = coords

in

Parcels/src/parcels/_core/utils/sgrid.py

Lines 457 to 466 in c5f6ff6

	xgcm_coords = {}
	for dim_dim_padding, axis in zip(dimensions, "XYZ", strict=False):
	xgcm_position = SGRID_PADDING_TO_XGCM_POSITION[dim_dim_padding.padding]
	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