Adding/augmenting functions to write a CF Coordinate Reference System (CRS) "container variable"

[durack1]

The DRCDP project has anticipated regional contributions (mainly for the CONUS) many of which in their native output include grid mapping criterion. While no obvious and extensive specification exists, some samples are available, including 3 defining a CRS WKT (CRS Well-Known Text format; also known as OGC WKT) within the CF Conventions that define the specifications for the "latitude_longitude", "transverse_mercator" using several common datums. The int crs variable written to a netCDF file has no associated data, but contains variable attributes which define characteristics of the grid mapping, .e.g., see here and example below

  int crs ;
    crs:grid_mapping_name = "latitude_longitude";
    crs:longitude_of_prime_meridian = 0.0 ;
    crs:semi_major_axis = 6378137.0 ;
    crs:inverse_flattening = 298.257223563 ;
    crs:crs_wkt =
     GEODCRS["WGS 84",
     DATUM["World Geodetic System 1984",
       ELLIPSOID["WGS 84",6378137,298.257223563,
         LENGTHUNIT["metre",1.0]]],
     PRIMEM["Greenwich",0],
     CS[ellipsoidal,3],
       AXIS["(lat)",north,ANGLEUNIT["degree",0.0174532925199433]],
       AXIS["(lon)",east,ANGLEUNIT["degree",0.0174532925199433]],
       AXIS["ellipsoidal height (h)",up,LENGTHUNIT["metre",1.0]]]

In some parallel discussions see PCMDI/obs4MIPs-cmor-tables#375, there may also be a need to account for "Equal-Area Scalable Earth" EASE grids (see here). Ongoing discussions about the DRCDP targeted MACAv3 dataset (thanks @khegewisch) have noted the following datums:

1. WGS84 (EPGS: 4326) is the latitude_longitude that CMIP, and Google Earth uses
2. EPSG:3857 a Spherical Mercator projection coordinate system popularized by web services such as Google Maps, Leaflet Maps and OpenStreetMap
3. Lambert Conformal Conic projection is used for WRF outputs (I've had to convert data in this projection to latitude_longitude before to do statistical downscaling on this in a hybrid method)
4. I have seen many more projections in requesting non-climate data using APIs from ESRI, USDA, USGS, EPA, etc, etc. Some of the different agencies have a favorite projection for their data and it doesn't always align with EPSG: 4326 or 3857.. but those are the most common.

CMOR3 has a couple of functions cmor.set_grid_mapping and cmor.grid that appear to provide some functionality, but currently do not support the latitude_longitude example. There is also a CMOR test which provides a polar_stereographic projection example. We'll need to do a little more investigation to more clearly define what projections/grid_mappings need support before acting to CMOR changes to enable an appropriate crs variable to be written.

And just to collapse email links into this issue, a couple of GIS software packages have vague docs pages on coordinate reference systems (the CRS acronym) QGIS, including a demo of how to use NCO tools to add the crs attribute after the fact, thanks StackExchange, also ESRI ArcGIS, GDAL "there is no universal way of storing georeferencing in NetCDF files", and finally rioxarray. It would be useful to ascertain whether a formal standard, other than the example defined in the CF Conventions (link up top).

Also adding the augmented CF to/from WKT/PROJ.4 translation spec - http://cfconventions.org/wkt-proj-4.html

ping @taylor13 @gleckler1 @paullric

[TomLav]

I will be interested in following and possibly contribute to this discussion. My main motivation is that the (observational) data I am responsible for are in polar-focused grids (mainly 6931 and 6932).

The pyresample python library, part of the pytroll ecosystem is another software handling map projections that might help convert to/from CF Conventions. pyresample uses pyproj under the hood.

It is noteworthy that having a crs variable in one's netCDF file is not enough for being CF. One must also correctly refer to it via :grid_mapping attributes, and one must also define coordinate variables using the correct :standard_name attribute (:standard_name = 'projection_x_coordinate' for many but not all projections). This is maybe already well known and covered by existing cmor functionalities, I do not know.

[durack1]

I will be interested in following and possibly contribute to this discussion.

@TomLav great, it would be really useful to have broad coverage on this so we can implement a solution that satisfies CF requirements, and also captures all required grid mapping flavours that CMOR users will need.

[durack1]

Just linking parallel CF discussions - https://github.com/orgs/cf-convention/discussions/404

[durack1]

Collating this comprehensive list of grid mappings, might be easier than I thought - the CF "Appendix H: Annotated Examples of Discrete Geometries" lists 16, both gridded, point data etc.

[TomLav]

@durack1, I know nothing of the CMOR software and if it is cumbersome to make it depend on other python modules.

But if CMOR can have a dependency on a fairly recent version of the pyproj module (>=3.0.0) then writing the crs variables and grid-related attributes in a CF-way should be rather straightforward, see https://pyproj4.github.io/pyproj/dev/build_crs_cf.html

[durack1]

@TomLav thanks for this, really useful. We've been trying to keep dependencies to a minimum, so will ponder this suggestion.

The CMOR software doesn't need to be built from source, you should be able to install it into a conda/mamba environment for your single (or multiple) dataset preparation - https://anaconda.org/conda-forge/cmor