Save crs_wkt on netCDF save, take 2 (#6519)

* Save crs_wkt, and update test outputs

* Extra test fixes

* What's new

* Update tests.

* What's New updates.

* Correct use of mocking.

---------

Co-authored-by: Will Benfold <[email protected]>

Author: trexfeathers

docs/src/whatsnew/latest.rst

@@ -39,6 +39,18 @@ This document explains the changes made to Iris for this release
    converted to :class:`~iris.coords.AuxCoord` if it is masked - the mask is
    now preserved when it was not previously.
 
+#. `@wjbenfold`_ and `@trexfeathers`_ added ``crs_wkt`` to the attributes when
+   saving a :class:`~iris.coord_systems.CoordSystem` to a NetCDF file. Note that
+   ``crs_wkt`` is considered *supplementary* by the CF conventions, with
+   ``grid_mapping`` being the primary source of information, and ``crs_wkt`` not
+   expected to contain conflicting information. Because of this, Iris generates
+   :class:`~iris.coord_systems.CoordSystem` exclusively from ``grid_mapping``
+   when loading, and writes a fresh ``crs_wkt`` whenever a
+   :class:`~iris.coord_systems.CoordSystem` is saved. If your use case goes
+   beyond the CF conventions, you can modify the save and load process for your
+   needs by using the `Ncdata`_ package.
+   See `CRS WKT in the CF Conventions`_ for more. (:issue:`3796`, :pull:`6519`)
+
 
 🐛 Bugs Fixed
 =============
@@ -105,3 +117,6 @@ This document explains the changes made to Iris for this release
 
 .. comment
     Whatsnew resources in alphabetical order:
+
+.. _CRS WKT in the CF Conventions: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#use-of-the-crs-well-known-text-format
+.. _Ncdata: https://github.com/pp-mo/ncdata

lib/iris/fileformats/netcdf/saver.py

@@ -2081,6 +2081,8 @@ def add_ellipsoid(ellipsoid):
                         category=iris.warnings.IrisSaveWarning,
                     )
 
+                cf_var_grid.crs_wkt = cs.as_cartopy_crs().to_wkt()
+
                 self._coord_systems.append(cs)
 
             # Refer to grid var

lib/iris/tests/results/integration/netcdf/aux_factories/TestAtmosphereSigma/save.cdl

@@ -17,6 +17,7 @@ variables:
 		rotated_latitude_longitude:grid_north_pole_latitude = 37.5 ;
 		rotated_latitude_longitude:grid_north_pole_longitude = 177.5 ;
 		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;
+		rotated_latitude_longitude:crs_wkt = "GEOGCRS[\"unnamed\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],DERIVINGCONVERSION[\"unknown\",METHOD[\"PROJ ob_tran o_proj=latlon\"],PARAMETER[\"o_lon_p\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],PARAMETER[\"o_lat_p\",37.5,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],PARAMETER[\"lon_0\",357.5,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	double time(time) ;
 		time:axis = "T" ;
 		time:units = "hours since 1970-01-01 00:00:00" ;

lib/iris/tests/results/integration/netcdf/aux_factories/TestHybridPressure/save.cdl

@@ -17,6 +17,7 @@ variables:
 		rotated_latitude_longitude:grid_north_pole_latitude = 37.5 ;
 		rotated_latitude_longitude:grid_north_pole_longitude = 177.5 ;
 		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;
+		rotated_latitude_longitude:crs_wkt = "GEOGCRS[\"unnamed\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],DERIVINGCONVERSION[\"unknown\",METHOD[\"PROJ ob_tran o_proj=latlon\"],PARAMETER[\"o_lon_p\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],PARAMETER[\"o_lat_p\",37.5,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],PARAMETER[\"lon_0\",357.5,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	double time(time) ;
 		time:axis = "T" ;
 		time:units = "hours since 1970-01-01 00:00:00" ;

lib/iris/tests/results/integration/netcdf/aux_factories/TestSaveMultipleAuxFactories/hybrid_height_and_pressure.cdl

@@ -17,6 +17,7 @@ variables:
 		rotated_latitude_longitude:grid_north_pole_latitude = 37.5 ;
 		rotated_latitude_longitude:grid_north_pole_longitude = 177.5 ;
 		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;
+		rotated_latitude_longitude:crs_wkt = "GEOGCRS[\"unnamed\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],DERIVINGCONVERSION[\"unknown\",METHOD[\"PROJ ob_tran o_proj=latlon\"],PARAMETER[\"o_lon_p\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],PARAMETER[\"o_lat_p\",37.5,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],PARAMETER[\"lon_0\",357.5,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	double time(time) ;
 		time:axis = "T" ;
 		time:units = "hours since 1970-01-01 00:00:00" ;

lib/iris/tests/results/integration/netcdf/general/TestPackedData/multi_packed_multi_dtype.cdl

@@ -17,6 +17,7 @@ variables:
 		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
 		latitude_longitude:longitude_of_prime_meridian = 0. ;
 		latitude_longitude:earth_radius = 6371229. ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	double time(time) ;
 		time:axis = "T" ;
 		time:bounds = "time_bnds" ;

lib/iris/tests/results/integration/netcdf/general/TestPackedData/multi_packed_single_dtype.cdl

@@ -17,6 +17,7 @@ variables:
 		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
 		latitude_longitude:longitude_of_prime_meridian = 0. ;
 		latitude_longitude:earth_radius = 6371229. ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	double time(time) ;
 		time:axis = "T" ;
 		time:bounds = "time_bnds" ;

lib/iris/tests/results/integration/netcdf/general/TestPackedData/single_packed_manual.cdl

@@ -16,6 +16,7 @@ variables:
 		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
 		latitude_longitude:longitude_of_prime_meridian = 0. ;
 		latitude_longitude:earth_radius = 6371229. ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	float latitude(latitude) ;
 		latitude:axis = "Y" ;
 		latitude:units = "degrees_north" ;

lib/iris/tests/results/integration/netcdf/general/TestPackedData/single_packed_signed.cdl

@@ -16,6 +16,7 @@ variables:
 		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
 		latitude_longitude:longitude_of_prime_meridian = 0. ;
 		latitude_longitude:earth_radius = 6371229. ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	float latitude(latitude) ;
 		latitude:axis = "Y" ;
 		latitude:units = "degrees_north" ;

lib/iris/tests/results/integration/netcdf/general/TestPackedData/single_packed_unsigned.cdl

@@ -16,6 +16,7 @@ variables:
 		latitude_longitude:grid_mapping_name = "latitude_longitude" ;
 		latitude_longitude:longitude_of_prime_meridian = 0. ;
 		latitude_longitude:earth_radius = 6371229. ;
+		latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6371229,0,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]],CS[ellipsoidal,2],AXIS[\"longitude\",east,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"latitude\",north,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]]" ;
 	float latitude(latitude) ;
 		latitude:axis = "Y" ;
 		latitude:units = "degrees_north" ;