Adopt WKT for CRS

[marqh]

Feature Request: Support crs_wkt grid_mapping syntax (CF1.8) in Iris.

As of version 1.8 of the CF Conventions, the use of Well Known Text for Coordinate Reference Systems encodings with standardised variable referencing is supported for CF netCDF files.

The adoption of the ISO/OGC standard for CRS definitions is a major step forward for CF, enabling enhanced interoperability and providing new capabilities.
It would be beneficial for Iris to support this capability effectively to enhance users' ability to effectively provide referencing information for coordinates.

The OGC publishes the joint ISO/OGC standard in the public domain at:
http://docs.opengeospatial.org/is/18-010r7/18-010r7.html

It is worth noting that recent version of Proj support WKT, including parsing of WKT CRS strings.

The use of WKT CRS strings is easily supported, as it just involves the use of the CF attribute name crs_wkt on a variable.
There are minor issues to note on escape characters, particularly the need to escape " characters within the string as \". This is likely easy to handle in software.

The challenge for Iris is the enhancements made to the grid_mapping data variable attribute to enable the references between the data variable, and by implication the associated coordinates.

This enhanced syntax is detailed in the CF Conventions:
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#grid-mappings-and-projections

The grid mapping variables are associated with the data and coordinate variables by the grid_mapping attribute. This attribute is attached to data variables so that variables with different mappings may be present in a single file. The attribute takes a string value with two possible formats. In the first format, it is a single word, which names a grid mapping variable. In the second format, it is a blank-separated list of words ": [ …​] [: …​]" , which identifies one or more grid mapping variables, and with each grid mapping associates one or more coordinatesVariables, i.e. coordinate variables or auxiliary coordinate variables.
Where an extended ": []" entity is defined, then the order of the references within the definition provides an explicit order for these coordinate value variables, which is used if they are to be combined into individual coordinate tuples.
This order is only significant if crs_wkt is also specified within the referenced grid mapping variable. Explicit 'axis order' is important when the grid_mapping_variable contains an attribute crs_wkt as it is mandated by the OGC CRS-WKT standard that coordinate tuples with correct axis order are provided as part of the reference to a Coordinate Reference System.

Iris does not recognise this syntax and does not interpret any coordinate reference system referencing where it is used.

The enhanced syntax has particular value to users using Iris to create data sets in CF-netCDF for sharing with wider communities.
The adoption of the CF syntax for these relations and their interpretation in the Iris data model brings valuable opportunities for improved referencing by coordinates and coordinate reference system definitions for data specification.

Please may we develop this ticket into a full enhancement proposal and aim to plan work to deliver this capability in an Iris release?

Tasks

Beta Give feedback

1. Save crs_wkt on netCDF save #4719
   Dragon Sub-Task 🦎 +1
   
2. Support loading/saving netCDF with multiple coordinate systems #3388
   Dragon Sub-Task 🦎 Feature: NetCDF + CF-conventions Release: Minor Type: Enhancement +4
3. Make modifying Cube CoordSystems easier #6030
   Dragon Sub-Task 🦎 +1

[marqh]

I have prepared a sample data file which illustrates the use of the grid_mapping attribute and the crs_wkt attribute.

This is a CDL file, to aid readability. The nD variables are left unspecified, saving space.

ncgen -k 4 -o CO2_GEO.nc CO2_GEO.cdl

will create a binary netCDF file with these unspecified variables encoded as missing data.

The CDL specification for the example file is:

netcdf CO2_GEO {
dimensions:
	longitude = 360 ;
	latitude = 181 ;
	levelist = 60 ;
	time = 3 ;
variables:
	float longitude(longitude) ;
		longitude:units = "degrees_east" ;
		longitude:standard_name = "longitude" ;
	float latitude(latitude) ;
		latitude:units = "degrees_north" ;
		latitude:standard_name = "latitude" ;
	int altitude(levelist, latitude, longitude) ;
                altitude:units = "metres" ;
                altitude:standard_name = "altitude" ;
        int levelist(levelist) ;
		levelist:long_name = "model_level_number" ;
	int time(time) ;
		time:units = "hours since 2000-01-01 00:00" ;
		time:standard_name = "time" ;
	short co2(time, levelist, latitude, longitude) ;
		co2:scale_factor = 0.000981685145029486 ;
		co2:add_offset = 403.192219379918 ;
		co2:_FillValue = -32767s ;
		co2:missing_value = -32767s ;
		co2:units = "kg kg**-1" ;
		co2:long_name = "Carbon Dioxide" ;
		co2:standard_name = "mass_fraction_of_carbon_dioxide_in_air" ;
                co2:grid_mapping = "WGS84_2D: latitude longitude EGM2008: altitude Gregorian: time" ;
	short lnsp(time, levelist, latitude, longitude) ;
		lnsp:add_offset = 11.2087164280841 ;
		lnsp:_FillValue = -32767s ;
		lnsp:missing_value = -32767s ;
		lnsp:long_name = "Logarithm of surface pressure" ;
                lnsp:grid_mapping = "WGS84_2D: latitude longitude EGM2008: altitude Gregorian: time" ;

        int WGS84_2D ;
                WGS84_2D: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,2],    AXIS[\"latitude\",north,ORDER[1]],    AXIS[\"longitude\",east,ORDER[2]],    ANGLEUNIT[\"degree\",0.01745329252],  ID[\"EPSG\",4326]]" ;
              WGS84_2D:grid_mapping_name = "latitude_longitude";
              WGS84_2D:longitude_of_prime_meridian = 0.0 ;
              WGS84_2D:semi_major_axis = 6378137.0 ;
              WGS84_2D:inverse_flattening = 298.257223563 ;

        int EGM2008;
                EGM2008:crs_wkt = "VERTCRS[\"EGM2008 height\",  VDATUM[\"EGM2008 geoid\"],  CS[vertical,1],    AXIS[\"gravity-related height (H)\",up],    LENGTHUNIT[\"metre\",1.0],  ID[\"EPSG\",3855]]" ;

        int Gregorian;
                Gregorian:crs_wkt = "TIMECRS[\"Calendar hours from 2000-01-01\",  TDATUM[\"01 January 2000\",TIMEORIGIN[2000-01-01T00:00Z]],  CS[TemporalCount,1],AXIS[\"Time\",future,TIMEUNIT[\"hour\"]]]" ;

// global attributes:
		:Conventions = "CF-1.8" ;

data:

 longitude = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
    18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 
    36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 
    54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 
    72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 
    90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 
    106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 
    120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 
    134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 
    148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 
    162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 
    176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 
    190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 
    204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 
    218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 
    232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 
    246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 
    260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 
    274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 
    288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 
    302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 
    316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 
    330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 
    344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 
    358, 359 ;

 latitude = 90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 
    74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 
    56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 
    38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 
    20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 
    -1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, -13, -14, -15, -16, 
    -17, -18, -19, -20, -21, -22, -23, -24, -25, -26, -27, -28, -29, -30, 
    -31, -32, -33, -34, -35, -36, -37, -38, -39, -40, -41, -42, -43, -44, 
    -45, -46, -47, -48, -49, -50, -51, -52, -53, -54, -55, -56, -57, -58, 
    -59, -60, -61, -62, -63, -64, -65, -66, -67, -68, -69, -70, -71, -72, 
    -73, -74, -75, -76, -77, -78, -79, -80, -81, -82, -83, -84, -85, -86, 
    -87, -88, -89, -90 ;

 levelist = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
    19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 
    37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 
    55, 56, 57, 58, 59, 60 ;

 time = 287088, 287100, 287112 ;
}

[pp-mo]

I agree that users will probably want this soon, and this is a really useful start.

However the new, multiple form of the "grid_mapping" attributes is I think better handled as a separate issue.
We do already have a ticket for that #3388

Let's by all means cross-couple them, but I think the effects relating to "crs_wkt" are only a minor aspect in that context.
Practically, I think it will be clearest to address #3388 first, and add crs_wkt support afterwards.

[marqh]

Hi @pp-mo
thank you for the input and cross referencing, i had not realised the issue overlap, which is now clear

I agree with your assessment:

Let's by all means cross-couple them, but I think the effects relating to "crs_wkt" are only a minor aspect in that context.
Practically, I think it will be clearest to address #3388 first, and add crs_wkt support afterwards.

The crs_wkt support is essentially just a string attribute. It will be fairly easy to address this as a simple addition, and can be done after #3388 quite sensibly. It is likely to be low impact, I suspect.

The capability that is missing in Iris and limiting here is the grid_mapping handling for extended cases.

Bearing this ticket in mind whilst addressing the enhanced grid_mapping syntax is useful, as the order of specified coordinates within an enhanced grid_mapping attribute is of key importance in CF1.8, but not clear in CF1.7.

So a solution that preserves the order of terms, and handles the references such that this order can be used is the important factor here.

a:grid_mapping = "WGS84_2D: latitude longitude != a:grid_mapping = "WGS84_2D: longitude latitude
(the latter is wrong, so i crossed it out)

So implementing #3388 in an unordered fashion could meet the scope of #3388 as currently stated, but would be a real limiting factor for addressing this capability request.

I hope this is useful confirmation and clarification of your perspective.

[github-actions]

In order to maintain a backlog of relevant issues, we automatically label them as stale after 500 days of inactivity.

If this issue is still important to you, then please comment on this issue and the stale label will be removed.

Otherwise this issue will be automatically closed in 28 days time.

[larsbarring]

I think that this would still be a useful extension of Iris, despite that the more general question dealt with in #3388 was closed.

[pp-mo]

I think that this would still be a useful extension of Iris, despite that the more general question dealt with in #3388 was closed.

@SciTools/peloton the real win would be #4643
but that is still rather a long way off.
Unless you have a really good driver case for this??

[trexfeathers]

We have [mistakenly] conducted a detailed conversation on the pull request, not the issue - see #4719.

It has basically been decided that this overlaps enough with #3388 that #3388 will also be required before this work is considered completed. See the tasklist at the top of this issue for more detail.