Error setting up a 2D fieldset from a curvilinear grid.

[imariaguiar]

Hello there,
I'm trying to set up a 2D fieldset with MOi hydrodynamics (psy4v3r1 daily) data following their tutorial (https://github.com/OceanParcels/UtrechtTeam/wiki/MOi-hydrodynamics-%28psy4v3r1-daily%29). I'm getting an error, and my guess is that it may be related to the coordinate file or the indices of the x and y dimensions. U and V data are cropped for a specific region, just like the coordinate file, but it seems like the indices doesn't match. I'm not sure if that's the issue, but I couldn't figure out how to fix it.
If anyone has had a similar problem, I'd appreciate some help.

Here's my code:

data_path = 'D:/data/moi_glo12/daily_mean/'
mesh_mask = data_path + 'coordinates.nc'

ufiles = sorted(glob(data_path+'psy4v3r1-U-.nc'))
vfiles = sorted(glob(data_path+'psy4v3r1-U-.nc'))

filenames = {'U': {'lon': mesh_mask, 'lat': mesh_mask, 'data': ufiles},
'V': {'lon': mesh_mask, 'lat': mesh_mask, 'data': vfiles}}
variables = {'U': 'vozocrtx', 'V': 'vomecrty'}
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'time': 'time_counter'},
'V': {'lon': 'glamf', 'lat': 'gphif', 'time': 'time_counter'}}

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True)

---

IndexError Traceback (most recent call last)
Cell In[80], line 1
----> 1 fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True)

File ~\anaconda3\envs\py3_parcels\lib\site-packages\parcels\fieldset.py:416, in FieldSet.from_netcdf(cls, filenames, variables, dimensions, indices, fieldtype, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, chunksize, **kwargs)
414 dFiles = fields[procvar].dataFiles
415 break
--> 416 fields[var] = Field.from_netcdf(paths, (var, name), dims, inds, grid=grid, mesh=mesh, timestamps=timestamps,
417 allow_time_extrapolation=allow_time_extrapolation,
418 time_periodic=time_periodic, deferred_load=deferred_load,
419 fieldtype=fieldtype, chunksize=varchunksize, dataFiles=dFiles, **kwargs)
421 u = fields.pop('U', None)
422 v = fields.pop('V', None)

File ~\anaconda3\envs\py3_parcels\lib\site-packages\parcels\field.py:390, in Field.from_netcdf(cls, filenames, variable, dimensions, indices, grid, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, **kwargs)
385 raise RuntimeError('Multiple files given but no time dimension specified')
387 if grid is None:
388 # Concatenate time variable to determine overall dimension
389 # across multiple files
--> 390 time, time_origin, timeslices, dataFiles = cls.collect_timeslices(timestamps, data_filenames,
391 _grid_fb_class, dimensions,
392 indices, netcdf_engine, netcdf_decodewarning)
393 grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
394 grid.timeslices = timeslices

File ~\anaconda3\envs\py3_parcels\lib\site-packages\parcels\field.py:259, in Field.collect_timeslices(timestamps, data_filenames, _grid_fb_class, dimensions, indices, netcdf_engine, netcdf_decodewarning)
256 time = time_origin.reltime(time)
258 if not np.all((time[1:] - time[:-1]) > 0):
--> 259 id_not_ordered = np.where(time[1:] < time[:-1])[0][0]
260 raise AssertionError(
261 'Please make sure your netCDF files are ordered in time. First pair of non-ordered files: %s, %s'
262 % (dataFiles[id_not_ordered], dataFiles[id_not_ordered + 1]))
263 return time, time_origin, timeslices, dataFiles

IndexError: index 0 is out of bounds for axis 0 with size 0

[erikvansebille]

Hi @imariaguiar, thanks for reporting. The error seems to related to the time dimension, but I can't really see from the error message what exactly goed wrong. Could you share what the time dimension of these files is like?

[imariaguiar]

Hello Erik, Thank you for getting back to me.

Regarding the U and V files, for each month of daily data, I divided it into two files because it was a bit large to download the entire month through Opendap. So, the time dimensions became as follows:

Dimensions:
(y: 240 x: 959 deptht: 28 time_counter: 15)
Dimensions:
(y: 240 x: 959 deptht: 28 time_counter: 16)
This is an example for the two files of January 2007, but I'm trying to read the entire data period.

And regarding the coordinate data, here are the dimensions:

xarray.Dataset
Dimensions:
y: 239 x: 959 z: 50 t: 1
Coordinates: (0)
Data variables: (21)
Indexes: (0)
Attributes:
file_name :
mesh_hgr.nc
TimeStamp :
03/02/2016 10:24:55

d['glamf'].t.values
array([0], dtype=int64)
d['gphif'].t.values
array([0], dtype=int64)

[erikvansebille]

Hmm, and what if you provide only one file (with 15 or 16 snapshots) as input? Does creating a fields work in that case?

[imariaguiar]

Hi Erik, I found where the problem was. There was only one file with incorrect dates, and it was disrupting the sequence of dates. Thank you very much for your help!

But taking advantage of the question about creating the fieldset, and to conclude the discussion, is it normal that when I call the FieldSet.from_nemo function and then request to plot the 'U' field, I get a plot of a rectilinear grid? I recently started importing the curvilinear grid into Parcels, so I'm not sure if it's correct this way.
Here is an example of the plot I am getting when I call fieldset.U.show()

[erikvansebille]

Dear @imariaguiar, since #1397 we don't support the built-in plotting scripts anymore (because they were indeed a bit buggy and very cumbersome to maintain), so I can't help you. Better to simply write your own plotting scripts using matplotlib, cartopy, etc.