convert.mitgcm_to_sgrid breaks on datasets_circulation_models["ds_MITgcm_netcdf"]

[erikvansebille]

As mentioned, in #2483, the following code currently breaks:

ds = datasets_circulation_models["ds_MITgcm_netcdf"]
ds = ds.rename({"X": "XG", "Y": "YG", "Z": "depth", "T": "time"})
coords = ds[["XG", "YG", "depth", "time"]]

ds = convert.mitgcm_to_sgrid(fields={"U": ds["U"], "V": ds["V"]}, coords=coords)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds)

with error

    def test_convert_mitgcm_offsets():
        ds = datasets_circulation_models["ds_MITgcm_netcdf"]
        ds = ds.rename({"X": "XG", "Y": "YG", "Z": "depth", "T": "time"})
        coords = ds[["XG", "YG", "depth", "time"]]
    
        ds = convert.mitgcm_to_sgrid(fields={"U": ds["U"], "V": ds["V"]}, coords=coords)
>       fieldset = parcels.FieldSet.from_sgrid_conventions(ds)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

tests/test_convert.py:50: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
src/parcels/_core/fieldset.py:311: in from_sgrid_conventions
    fields["U"] = Field("U", ds["U"], grid, XLinear)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
src/parcels/_core/field.py:106: in __init__
    assert_all_field_dims_have_axis(data, grid.xgcm_grid)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

da = <xarray.DataArray 'U' (time: 13, depth: 90, lat: 60, Xp1: 31)> Size: 9MB
array([[[[0.8566365 , 0.50562197, 0.5710538 ,... 248B 0.0 136.7 273.3 ... 3.827e+03 3.963e+03 4.1e+03
Attributes:
    units:        m/s
    coordinates:  XU YU RC iter
xgcm_grid = <xgcm.Grid>
X Axis (not periodic, boundary='fill'):
  * left     lon
Y Axis (not periodic, boundary='fill'):
  * left ...'):
  * center   depth --> left
  * left     depth --> center
T Axis (not periodic, boundary='fill'):
  * center   time

    def assert_all_field_dims_have_axis(da: xr.DataArray, xgcm_grid: xgcm.Grid) -> None:
        ax_dims = [(get_axis_from_dim_name(xgcm_grid.axes, dim), dim) for dim in da.dims]
        for dim in ax_dims:
            if dim[0] is None:
>               raise ValueError(
                    f'Dimension "{dim[1]}" has no axis attribute. '
                    f'HINT: You may want to add an {{"axis": A}} to your DataSet["{dim[1]}"], where A is one of "X", "Y", "Z" or "T"'
                )
E               ValueError: Dimension "Xp1" has no axis attribute. HINT: You may want to add an {"axis": A} to your DataSet["Xp1"], where A is one of "X", "Y", "Z" or "T"

src/parcels/_core/xgrid.py:56: ValueError

This is despite an axis attribute X being added to Xp1 in

Parcels/src/parcels/convert.py

Lines 335 to 336 in f6e1e53

	ds = _set_axis_attrs(ds, _MITGCM_AXIS_VARNAMES)
	ds = _maybe_swap_depth_direction(ds)

Parcels/src/parcels/convert.py

Lines 45 to 59 in f6e1e53

	_MITGCM_AXIS_VARNAMES = {
	"XC": "X",
	"XG": "X",
	"Xp1": "X",
	"lon": "X",
	"YC": "Y",
	"YG": "Y",
	"Yp1": "Y",
	"lat": "Y",
	"Zu": "Z",
	"Zl": "Z",
	"Zp1": "Z",
	"time": "T",
	}

@VeckoTheGecko, could you have a look what's going on and how to fix?

[VeckoTheGecko]

The error message is incorrect.

---

import parcels
import parcels.convert as convert
from parcels._datasets.structured.circulation_models import datasets as datasets_circulation_models
import cf_xarray

ds = datasets_circulation_models["ds_MITgcm_netcdf"]
ds = ds.rename({"X": "XG", "Y": "YG", "Z": "depth", "T": "time"})
coords = ds[["XG", "YG", "depth", "time"]]
breakpoint()
ds = convert.mitgcm_to_sgrid(fields={"U": ds["U"], "V": ds["V"]}, coords=coords)
breakpoint()
fieldset = parcels.FieldSet.from_sgrid_conventions(ds)

At the first breakpoint

>>> ds
Dimensions:  (time: 13, depth: 90, YG: 60, Xp1: 31, Yp1: 61, XG: 30, Zl: 90)
Coordinates:
  * time     (time) float64 104B 0.0 1.0 2.0 3.0 4.0 ... 8.0 9.0 10.0 11.0 12.0
  * depth    (depth) float64 720B -25.0 -80.9 -136.8 ... -4.944e+03 -5e+03
  * YG       (YG) float64 480B 500.0 576.3 652.5 ... 4.847e+03 4.924e+03 5e+03
  * Xp1      (Xp1) float64 248B 0.0 136.7 273.3 ... 3.827e+03 3.963e+03 4.1e+03
  * Yp1      (Yp1) float64 488B 0.0 75.0 150.0 ... 4.35e+03 4.425e+03 4.5e+03
  * XG       (XG) float64 240B 500.0 655.2 810.3 ... 4.69e+03 4.845e+03 5e+03
  * Zl       (Zl) float64 720B 0.0 -50.56 -101.1 ... -4.449e+03 -4.5e+03
Data variables:
    U        (time, depth, YG, Xp1) float32 9MB 0.2078 0.4842 ... 0.5442 0.1105
    V        (time, depth, Yp1, XG) float32 9MB 0.6604 0.7276 ... 0.7548 0.8881
    W        (time, Zl, YG, XG) float32 8MB 0.5705 0.2498 ... 0.7213 0.7758
    Temp     (time, depth, YG, XG) float32 8MB 0.3924 0.4062 ... 0.8701 0.2597
>>> ds.cf.axes
{}

At the second breakpoint

>>> ds
Dimensions:  (time: 13, depth: 90, lat: 60, Xp1: 31, Yp1: 61, lon: 30)
Coordinates:
  * time     (time) float64 104B 0.0 1.0 2.0 3.0 4.0 ... 8.0 9.0 10.0 11.0 12.0
  * depth    (depth) float64 720B -5e+03 -4.944e+03 -4.888e+03 ... -80.9 -25.0
  * lat      (lat) float64 480B 500.0 576.3 652.5 ... 4.847e+03 4.924e+03 5e+03
  * Xp1      (Xp1) float64 248B 0.0 136.7 273.3 ... 3.827e+03 3.963e+03 4.1e+03
  * Yp1      (Yp1) float64 488B 0.0 75.0 150.0 ... 4.35e+03 4.425e+03 4.5e+03
  * lon      (lon) float64 240B 500.0 655.2 810.3 ... 4.69e+03 4.845e+03 5e+03
Data variables:
    U        (time, depth, lat, Xp1) float32 9MB 0.5809 0.0794 ... 0.9262 0.2308
    V        (time, depth, Yp1, lon) float32 9MB 0.09289 0.01039 ... 0.3083
    grid     int64 8B 0
>>> ds.cf.axes
{'X': ['Xp1', 'lon'], 'Y': ['Yp1', 'lat'], 'Z': ['depth'], 'T': ['time']}
>>> ds.grid.attrs
{'cf_role': 'grid_topology', 'topology_dimension': 2, 'node_dimensions': 'lon lat', 'face_dimensions': 'XC:lon (padding:high) YC:lat (padding:high)', 'node_coordinates': 'lon lat', 'vertical_dimensions': 'depth:depth (padding:high)'}

I'm seeing XC and YC mentioned in the SGRID metadata but not present in the data - is that it?

[VeckoTheGecko]

The error message is incorrect.

I think this happened because before we were pulling from this axis metadata, where now I think we're pulling from the SGRID conventions

[erikvansebille]

OK, so any way to fix it? Should we rename X to XC? If xo, under the hood, or ask users to do that?