Documentation on dimensions and indexing in Parcels

Parcels supports Fields on any curvilinear Grid. For velocity Fields (U, V and W), there are some additional restrictions if the Grid is discretized as an Arakawa B- or Arakawa C-grid. That is because under the hood, Parcels uses a specific interpolation scheme for each of these grid types. This is described in Section 2 of Delandmeter and Van Sebille (2019) and summarized below.

The summary is: for Arakawa B-and C-grids, Parcels requires the locations of the corner points ($f$-points) of the grid cells for the dimensions dictionary of velocity Fields. In other words, on an Arakawa C-grid, the [k, j, i] node of U will not be located at position [k, j, i] of U.grid.

Barycentric coordinates and indexing in Parcels

arakawa C-grids

In a regular grid (also called an Arakawa A-grid), the velocities (U, V and W) and tracers (e.g. temperature) are all located in the center of each cell. But hydrodynamic model data is often supplied on staggered grids (e.g. for NEMO, POP and MITgcm), where U, V and W are shifted with respect to the cell centers.

To keep it simple, let's take the case of a 2D Arakawa C-grid. Zonal (U) velocities are located at the east and west faces of each cell and meridional (V) velocities at the north and south. The following diagram shows a comparison of 4x4 A- and C-grids.

Arakawa Grid layouts

Note that different models use different indices to relate the location of the staggered velocities to the cell centers. The default C-grid indexing notation in Parcels is the same as for NEMO (middle panel): U[j, i] is located between the cell corners [j-1, i] and [j, i], and V[j, i] is located between cell corners [j, i-1] and [j, i].

Now, as you've noticed on the grid illustrated on the figure, there are 4x4 cells. The grid providing the cell corners is a 5x5 grid, but there are only 4x5 U nodes and 5x4 V nodes, since the grids are staggered. This implies that the first row of U data and the first column of V data is never used (and do not physically exist), but the U and V fields are correctly provided on a 5x5 table. If your original data are provided on a 4x5 U grid and a 5x4 V grid, you need to regrid your table to follow Parcels notation before creating a FieldSet!

MITgcm uses a different indexing: U[j, i] is located between the cell corners [j, i] and [j+1, i], and V[j, i] is located between cell corners [j, i] and [j, i+1]. If you use xmitgcm to write your data to a NetCDF file, U and V will have the same dimensions as your grid (in the above case 4x4 rather than 5x5 as in NEMO), meaning that the last column of U and the last row of V are omitted. In MITgcm, these either correspond to a solid boundary, or in the case of a periodic domain, to the first column and row respectively. In the latter case, and assuming your model domain is periodic, you can use the FieldSet.add_periodic_halo method to make sure particles can be correctly interpolated in the last zonal/meridional cells.

Parcels can take care of loading C-grid data for you, through the general FieldSet.from_c_grid_dataset method (which takes a gridindexingtype argument), and the model-specific methods FieldSet.from_nemo and FieldSet.from_mitgcm.

The only information that Parcels needs is the location of the 4 corners of the cell, which are called the $f$-node. Those are the ones provided by U.grid (and are equal to the ones in V.grid). Parcels does not need the exact location of the U and V nodes, because in C-grids, U and V nodes are always located in the same position relative to the $f$-node.

Importantly, the interpolation of velocities on Arakawa C-grids is done in barycentric coordinates: $\xi$, $\eta$ and $\zeta$ instead of longitude, latitude and depth. These coordinates are always between 0 and 1, measured from the corner of the grid cell. This is more accurate than simple linear interpolation on curvilinear grids.

When calling FieldSet.from_c_grid_dataset() or a method that is based on it like FieldSet.from_nemo(), Parcels automatically sets the interpolation method for the U, V and W Fields to cgrid_velocity. With this interpolation method, the velocity interpolation follows the description in Section 2.1.2 of Delandmeter and Van Sebille (2019).

NEMO Example
In [1]:
from parcels import FieldSet
from glob import glob
import numpy as np
from os import path
from datetime import timedelta as delta

Let's see how this works. We'll use the NemoNorthSeaORCA025-N006_data data, which is on an arakawa C-grid. If we create the FieldSet using the coordinates in the U, V and W files, we get an Error as seen below:

In [2]:
data_path = 'NemoNorthSeaORCA025-N006_data/'
ufiles = sorted(glob(data_path+'ORCA*U.nc'))
vfiles = sorted(glob(data_path+'ORCA*V.nc'))
wfiles = sorted(glob(data_path+'ORCA*W.nc'))

filenames = {'U': ufiles, 'V': vfiles, 'W': wfiles}
variables = {'U': 'uo', 'V': 'vo', 'W': 'wo'}
dimensions = {'U': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthu', 'time': 'time_counter'},
              'V': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthv', 'time': 'time_counter'},
              'W': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthw', 'time': 'time_counter'}}

fieldset = FieldSet.from_nemo(filenames, variables, dimensions)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-551c32df110d> in <module>
     10               'W': {'lon': 'nav_lon', 'lat': 'nav_lat', 'depth': 'depthw', 'time': 'time_counter'}}
     11 
---> 12 fieldset = FieldSet.from_nemo(filenames, variables, dimensions)

~/Codes/PARCELScode/parcels/fieldset.py in from_nemo(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, field_chunksize, **kwargs)
    463         fieldset = cls.from_c_grid_dataset(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
    464                                            allow_time_extrapolation=allow_time_extrapolation, tracer_interp_method=tracer_interp_method,
--> 465                                            field_chunksize=field_chunksize, **kwargs)
    466         if hasattr(fieldset, 'W'):
    467             fieldset.W.set_scaling_factor(-1.)

~/Codes/PARCELScode/parcels/fieldset.py in from_c_grid_dataset(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, field_chunksize, **kwargs)
    524 
    525         if 'U' in dimensions and 'V' in dimensions and dimensions['U'] != dimensions['V']:
--> 526             raise ValueError("On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. "
    527                              "See also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb")
    528         if 'U' in dimensions and 'W' in dimensions and dimensions['U'] != dimensions['W']:

ValueError: On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U and V. See also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb

We can still load the data this way, if we use the FieldSet.from_netcdf() method. But because this method assumes the data is on an Arakawa-A grid, this will mean wrong interpolation.

In [3]:
fieldsetA = FieldSet.from_netcdf(filenames, variables, dimensions)

Instead, we need to provide the coordinates of the $f$-points. In NEMO, these are called glamf, gphif and depthw (in MITgcm, these would be called XG, YG and Zl). The first two are in the coordinates.nc file, the last one is in the W file.

In [4]:
mesh_mask = data_path + 'coordinates.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}}

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'}}

fieldsetC = FieldSet.from_nemo(filenames, variables, dimensions)
WARNING: File NemoNorthSeaORCA025-N006_data/coordinates.nc could not be decoded properly by xarray (version 0.15.1).
         It will be opened with no decoding. Filling values might be wrongly parsed.
WARNING: Casting lon data to np.float32
WARNING: Casting lat data to np.float32
WARNING: Trying to initialize a shared grid with different chunking sizes - action prohibited. Replacing requested field_chunksize with grid's master chunksize.

We can plot the different grid points to see that indeed, the longitude and latitude of fieldsetA.U and fieldsetA.V are different.

In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa

fig, ax = plt.subplots()
nind = 3
ax1 = ax.plot(fieldsetA.U.grid.lon[:nind, :nind], fieldsetA.U.grid.lat[:nind, :nind], '.r', label='U points assuming A-grid')
ax2 = ax.plot(fieldsetA.V.grid.lon[:nind, :nind], fieldsetA.V.grid.lat[:nind, :nind], '.b', label='V points assuming A-grid')

ax3 = ax.plot(fieldsetC.U.grid.lon[:nind, :nind], fieldsetC.U.grid.lat[:nind, :nind], '.k', label='C-grid points')

ax.legend(handles=[ax1[0], ax2[0], ax3[0]], loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()

Further information about the NEMO C-grids is available in the NEMO 3D tutorial.

Arakawa B-grids

Interpolation for Arakawa B-grids is detailed in Section 2.1.3 of Delandmeter and Van Sebille (2019). Again, the dimensions that need to be provided are those of the barycentric cell edges (i, j, k).

To load B-grid data, you can use the method FieldSet.from_b_grid_dataset, or specifically in the case of POP-model data FieldSet.from_pop.

3D C- and B-grids

For 3D C-grids and B-grids, the idea is the same. It is important to follow the indexing notation, which is defined in Parcels and in Delandmeter and Van Sebille (2019). In 3D C-grids, the vertical (W) velocities are at the top and bottom. Note that in the vertical, the bottom velocity is often omitted, since a no-flux boundary conditions implies a zero vertical velocity at the ocean floor. That means that the vertical dimension of W often corresponds to the amount of vertical levels.