FieldSet.from_netcdf when the NetCDF files have depth dimension

[albertotb]

Parcels version

3.1.2

Description

I'm trying to update from OceanParcels v2.3.0 to v3.1.2. I was following the example to create a FieldSet from netCDF files: https://docs.oceanparcels.org/en/latest/examples/parcels_tutorial.html#Reading-in-data-from-arbritrary-NetCDF-files

However, my files do have a depth dimension:

Otherwise, they have a similar structure to the example ones:

However, when trying to create a FieldSet I'm getting this error:

fieldset = parcels.FieldSet.from_netcdf(
    {
        "U": ["./test/data/uv/2020-03-20.nc", "./test/data/uv/2020-03-21.nc"],
        "V": ["./test/data/uv/2020-03-20.nc", "./test/data/uv/2020-03-21.nc"],
    },
    variables={"U": "uo", "V": "vo"},
    dimensions={
        "time": "time",
        "depth": "depth",
        "lat": "latitude",
        "lon": "longitude",
    },
)

AssertionError: Field uo expecting a data shape of [tdim=2, zdim=26, ydim=241, xdim=241] but got shape (1, 26, 241, 241). Flag transpose=True could help to reorder the data.

I tried the suggestion of using transpose=True, but no luck. For some reason, it seems to expect 2 values in the time dimension, but each NetCDF file only has one, you need to combine them both. I went ahead and removed the depth dimension from the files:

And that seems to work:

fieldset = parcels.FieldSet.from_netcdf(
    {
        "U": [
            "./test/data/uv/2020-03-20_nodepth.nc",
            "./test/data/uv/2020-03-21_nodepth.nc",
        ],
        "V": [
            "./test/data/uv/2020-03-20_nodepth.nc",
            "./test/data/uv/2020-03-21_nodepth.nc",
        ],
    },
    variables={"U": "uo", "V": "vo"},
    dimensions={
        "time": "time",
        "lat": "latitude",
        "lon": "longitude",
    },
)

I'm not sure if this is a bug, but this same code used to work in OceanParcels v2.3.0 with both depth and time dimension.

Code sample

I know that to run the example you need the data files. If required, I can try and upload them somewhere, because they are too big to be attached directly in Github.

fieldset = parcels.FieldSet.from_netcdf(
    {
        "U": ["./test/data/uv/2020-03-20.nc", "./test/data/uv/2020-03-21.nc"],
        "V": ["./test/data/uv/2020-03-20.nc", "./test/data/uv/2020-03-21.nc"],
    },
    variables={"U": "uo", "V": "vo"},
    dimensions={
        "time": "time",
        "depth": "depth",
        "lat": "latitude",
        "lon": "longitude",
    },
)

AssertionError: Field uo expecting a data shape of [tdim=2, zdim=26, ydim=241, xdim=241] but got shape (1, 26, 241, 241). Flag transpose=True could help to reorder the data.

[VeckoTheGecko]

I would need to have a better look at the actual datasets that you're using in order to get a better idea. Would you be able to provide me the data files so that I can have a look? (via https://wetransfer.com/ perhaps)

I think it might be something to do with the time dimension, but that's all I can gather without the data files to look at.

PS: Just migrating this to a Q&A since this is more of a usage question than a bug :)

[albertotb]

@VeckoTheGecko Thanks, here you have the data. I labelled it a bug because it used to work, but maybe it needs to be changed in newer versions:

data.zip

(I reduced the lat/lons so it fits here compressed)

[albertotb]

This is weird. With exactly the same data files, this works:

with xr.open_mfdataset(
    [
        "./test/data/uv/2020-03-20.nc",
        "./test/data/uv/2020-03-21.nc",
    ],
    combine="by_coords",
) as ds_:
    fieldset = parcels.FieldSet.from_xarray_dataset(
        ds_,
        variables={"U": "uo", "V": "vo"},
        dimensions={
            "time": "time",
            "depth": "depth",
            "lat": "latitude",
            "lon": "longitude",
        },
    )

But this doesn't

fieldset = parcels.FieldSet.from_netcdf(
    {
        "U": ["./test/data/uv/2020-03-20.nc", "./test/data/uv/2020-03-21.nc"],
        "V": ["./test/data/uv/2020-03-20.nc", "./test/data/uv/2020-03-21.nc"],
    },
    variables={"U": "uo", "V": "vo"},
    dimensions={
        "time": "time",
        "depth": "depth",
        "lat": "latitude",
        "lon": "longitude",
    },
)

I (maybe incorrectly) assume that both things should be equivalent, right?

[VeckoTheGecko]

Hmm, yes this is a bug. For some reason having the depth dimension there is causing an issue.

I (maybe incorrectly) assume that both things should be equivalent, right?

Yes, these both should be equivalent.

Internally in v3 of Parcels we don't use mf_dataset to open the NetCDF files. In v4 we will be working very closely with xarray for fieldset loading. Going to make an issue so that we can track this. I think it would be good if you go ahead with using from_xarray_dataset for the timebeing.

Thanks for the report and minimal examples!