Ocean grid compatibility with the data manager

[jkrasting]

The tropical sea level POD submitted by @chiaweh2 is running into a problem with the new data manager. The error is:

ERROR: cf_xarray error: couldn't assign ['nlat', 'nlon'] to axes for tauvo(assigned axes: {'X': ['lon'], 'Y': ['lat'], 'Z': [], 'T': ['time']})

This seems related to the 2D lat/lon coordinate variables which are defined by the 1D nlat/nlon index dimensions. This is a very common data format in the ocean model realm.

@tsjackson-noaa - should the new data manager be able to handle this?

netcdf CESM2_omip1_r1i1p1f1_gn.tauuo.mon {

dimensions:
        time = UNLIMITED ; // (3720 currently)
        nlat = 384 ;
        nlon = 320 ;
        d2 = 2 ;
        vertices = 4 ;

variables:
        float tauuo(time, nlat, nlon) ;
                tauuo:_FillValue = 1.e+20f ;
                tauuo:cell_measures = "--OPT" ;
                tauuo:cell_methods = "time: mean" ;
                tauuo:comment = "TAUX" ;
                tauuo:coordinates = "time lat lon" ;
                tauuo:description = "This is the stress on the liquid ocean from overlying atmosphere, sea ice, ice shelf, etc." ;
                tauuo:frequency = "mon" ;
                tauuo:id = "tauuo" ;
                tauuo:long_name = "Surface Downward X Stress" ;
                tauuo:mipTable = "Omon" ;
                tauuo:missing_value = 1.e+20 ;
                tauuo:out_name = "tauuo" ;
                tauuo:positive = "down" ;
                tauuo:prov = "Omon ((isd.003))" ;
                tauuo:realm = "ocean" ;
                tauuo:standard_name = "surface_downward_x_stress" ;
                tauuo:time = "time" ;
                tauuo:time_label = "time-mean" ;
                tauuo:time_title = "Temporal mean" ;
                tauuo:title = "Surface Downward X Stress" ;
                tauuo:type = "real" ;
                tauuo:units = "N m-2" ;
                tauuo:variable_id = "tauuo" ;

        double lat(nlat, nlon) ;
                lat:axis = "Y" ;
                lat:bounds = "lat_bnds" ;
                lat:standard_name = "latitude" ;
                lat:title = "Latitude" ;
                lat:type = "double" ;
                lat:units = "degrees_north" ;
                lat:valid_max = 90. ;
                lat:valid_min = -90. ;

        double lon(nlat, nlon) ;
                lon:axis = "X" ;
                lon:bounds = "lon_bnds" ;
                lon:standard_name = "longitude" ;
                lon:title = "Longitude" ;
                lon:type = "double" ;
                lon:units = "degrees_east" ;
                lon:valid_max = 360. ;
                lon:valid_min = 0. ;

        int nlat(nlat) ;
                nlat:long_name = "cell index along second dimension" ;
                nlat:units = "1" ;
        int nlon(nlon) ;
                nlon:long_name = "cell index along first dimension" ;
                nlon:units = "1" ;


        float lat_bnds(nlat, nlon, vertices) ;
                lat_bnds:units = "degrees_north" ;
        float lon_bnds(nlat, nlon, vertices) ;
                lon_bnds:units = "degrees_east" ;

        double time(time) ;
                time:axis = "T" ;
                time:bounds = "time_bnds" ;
                time:standard_name = "time" ;
                time:title = "time" ;
                time:type = "double" ;
                time:units = "days since 0001-01-01 00:00:00" ;
                time:calendar = "365_day" ;
        double time_bnds(time, d2) ;
                time_bnds:calendar = "noleap" ;
                time_bnds:units = "days since 0001-01-01 00:00:00" ;

[tsjackson-noaa]

Hi @jkrasting,
The short answer to your question is "no." This coordinate configuration is covered by the CF conventions but doesn't appear to be supported by cf_xarray -- I'm inferring this from the absence of language covering this case in their docs, and the fact that the framework is passing through the error raised by cf_xarray. My broad understanding is that this aspect of cf_xarray was based on MetPy's functionality, which might make it less surprising that they haven't thought of oceanic applications.

At least at this point, it appears that we'll need to roll our own solution, as I'm not aware of any project that provides (better) CF convention support within the xarray ecosystem.

[jkrasting]

Thanks for the assessment @tsjackson-noaa. Cecilia's sea ice PODs will probably run into the same issue since the sea ice model shares the same grid as the ocean. We'll have to strategize something for this round of PODs.

[mgrover1]

There is a conversation ongoing within the cf-xarray repo xarray-contrib/cf-xarray#204 related to this... might be helpful to collaborate on this.

[dcherian]

Can you post a minimal example illustrating the problem with cf_xarray? It's not clear what you're trying to do.

doesn't appear to be supported by cf_xarray

I'm not sure what this means.

EDIT: Ignore the rest, I misread the repr.

This kind of dataset has the issue that nlon,nlat don't have anything except a size associated with them. There is no underlying variable, so there are no attributes, so there's nothing for cf_xarray to interpret. We discussed this for the ROMS ocean model here (xarray-contrib/cf-xarray#84) where xi_rho ~ nlon and eta_rho ~ nlat.

In practice, I've assigned a ds["nlon"] = ("nlon", np.arange(ds.sizes("nlon")), {"axis": "X"}). This lets me do ds.cf.mean("X") ~ ds.mean("nlon")

There is a similar issue for models with a parametric vertical coordinate (xarray-contrib/cf-xarray#87) where s_rho ~ nlon-like though for that one we could merge xarray-contrib/cf-xarray#155

As @mgrover1 points out, the only solution I can see is to support sgrid which explicitly labels the nlon, nlat dimensions and the lat, lon coordinates in a standard way.

Very curious to hear if you have an alternate solution.

[dcherian]

Sorry, totally mistaken. You do have an nlon, nlat variable! (unlike most model datasets I work with)

        int nlat(nlat) ;
                nlat:long_name = "cell index along second dimension" ;
                nlat:units = "1" ;
        int nlon(nlon) ;
                nlon:long_name = "cell index along first dimension" ;
                nlon:units = "1" ;

So now I don't know what you're trying to do :). Is "cell index along first dimension" a standard CMOR/OMIP thing? If so we could add it to our autoguessing criteria for axis="X" etc. Would that help?

[wrongkindofdoctor]

@chiaweh2 Do you have any insight re: @dcherian's question about the nlat and nlon descriptions?

[chiaweh2]

@wrongkindofdoctor Yes I think the description @dcherian suggested sounds reasonable. I have looked at CESM2 and GFDL-OM4 output in CMIP6-OMIP so far. I don't think they have a consistent long_name.
CESM2 uses "nlon", "nlat".
Attributes:
long_name: cell index along first dimension
units: 1

GFLD-OM4 uses "x" and "y" for the cell index:
Attributes:
long_name: x coordinate of projection
units: degrees
cartesian_axis: X
standard_name: projection_x_coordinate
axis: X

[jkrasting]

Yes, I also agree this could help.--