Transposed fields in FieldSet from NetCDF

[MiriamSterl]

I'm trying to create a FieldSet from a NetCDF file:

filenames = {'U': "output_layer1.nc",
             'V': "output_layer1.nc"}
variables = {'U': 'u', 'V': 'v'}
dimensions = {'time': 't', 'lon': 'x', 'lat': 'y'}
fieldset = FieldSet.from_netcdf(filenames,variables,dimensions,mesh='flat',allow_time_extrapolation=True)
fieldset.add_periodic_halo(zonal=True,meridional=True)

When I then do fieldset.U.show() as a test, I see that the velocity field is transposed. Is there a way to read in U and V in transposed form? (Of course I can change how they are saved in the netCDF file, but I'm hoping for a more elegant solution.) I saw that it is possible for a Field (setting transpose=True), but not sure how it works for a FieldSet.

[VeckoTheGecko]

Its a bit difficult to diagnose the problem without having the file output_layer1.nc to look at. From my understanding, the transpose=True flag is mainly to do with the internal ordering of the dataset (ie. how the NetCDF file appears on disk). Parcels is sensitive to this ordering in that it needs to be ["time", "y", "x"] as mentioned here in the documentation. Normally incorrect ordering is picked up by parcels on read-in, spitting back an informative error message, but your domain could be perfectly square and hence snuck through.

If you plot by loading the dataset in xarray, do you get the same problem?

import xarray as xr
ds = xr.load_dataset('test.nc')
ds.u.plot(x='LON', y='LAT') # LON and LAT here are the dimensions, not the meshgrid variables lon and lat

If yes, then the issue was with how output_layer1.nc was created.
If no (which I think is more likely), then the issue is with how parcels is reading it in. This can be resolved by transposing the dataset. This could be as simple as setting transpose=True for Fieldset.from_netcdf:

fieldset = FieldSet.from_netcdf(filenames,variables,dimensions,mesh='flat',allow_time_extrapolation=True, transpose=True)

or by using xarray to transpose and save back to .nc, or load directly from the xarray dataset using Fieldset.from_xarray_dataset. I have example code in a related issue that shows how this can be done using xarray.

Feel free to upload the dataset to WeTransfer and post the link if you'd like me to look at the dataset.

[MiriamSterl]

Thank you! My domain is indeed square, so probably that's why Parcels isn't picking it up. I tried to add transpose=True in the Fieldset.from_netcdf function, but that didn't do the trick for me. I've now just saved the velocity fields in transposed form in the .nc file and now it works fine.