Incorporation of sgrid conventions

[hetland]

I would like to propose that parcels accept Arakawa C-grids that follow the sgrid conventions. Currently, curvilinear grids in parcels follow a convention that is used by NEMO. Unfortunately, this is not a common convention in other C-grid based models (ROMS and MITGCM, in particular). This would allow a wider variety of models to be easily used, without cumbersome pre-processing.

[erikvansebille]

Hi @hetland, thanks for the question. In fact, from the latest v2.2.1, Parcels supports MITgcm grid indexing, using FieldSet.from_mitgcm(), or by setting Field.gridindexingtype='mitgcm', see also #897.

Is that what you were after?

[hetland]

I did not know about the support for MITgcm indexing. This will be very useful.

This specifically is the part that I was after, but generally I still think considering s-grid support would make it so that a wider range of models could be supported. S-grid conventions would let parcels identify the indexing convention from the Dataset attributes. It would also allow for any type of model output. Basically, the node and face configuration can be specified per dimension, for so that arbitrary (though still structured) grids can be used. This removes the requirement to specify the particular grid configuration for each specific model output type.

[erikvansebille]

Interesting idea, @hetland! But how would you identify the indexing convention from the Dataset attributes? Is there some CF-convention for indexing that we can query from the attributes? Could you give an example?

[hetland]

The S-grid conventions, linked above, describe how this would work: https://sgrid.github.io/sgrid/

A simple example on that page looks like this:

dimensions:
    time = UNLIMITED ;
    inode = 10 ;
    jnode = 20 ;
    icell = 9 ;
    jcell = 19 ;

variables:
    char time(time) ;
        time:standard_name = "time" ;
        time:long_name = "time" ;
        time:units = "seconds since 2015-01-01 00:00:00" ;
    float u(time, jcell, inode) ;
        u:description = "x-velocity" ;
        u:units = "m s-1" ;
        u:grid = "MyGrid" ;
        u:location = "edge1" ;
    float v(time, jnode, icell) ;
        u:description = "y-velocity" ;
        u:units = "m s-1" ;
        u:grid = "MyGrid" ;
        u:location = "edge2" ;
    float c(time, jcell, icell) ;
        c:description = "some concentration" ;
        c:grid = "MyGrid" ;
        c:location = "face" ;
    float node_lat(jnode, inode) ;
        node_lat:standard_name = "latitude" ;
        node_lat:units = "degree_north" ;
    float node_lon(jnode, inode) ;
        node_lon:standard_name = "longitude" ;
        node_lon:units = "degree_east" ;

    int MyGrid ;
        grid:cf_role = grid_topology
        grid:topology_dimension = 2 ;
        grid:node_dimensions = "inode jnode" ;
        grid:face_dimensions = "icell: inode (padding: none) jcell: jnode (padding: none)" ;
        grid:node_coordinates = "node_lon node_lat" ;

[hetland]

So, I've finally been able to get back to this work, and have run into this issue:

NotImplementedError: Curvilinear Grids are not implemented for mitgcm-style grid indexing.
If you have a use-case for this, please let us know by filing an Issue on github

I have a use case for this, and so I'm bringing this point up for discussion. I think one solution would be my original suggestion for generalized support for structured grids, the sgrid conventions.

In the meantime, does anybody have any suggestions for workarounds, preferably that work with xarray datasets (read in from a zarr file)?

[hetland]

I tried to comment this code out, and I can get it to work, but am not convinced that it is working properly. I would be very happy to get some suggestions for a hack. I can provide code for what I'm doing, but thought it might be good to create a very simple test case, based on an analytically based xarray dataset.

Here is some code that might be the foundation for a nice test case. A good (simple) analytic velocity would be u (in the C-grid reference frame, not eastward) velocity could be a constant, v (ditto) could be zero, and the particles should follow the arc of the grid.

The grid object from pygridgen makes it straightforward to get all the components of the grid, like the position on the u, v, and psi points. Let me know if there is interest in using this grid to create a geographically based xarray dataset with time and velocity included.

import numpy as np
import matplotlib.pyplot as plt
from pygridgen import grid

R1 = 10.e5       # Inner radius of arc [m]
R2 = 12.e5        # Outer radius of arc [m]
theta = 60.0     # Sweep of arc [deg]
nR = 20          # number of grid points in radial direction
nT = 100         # number of grid points in theta direction

phi, r = np.meshgrid( np.linspace(-theta/2., theta/2, nT) * np.pi / 180.0, 
                      np.linspace(R1, R2, nR) )

x = r * np.sin(phi)
y = r * np.cos(phi)

x -= x.mean()
y -= y.mean()

grd = grid.CGrid(x, y)

fig = plt.figure(figsize=(15, 5))
ax = fig.add_subplot(111)
plt.plot(grd.x_vert/1e3, grd.y_vert/1e3, '-k');
plt.plot(grd.x_vert.T/1e3, grd.y_vert.T/1e3, '-k');
ax.set_aspect(1.0)

[hetland]

Here is some code to make a very basic dataset. I would like to be able to load this into parcels using FieldSet.from_xarray_dataset. These are based on the grid as defined above.

I have two examples. The first uses the velocities native to the C-Grid, that is u and v have different shapes, and are not co-located. Note in this first example, interp_method='cgrid_velocity' is not possible, since the velocities have different dimensions.

import pandas as pd

eta_rho, xi_u = grd.x_u.shape
eta_v, xi_rho = grd.x_v.shape

periods = 30
time = pd.date_range("2000-01-01", periods=periods)

u = np.ones((periods, eta_rho, xi_u)) * 0.05
v = np.zeros((periods, eta_v, xi_rho))

ds = xr.Dataset(
    data_vars=dict(
        u=(['time', 'eta_rho', 'xi_u'], u),
        v=(['time', 'eta_v', 'xi_rho'], v),
    ),
    coords=dict(
        x_u=(['eta_rho', 'xi_u'], grd.x_u),
        y_u=(['eta_rho', 'xi_u'], grd.y_u),
        x_v=(['eta_v', 'xi_rho'], grd.x_v),
        y_v=(['eta_v', 'xi_rho'], grd.y_v),
        time=("time", time),
    ),
)


variables = {"U": "u", "V": "v"}
dimensions = {
    "U": {"lon": "x_u", "lat": "y_u", "time": "time"},
    "V": {"lon": "x_v", "lat": "y_v", "time": "time"},
}

fieldset = FieldSet.from_xarray_dataset(ds, variables, dimensions, 
                                        mesh='flat',
                                        gridindexingtype='mitgcm')


pset = ParticleSet(fieldset=fieldset, 
                   pclass=JITParticle,
                   lon=grd.x_rho[:, 5],  
                   lat=grd.y_rho[:, 5])   

output_file = pset.ParticleFile(name="foo.nc", outputdt=timedelta(hours=1))

pset.execute(AdvectionRK4,                 # the kernel (which defines how particles move)
             runtime=timedelta(days=15),    # the total length of the run
             dt=timedelta(minutes=5),      # the timestep of the kernel
             output_file=output_file,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}
            )

pset.show(field=fieldset.U)

This case fails, as the constant u velocity is interpreted as a constant eastward velocity, not as a constant velocity aligned with the grid.

The second case uses co-located velocity points, which succeeds at least in advecting the parcels along the grid, not simply to the east:

import pandas as pd

eta_psi, xi_psi = grd.x_psi.shape

periods = 30
time = pd.date_range("2000-01-01", periods=periods)

u = np.ones((periods, eta_psi, xi_psi)) * 0.05
v = np.zeros((periods, eta_psi, xi_psi))

ds = xr.Dataset(
    data_vars=dict(
        u=(['time', 'eta_psi', 'xi_psi'], u),
        v=(['time', 'eta_psi', 'xi_psi'], v),
    ),
    coords=dict(
        x_psi=(['eta_psi', 'xi_psi'], grd.x_psi),
        y_psi=(['eta_psi', 'xi_psi'], grd.y_psi),
        time=("time", time),
    ),
)

variables = {"U": "u", "V": "v"}
dimensions = {
    "U": {"lon": "x_psi", "lat": "y_psi", "time": "time"},
    "V": {"lon": "x_psi", "lat": "y_psi", "time": "time"},
}

fieldset = FieldSet.from_xarray_dataset(ds, variables, dimensions, 
                                        mesh='flat',
                                        interp_method='cgrid_velocity',
                                        gridindexingtype='mitgcm')


pset = ParticleSet(fieldset=fieldset, 
                   pclass=JITParticle,
                   lon=grd.x_rho[:, 5],  
                   lat=grd.y_rho[:, 5])   

output_file = pset.ParticleFile(name="foo.nc", outputdt=timedelta(hours=1))

pset.execute(AdvectionRK4,                 # the kernel (which defines how particles move)
             runtime=timedelta(days=15),    # the total length of the run
             dt=timedelta(minutes=5),      # the timestep of the kernel
             output_file=output_file,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}
            )

pset.show(field=fieldset.U)

Though the second method seems to work, the part I still don't understand is how interp_method='cgrid_velocity' works if the points are co-located? This seems to contradict the documentation on c-grids. Are the points secretly staggered, or is it assumed they are interpolated to psi-points (the interior grid vertices) and thus co-located?

[erikvansebille]

Thanks @hetland, for providing this extensive code. It's in idea similar to how we test NEMO in this notebook. I have a lot of teaching prep to do tomorrow, but may have some time later in the week to look at this. Stay tuned!

[hetland]

Thanks for taking the time to look at this, @erikvansebille.

I've got a hack that works now. I needed to trim my velocity fields such that the grids are the same shape as the f points (the psi points in my ROMS-centric notation above). This is actually a bit of a pain in xarray, and there is no straightforward way to do this currently. My current dataset is small enough that I just read in the values as a numpy array, and define new f-point (or psi-point) dimensions; generally a poor solution.

I believe that the problem I am having is based on the fact that the velocities currently must have the same shape as the f-grid. However, I don't think this is common, and is certainly not a natural way to define velocities on a Cgrid. The relative positions is the focus of the differences between NEMO and MITGCM indexing, but in both cases, Parcels requires the shape of the velocity field must be identical. In ROMS (and I believe in MITGCM) the u and v fields have a different shape, as you can see in this figure of a whole ROMS domain. I think this is a natural way for velocities to be defined, otherwise, e.g., the eastern and western boundaries are not the same relative to a tracer field.

So, to bring this back to the s-grid conventions, being able to define a CGrid velocity field where u and v have different shapes will allow less ad hoc modification of the raw model output to input into parcels.

[daanreijnders]

Hi @hetland,

I'm also using Parcels in combination with MITgcm and was involved with adding the support over MITgcm's indexing within Parcels, so I think I can shed some light on this.

Of the C-grid models, NEMO was the first to be supported by Parcels. While in C-grids the U and V velocities are not co-located, NEMO stores the U and V outputs in arrays of the same shape, with a row (/column) of zeros prepended so that U has the same shape as V (/V as U). C-grid support in Parcels was then implemented to play nice with this specific type of data shape and indexing. Parcels needs to be fed the U and V data (same shape) and the f-point coordinates. Parcels doesn't assume U and V are collocated, but works out the correct positioning of U and V through the f-point coordinates.

Then there's MITgcm. The shape of MITgcm output depends on the method you store your files with. If you use the mnc module, MITgcm outputs netCDF instead of its default bespoke 'mds' binary files. In mnc output, U and V have unequal sizes: the amount of U and V points in the output corresponds to the actual amount of U and V points in the C-grid stencil. Another way to generate MITgcm output is using xmitgcm, which reads the 'mds' files into xarray, allowing netCDF export in turn. xmitgcm stores U and V in equal shapes, but instead of of adding a row/column of zeros, it trims the data by a row/column. This way of storing data doesn't actually lose information though: since MITgcm natively assumes periodic boundaries, the first column of U is the same as the last (similar for V; closed domains are constructed by setting cells at the boundary to zero depth).

When I started to work with MITgcm, there were two issues: I used mnc, so U and V were of different shape (not supported by Parcels) and the NEMO-like indexing assumption didn't work with MITgcm data, due to the row/column of zeros being prepended in NEMO for U/V. We fixed the second issue within Parcels itself, so that MITgcm-style indexing was allowed (see #897). I solved the first issue for myself by appending a row/column of zeros to my netCDF files through a simple script, but I eventually started using xmitgcm to circumvent this issue altogether (xmitgcm also came with other benefits like the xarray interface and strict CF-compliance).

So: the way Parcels handles C-grid data is historically motivated ('NEMO first') and works fine with MITgcm if xmitgcm is used for storing output. Changing Parcels so that it takes in arbitrary C-grid data should be possible, of course. I think the first step towards that would be to drop the assumption that U and V have the same size and inferring the locations of the f-points from the coordinates of U and V that are readily available. The type of indexing (mitgcm versus nemo) can then also be generalized in accordance to sgrid-conventions. Parsing the type of indexing from netCDF also shouldn't be that hard, I believe. But altogether these steps are still quite the coding project.