Storing the fieldset components (U, V, W, S, T) as netcdfs

[ischiller-weiss]

Question

Hello!

I am using C-gridded ocean model output and I am hoping to save the FieldsetC (U, V, W, S, T) components as individual netcdf files so one doesn't need to load the FieldsetC every time when running a parcels experiment.

I think this is quite a basic thing, but when writing the FieldsetC.write, there is an error that's thrown. My question is whether, when loading the entire Fieldset, can one store either the whole Fieldset or Fieldset.U, V, W, etc. as individual files and use them later?

Thank you so much! Please let me know if my question or anything is unclear as well.

Supporting code/error messages

ufiles = sorted(glob(f"{inpath}/U/*.nc"))[2:max_ind]
#ufiles = [f for f in ufiles if '2019' in f]

vfiles = sorted(glob(f"{inpath}/V/*.nc"))[2:max_ind] #[-366*3:] #[:5]#[:360]
#vfiles = [f for f in vfiles if '2019' in f]

wfiles = sorted(glob(f"{inpath}/W/*.nc"))[2:max_ind] #[-366*3:] #[:5]#[:360] #-366*3:
#wfiles = [f for f in wfiles if '2019' in f]

sfiles = sorted(glob(f"{inpath}/S/*.nc"))[2:max_ind] #[-366*3:]#[:5]
#sfiles = [f for f in sfiles if '2019' in f]

tfiles = sorted(glob(f"{inpath}/T/*.nc"))[2:max_ind] #[-366*3:]#[:5] 
#tfiles = [f for f in tfiles if '2019' in f]

variables = {'U': 'vozocrtx', 'V': 'vomecrty', 'W': 'vovecrtz', 'S': 'vosaline', 'T': 'votemper'}

filenames = {
    "U": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": ufiles}, # must use same everywhere w,files. but w depth is 0, northeast corner of T grid is lid
    "V": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": vfiles},
    "W": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": wfiles},
    "S": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": sfiles},
    "T": {"lon": wfiles[0], "lat": wfiles[0], "depth": wfiles[0], "data": tfiles}
}


c_grid_dimensions = {
    "lon": "nav_lon",
    "lat": "nav_lat",
    "depth": "depthw",
    "time": "time_counter",
}

dimensions = {
    "U": c_grid_dimensions,
    "V": c_grid_dimensions,
    "W": c_grid_dimensions,
    "S": c_grid_dimensions,
    "T": c_grid_dimensions
}

with warnings.catch_warnings():
    warnings.simplefilter("ignore", parcels.FileWarning)
    fieldsetC = parcels.FieldSet.from_nemo(filenames, variables, dimensions,
        allow_time_extrapolation=True)

# as an example: 

fieldsetC.U.write('fieldsetC_U_test.nc')

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[53], line 1
----> 1 fieldsetC.U.write('fieldsetC_U_test.nc')

File ~/miniconda3/envs/py3_std_maps_2023-11-20/lib/python3.12/site-packages/parcels/field.py:1787, in Field.write(self, filename, varname)
   1784 attrs = {"units": "seconds since " + str(self.grid.time_origin)} if self.grid.time_origin.calendar else {}
   1785 time_counter = xr.DataArray(self.grid.time, dims=["time_counter"], attrs=attrs)
   1786 vardata = xr.DataArray(
-> 1787     self.data.reshape((self.grid.tdim, self.grid.zdim, self.grid.ydim, self.grid.xdim)),
   1788     dims=["time_counter", vname_depth, "y", "x"],
   1789 )
   1790 # Create xarray Dataset and output to netCDF format
   1791 attrs = {"parcels_mesh": self.grid.mesh}

AttributeError: 'DeferredArray' object has no attribute 'reshape'

[erikvansebille]

Hi @ischiller-weiss, thanks for using Parcels. However, I'm not sure why you would want to write your Field. The Fieldwrite() method is normally only used in particle-field interaction, like in this tutorial. Otherwise, I've never seen a need to write your hydrodynamic file.

If you want to separate out the different variables of your model files, I suggest using the nco package or xarray in Python. I wouldn't use Parcels for that...

[ischiller-weiss]

Ok, thank you! Hm I think the issue for me is that the Fieldset takes about 1 hour to load all the gridded fields prior to running a parcels expt, so I wasn't sure if there was a way to store the Fieldset in some way since it's just based on the input data...

[erikvansebille]

Ah, we've previously experienced that the opening of all files to read in the time-information can take very long, when constructing a FieldSet with many files.

If that's indeed the bottleneck in your case too, then a workaround is to explicitly provide a lumpy array of the timestamps in each file. This means that Parcels doesn't have to open them anymore and creating a FieldSet could go much faster.

See below an example code we use for the MOi GLO12 data:

data_path = '/storage/shared/oceanparcels/input_data/MOi/'
ufiles = sorted(glob(data_path+'GLO12/psy4v3r1-daily_U_*.nc'))
vfiles = [f.replace('_U_', '_V_') for f in ufiles]
wfiles = [f.replace('_U_', '_W_') for f in ufiles]
mesh_mask = data_path + 'domain_ORCA0083-N006/PSY4V3R1_mesh_hgr.nc'

# Create list of timestamps for _much_ faster loading of the fieldset
timestamps = np.arange(np.datetime64(ufiles[0][-13:-3]),
                       np.datetime64(ufiles[-1][-13:-3])+np.timedelta64(1, 'h'),  # adding an hour to include the last day
                       np.timedelta64(1, 'D'))
assert len(ufiles) == len(vfiles) == len(timestamps), 'Different number of U, V files and timestamps'
timestamps = [[t] for t in timestamps]  # convert to a list of lists

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}}
variables = {'U': 'vozocrtx', 'V': 'vomecrty'}
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'},
              'V': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'}}
fieldset = FieldSet.from_nemo(filenames, variables, dimensions, timestamps=timestamps)

[ischiller-weiss]

Thank you so much! Wow yes, it's much much faster. I really appreciate all the help!