Error concatenation

[imariaguiar]

Question

Hello, I'm creating a fieldset with data from Mercator Ocean International (GLORYS reanalysis) which has a curvilinear grid. When I create the fieldset, I do not receive any errors, but when I run a simulation after the particleset and execution, I receive the following error:

PSET AND EXECUTION CODE:

start_date = datetime.date(1993, 1, 2)
npart = 10
pset = ParticleSet.from_list(fieldset=fieldset, 
                             pclass=ScipyParticle,
                             lon=list(np.linspace(-38.9, -38.5, npart)) + list(np.linspace(-43.2, -42.9, npart)),  
                             lat=list(np.full(npart, -13.0)) + list(np.full(npart, -22.9)),  
                             time=start_date)

output_file = pset.ParticleFile(name='bg_bts_advection.zarr', outputdt=delta(hours=24))

pset.execute(AdvectionRK4,
             runtime=delta(days=365),
             dt=delta(hours=6),
             output_file=output_file)

Error: ValueError                                Traceback (most recent call last)
Cell In[4], line 16
      8 pset = ParticleSet.from_list(fieldset=fieldset, 
      9                              pclass=ScipyParticle,
     10                              lon=list(np.linspace(-38.9, -38.5, npart)) + list(np.linspace(-43.2, -42.9, npart)),  
     11                              lat=list(np.full(npart, -13.0)) + list(np.full(npart, -22.9)),  
     12                              time=start_date)
     14 output_file = pset.ParticleFile(name='bg_bts_advection.zarr', outputdt=delta(hours=24))
---> 16 pset.execute(AdvectionRK4,
     17              runtime=delta(days=365),
     18              dt=delta(hours=6),
     19              output_file=output_file)

File ~\anaconda3\envs\parcels\Lib\site-packages\parcels\particleset.py:932, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt)
    930     next_output = time + np.infty if dt > 0 else - np.infty
    931 next_callback = time + callbackdt if dt > 0 else time - callbackdt
--> 932 next_input = self.fieldset.computeTimeChunk(time, np.sign(dt)) if self.fieldset is not None else np.inf
    934 if output_file:
    935     logger.info(f'Output files are stored in {output_file.fname}.')

File ~\anaconda3\envs\parcels\Lib\site-packages\parcels\fieldset.py:1198, in FieldSet.computeTimeChunk(self, time, dt)
   1196             fb.close()
   1197         fb = None
-> 1198     data = f.computeTimeChunk(data, tind)
   1199 data = f.rescale_and_set_minmax(data)
   1201 if (isinstance(f.data, DeferredArray)):

File ~\anaconda3\envs\parcels\Lib\site-packages\parcels\field.py:1488, in Field.computeTimeChunk(self, data, tindex)
   1486 elif len(buffer_data.shape) == 3:
   1487     buffer_data = lib.reshape(buffer_data, sum(((buffer_data.shape[0], 1, ), buffer_data.shape[1:]), ()))
-> 1488 data = self.data_concatenate(data, buffer_data, tindex)
   1489 self.filebuffers[tindex] = filebuffer
   1490 return data

File ~\anaconda3\envs\parcels\Lib\site-packages\parcels\field.py:1446, in Field.data_concatenate(self, data, data_to_concat, tindex)
   1444 lib = np if isinstance(data, np.ndarray) else da
   1445 if tindex == 0:
-> 1446     data = lib.concatenate([data_to_concat, data[tindex+1:, :]], axis=0)
   1447 elif tindex == 1:
   1448     data = lib.concatenate([data[:tindex, :], data_to_concat], axis=0)

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 3, the array at index 0 has size 839 and the array at index 1 has size 838

When I print the dimensions of each variable set, I get this result:

ds_u = xr.open_dataset(ufiles[0])
ds_v = xr.open_dataset(vfiles[0])
ds_w = xr.open_dataset(wfiles[0])
mask = xr.open_dataset(mesh_mask)

print(ds_u.dims)
print(ds_v.dims)
print(ds_w.dims)
print(mask.dims)

FrozenMappingWarningOnValuesAccess({'y': 313, 'x': 839, 'deptht': 30, 'time_counter': 31})
FrozenMappingWarningOnValuesAccess({'y': 312, 'x': 838, 'deptht': 30, 'time_counter': 31})
FrozenMappingWarningOnValuesAccess({'y': 313, 'x': 838, 'depthw': 30, 'time_counter': 31})
FrozenMappingWarningOnValuesAccess({'y': 313, 'x': 838, 'z': 30, 't': 1})

Does anyone have any idea what the problem might be? The data is daily from 1993 to 1994, and each file represents one month of each year, for example: GLORYS-U-1993-01, GLORYS-U-1993-02... same for the V and W components.

By the way, this is how I am creating my fieldset:

data_path = 'C:/Users/55859/Desktop/curso_ieapm/velocity_data/' 
mesh_mask = data_path + 'mask.nc'

ufiles = sorted(glob(data_path+'GLORYS-U-*.nc'))
vfiles = sorted(glob(data_path+'GLORYS-V-*.nc'))
wfiles = sorted(glob(data_path+'GLORYS-W-*.nc'))

filenames = {'U': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles[0],'data': ufiles},
             'V': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles[0],'data': vfiles},
             'W': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles[0],'data': wfiles}}

variables = {'U': 'vozocrtx', 'V': 'vomecrty', 'W':'vovecrtz'}
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif','depth': 'depthw', 'time': 'time_counter'},
              'V': {'lon': 'glamf', 'lat': 'gphif','depth': 'depthw', 'time': 'time_counter'},
              'W': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'}}

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

[erikvansebille]

Thanks for reporting, @imariaguiar! This error occurs because the grid for the u field is one larger than for the v field (312 vs 312 and 839 vs 838 grid points). Do you know why that is?

[imariaguiar]

I thought this difference in dimension was linked to the curvilinear grid, because I downloaded and cropped the data to the same latitudes and longitudes. Even the dimension of the W component is different, and I don't understand why. But thank you very much for the response, I'll check into the reason for this.