Loading data from GEFS (FV3) tiled restart files with 2d-coords

[StevePny]

Hi,

I'd like to read the 6-tiled restart files coming from the NOAA GEFS ensemble, which are output from the FMS-based FV3 cubed sphere atmosphere model, for example ensemble member 1 (out of 30) from 2/8/24 is here:
https://noaa-gefs-pds.s3.amazonaws.com/index.html#gefs.20240208/00/atmos/init/p01/

There are 6 files for the atmosphere and 6 files for the surface:

Typically, I believe this would be done with the GFDL command-line utility mppnccombine (from https://github.com/NOAA-GFDL/FRE-NCtools), however that has not worked for me (I think possibly because I'm using a mac M1, or that I'm building FRE-NCtools in docker).

I'd like to be able to read the files with Python using xarray and then process the data as a single Dataset object. Ideally this would be via a function like xarray.open_mfdataset, so it still has lazy access to all data files seamlessly.

The netcdf data header looks like this (in this case I merged the 30 ensemble members for each tile into one netcdf file):

Each variable has a 2D coordinate (geolon*, geolat*), each of which is defined by different dimension variables (lon,lat,lonp,latp,lev,levp).

It appears that upon reading in, xarray can identify from the netcdf file that the geolon*/geolat* variables are the coordinates:

And, I am able to plot the global map of the data using the cartopy projection and matplotlib transform option. Here I provide "geolon" and "geolat" to the plotting routine so that it knows where to put the values in the final map plot. (In the case below, since I'm reading in the entire ensemble for each tile I'm plotting the mean here. Also, here I renamed lon,lat,lev,record to longitude,latitude,level,number):

The result looks like this:

Just for reference, without fixing the colorbars you can see the tile locations more clearly, which shows where the tiles are separated:

However, again, I would like to perform further operations on the dataset, which include subsetting local patches that may overlap tile boundaries. I have not yet found anywhere where combining such tile files has been done with xarray, and the built in tools that I've tried to use (combine_nested, combine_by_coords, combine_first, open_mfdataset, concat, merge) have not worked. While I'm sure there is a very manual way to reconstruct a new Dataset object using these inputs, I'd like to avoid reinventing anything that might already work. Is there a way that I'm missing, using built-in xarray tools, to merge these 6 files either into a single Dataset object, or to merge them into a new file (similarly to mppnccombine ), so that in either case it can take advantage of the xarray lazy loading while working with the full dataset?

[spencerkclark]

I think it depends on what you mean by "combine." The topology of the cubed sphere is such that you cannot arrange the tiles cleanly on a logically rectangular grid. The best you can do is something like what the to_cross function defined here does, but that obviously fills much of the rectangle with NaNs. See for example the plot below:

Instead I typically find simply concatenating the tiles along a new "tile" dimension is sufficient for most of the operations I am interested in (this is something that can be done out of the box with concat, combine_nested, or open_mfdataset). While it does not completely solve plotting without another helper function, you can still compute regional averages that cross tiles (masking by latitude and longitude) fairly straightforwardly.

Note that mppnccombine has a different purpose. It is meant for combining I/O subdomains into individual tiles (which are cleanly logically rectangular) rather than combining the full tiles themselves into some other shape like to_cross. It is more like a special case of combine_by_coords.

[StevePny]

I think it depends on what you mean by "combine."

I suppose what I meant here was to have a single Dataset that allows one to access the field values with 'lat' and 'lon' coordinates, despite the fact that the multi-tile grid is not logically rectangular and regardless of the underlying data structure, similar to the transform functionality from cartopy/matplotlib mentioned above.

I tried adding a tile dimension to each dataset and running:

ds = xr.combine_nested(datasets=[ds1,ds2,ds3,ds4,ds5,ds6],concat_dim='tile')

With this approach I run out of memory (as I mentioned above I'm working with the full 30-member ensemble). I assume a copy is being done.

This is closer to what I was looking for:

ds = xr.open_mfdataset(paths="gfs_data.tile?.ens.nc", combine='nested', concat_dim='tile')

It still doesn't get me the functionality I'd like though, so I may end up having to just copy the fields into a separate array for now.

[TomNicholas]

Xgcm has some limited support for some of the functionality you are after, because it was made to handle the lat-lon-cap geometry of MITgcm. You could use the face connections spec to tell xgcm how your different rectangular subarrays are connected (your tile dimension is the face dimension in xGCM-speak).

It wouldn't replace your open_mfdataset call, but you could use this with dask to lazily combine your dataset together, organized via an accompanying xgcm.Grid object.

@jbusecke how far did we actually get with supporting operations on face-connected grids? Did we implemented it for the pad method? I can't remember...

But honestly the general version of this is a massive rabbit hole, and I'm not aware of a package that supports applying operations across these boundaries neatly out of the box.