Merge pull request #6537 from SciTools/FEATURE_wkt

Merge `FEATURE_wkt` into `main`
Provides support for:
  - Generation of `crs_wkt` (well known text) attribute on netCDF output
  - Reading and writing of multiple coordinate systems via extended grid mapping support.

Author: ukmo-ccbunney

docs/src/further_topics/netcdf_io.rst

@@ -235,3 +235,160 @@ Worked example:
     >>> my_coord.ignore_axis = True
     >>> print(guess_coord_axis(my_coord))
     None
+
+Multiple Coordinate Systems and Ordered Axes
+--------------------------------------------
+
+In a CF compliant NetCDF file, the coordinate variables associated with a
+data variable can specify a specific *coordinate system* that defines how
+the coordinate values relate to physical locations on the globe. For example,
+a coordinate might have values with units of metres that should be referenced
+against a *Transverse Mercator* projection with a specific origin. This
+information is not stored on the coordinate itself, but in a separate
+*grid mapping* variable. Furthermore, the grid mapping for a set of
+coordinates is associated with the data variable (not the coordinates
+variables) via the ``grid_mapping`` attribute.
+
+For example, a temperature variable defined on a *rotated pole* grid might
+look like this in a NetCDF file (extract of relevant variables):
+
+.. code-block:: text
+
+  float T(rlat,rlon) ;
+    T:long_name = "temperature" ;
+    T:units = "K" ;
+    T:grid_mapping = "rotated_pole" ;
+
+  char rotated_pole ;
+    rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
+    rotated_pole:grid_north_pole_latitude = 32.5 ;
+    rotated_pole:grid_north_pole_longitude = 170. ;
+
+  float rlon(rlon) ;
+    rlon:long_name = "longitude in rotated pole grid" ;
+    rlon:units = "degrees" ;
+    rlon:standard_name = "grid_longitude";
+
+  float rlat(rlat) ;
+    rlat:long_name = "latitude in rotated pole grid" ;
+    rlat:units = "degrees" ;
+    rlat:standard_name = "grid_latitude";
+
+
+Note how the ``rotated pole`` grid mapping (coordinate system) is referenced
+from the data variable ``T:grid_mapping = "rotated_pole"`` and is implicitly
+associated with the dimension coordinate variables ``rlat`` and ``rlon``.
+
+
+Since version `1.8 of the CF Conventions
+<https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#grid-mappings-and-projections>`_
+, there has been support for a more explicit version of the ``grid_mapping``
+attribute. This allows for **multiple coordinate systems** to be defined for
+a data variable and individual coordinates to be explicitly associated with
+a coordinate system. This is achieved by use of an **extended syntax** in the
+``grid_mapping`` variable of a data variable:
+
+
+.. code-block:: text
+
+  <grid_mapping_var>: <coord_var> [<coord_var>] [<grid_mapping_var>: <coord_var> ...]
+
+where each ``grid_mapping_var`` identifies a grid mapping variable followed by
+the list of associated coordinate variables (``coord_var``). Note that with
+this syntax it is possible to specify multiple coordinate systems for a
+data variable.
+
+For example, consider the following *air pressure* variable that is
+defined on an *OSGB Transverse Mercator grid*:
+
+.. code-block:: text
+
+    float pres(y, x) ;
+        pres:standard_name = "air_pressure" ;
+        pres:units = "Pa" ;
+        pres:coordinates = "lat lon" ;
+        pres:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+
+    double x(x) ;
+        x:standard_name = "projection_x_coordinate" ;
+        x:units = "m" ;
+
+    double y(y) ;
+        y:standard_name = "projection_y_coordinate" ;
+        y:units = "m" ;
+
+    double lat(y, x) ;
+        lat:standard_name = "latitude" ;
+        lat:units = "degrees_north" ;
+
+    double lon(y, x) ;
+        lon:standard_name = "longitude" ;
+        lon:units = "degrees_east" ;
+
+    int crsOSGB ;
+        crsOSGB:grid_mapping_name = "transverse_mercator" ;
+        crsOSGB:semi_major_axis = 6377563.396 ;
+        crsOSGB:inverse_flattening = 299.3249646 ;
+        <snip>
+
+    int crsWGS84 ;
+        crsWGS84:grid_mapping_name = "latitude_longitude" ;
+        crsWGS84:longitude_of_prime_meridian = 0. ;
+        <snip>
+
+
+The dimension coordinates ``x`` and ``y`` are explicitly defined on
+an a *transverse mercator* grid via the ``crsOSGB`` variable.
+
+However, with the extended grid syntax, it is also possible to define
+a second coordinate system on a standard **latitude_longitude** grid
+and associate it with the auxiliary ``lat`` and ``lon`` coordinates:
+
+::
+
+    pres:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
+
+
+Note, the *order* of the axes in the extended grid mapping specification is
+significant, but only when used in conjunction with a
+`CRS Well Known Text (WKT)`_ representation of the coordinate system where it
+should be consistent with the ``AXES ORDER`` specified in the ``crs_wkt``
+attribute.
+
+
+Effect on loading
+^^^^^^^^^^^^^^^^^
+
+When Iris loads a NetCDF file that uses the extended grid mapping syntax
+it will generate an :class:`iris.coord_systems.CoordSystem` for each
+coordinate system listed and attempt to attach it to the associated
+:class:`iris.coords.Coord` instances on the cube. Currently, Iris considers
+the ``crs_wkt`` supplementary and builds coordinate systems exclusively
+from the ``grid_mapping`` attribute.
+
+The :attr:`iris.cube.Cube.extended_grid_mapping` property will be set to
+``True`` for cubes loaded from NetCDF data variables utilising the extended
+``grid_mapping`` syntax.
+
+Effect on saving
+^^^^^^^^^^^^^^^^
+
+To maintain existing behaviour, saving an :class:`iris.cube.Cube` to
+a netCDF file will default to the "simple" grid mapping syntax, unless
+the cube was loaded from a file using the extended grid mapping syntax.
+If the cube contains multiple coordinate systems, only the coordinate
+system of the dimension coordinate(s) will be specified.
+
+To enable saving of multiple coordinate systems with ordered axes,
+set the :attr:`iris.cube.Cube.extended_grid_mapping` to ``True``.
+This will generate a ``grid_mapping`` attribute using the extended syntax
+to specify all coordinate systems on the cube. The axes ordering of the
+associated coordinate variables will be consistent with that of the
+generated ``crs_wkt`` attribute.
+
+Note, the ``crs_wkt`` attribute will only be generated when the
+extended grid mapping is also written, i.e. when
+``Cube.extended_grid_mapping=True``.
+
+
+.. _CRS Well Known Text (WKT): https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#use-of-the-crs-well-known-text-format

docs/src/whatsnew/latest.rst

@@ -40,6 +40,31 @@ This document explains the changes made to Iris for this release
    now preserved when it was not previously. See also: :ref:`load-problems`.
    (:pull:`6465`, :pull:`6529`)
 
+#. `@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`)
+
+#. `@ukmo-ccbunney`_ and `@trexfeathers`_ added support for
+   **multiple coordinate systems** and **ordered coordinates** when loading
+   and saving NetCDF files.
+   This allows for coordinates to be explicitly associated with a coordinate
+   system via an extended syntax in the ``grid_mapping`` attribute of a NetCDF
+   data variable. This extended syntax also supports specification of multiple
+   coordinate systems per data variable. Setting the property
+   ``cube.extended_grid_mapping = True`` will enable extended grid mapping
+   syntax when saving a NetCDF file and also generate an associated **well known
+   text** attribute (``crs_wkt``; as described in :issue:`3796`).
+   See `CRS Grid Mappings and Projections`_ for more information.
+   (:issue:`3388`:, :pull:`6536`:)
+
 #. `@ESadek-MO`_ made MeshCoords immutable. :class:`iris.MeshCoord`s are now updated automatically when
    changing the attached mesh. All changes to the :class:`iris.MeshCoord` should instead be done to
    the relevant :class:`iris.Coord` located on the attached :class:`iris.MeshXY`. This change also affects 
@@ -166,4 +191,7 @@ 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
+.. _CRS Grid Mappings and Projections: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#grid-mappings-and-projections
+.. _Ncdata: https://github.com/pp-mo/ncdata
 .. _Trusted Publishing: https://docs.pypi.org/trusted-publishers/

lib/iris/cube.py

@@ -2461,6 +2461,32 @@ def coord_system(
 
         return result
 
+    def coord_systems(self) -> list[iris.coord_systems.CoordSystem]:
+        """Return a list of all coordinate systems used in cube coordinates."""
+        # Gather list of our unique CoordSystems on cube:
+        coord_systems = ClassDict(iris.coord_systems.CoordSystem)
+        for coord in self.coords():
+            if coord.coord_system:
+                coord_systems.add(coord.coord_system, replace=True)
+
+        return list(coord_systems.values())
+
+    @property
+    def extended_grid_mapping(self) -> bool:
+        """Return True if a cube will use extended grid mapping syntax to write axes order in grid_mapping.
+
+        Only relevant when saving a cube to NetCDF file format.
+
+        For more details see "Grid Mappings and Projections" in the CF Conventions document:
+        https://cfconventions.org/cf-conventions/conformance.html
+        """
+        return self.attributes.get("iris_extended_grid_mapping", False)
+
+    @extended_grid_mapping.setter
+    def extended_grid_mapping(self, ordered: bool) -> None:
+        """Set to True to enable extended grid mapping syntax."""
+        self.attributes["iris_extended_grid_mapping"] = ordered
+
     def _any_meshcoord(self) -> MeshCoord | None:
         """Return a MeshCoord if there are any, else None."""
         mesh_coords = self.coords(mesh_coords=True)

lib/iris/exceptions.py

@@ -172,3 +172,9 @@ def __str__(self):
             "operations are not currently supported."
         )
         return msg.format(super().__str__())
+
+
+class CFParseError(IrisError):
+    """Raised when a string associated with a CF defined syntax could not be parsed."""
+
+    pass

lib/iris/fileformats/_nc_load_rules/actions.py

@@ -203,7 +203,10 @@ def action_provides_grid_mapping(engine, gridmapping_fact):
 
         def build_outer(engine_, cf_var_):
             coordinate_system = builder(engine_, cf_var_)
-            engine_.cube_parts["coordinate_system"] = coordinate_system
+            # We can now handle more than one coordinate_system, so store as dictionary:
+            engine_.cube_parts["coordinate_systems"][cf_var_.cf_name] = (
+                coordinate_system
+            )
 
         # Part 1 - only building - adding takes place downstream in
         #  helpers.build_and_add_dimension/auxiliary_coordinate().
@@ -218,11 +221,8 @@ def build_outer(engine_, cf_var_):
             ),
         )
 
-        # Check there is not an existing one.
-        # ATM this is guaranteed by the caller, "run_actions".
-        assert engine.fact_list("grid-type") == []
-
-        engine.add_fact("grid-type", (grid_mapping_type,))
+        # Store grid-mapping name along with grid-type to match them later on
+        engine.add_fact("grid-type", (var_name, grid_mapping_type))
 
     else:
         message = "Coordinate system not created. Debug info:\n"
@@ -343,7 +343,35 @@ def action_build_dimension_coordinate(engine, providescoord_fact):
         # Non-conforming lon/lat/projection coords will be classed as
         # dim-coords by cf.py, but 'action_provides_coordinate' will give them
         # a coord-type of 'miscellaneous' : hence, they have no coord-system.
-        coord_system = engine.cube_parts.get("coordinate_system")
+        #
+        # At this point, we need to match any "coordinate_system" entries in
+        # the engine to the coord we are building. There are a couple of cases here:
+        #  1. Simple `grid_mapping = crs` is used, in which case
+        #     we should just apply that mapping to all dim coords.
+        #  2. Extended `grid_mapping = crs: coord1 coord2 crs: coord3 coord4`
+        #     is used in which case we need to match the crs to the coord here.
+
+        # We can have multiple coordinate_system, so now stored as a list (note plural key)
+        coord_systems = engine.cube_parts.get("coordinate_systems")
+        coord_system = None
+
+        if len(coord_systems):
+            # Find which coord system applies to this coordinate.
+            cs_mappings = engine.cube_parts.get("coordinate_system_mappings")
+            if cs_mappings and coord_systems:
+                if len(coord_systems) == 1 and None in cs_mappings:
+                    # Simple grid mapping (a single coord_system with no explicit coords)
+                    # Applies to spatial DimCoord(s) only. In this case only one
+                    # coordinate_system will have been built, so just use it.
+                    (coord_system,) = coord_systems.values()
+                    (cs_name,) = cs_mappings.values()
+                else:
+                    # Extended grid mapping, e.g.
+                    #  `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
+                    # We need to search for coord system that references our coordinate.
+                    if cs_name := cs_mappings.get(cf_var.cf_name):
+                        coord_system = coord_systems.get(cs_name, None)
+
         # Translate the specific grid-mapping type to a grid-class
         if coord_system is None:
             succeed = True
@@ -352,8 +380,13 @@ def action_build_dimension_coordinate(engine, providescoord_fact):
             # Get a grid-class from the grid-type
             # i.e. one of latlon/rotated/projected, as for coord_grid_class.
             gridtypes_factlist = engine.fact_list("grid-type")
-            (gridtypes_fact,) = gridtypes_factlist  # only 1 fact
-            (cs_gridtype,) = gridtypes_fact  # fact contains 1 term
+
+            # potentially multiple grid-type facts; find one for CRS varname
+            cs_gridtype = None
+            for fact_cs_name, fact_cs_type in gridtypes_factlist:
+                if fact_cs_name == cs_name:
+                    cs_gridtype = fact_cs_type
+
             if cs_gridtype == "latitude_longitude":
                 cs_gridclass = "latlon"
             elif cs_gridtype == "rotated_latitude_longitude":
@@ -446,6 +479,7 @@ def action_build_auxiliary_coordinate(engine, auxcoord_fact):
     """Convert a CFAuxiliaryCoordinateVariable into a cube aux-coord."""
     (var_name,) = auxcoord_fact
     rule_name = "fc_build_auxiliary_coordinate"
+    cf_var = engine.cf_var.cf_group[var_name]
 
     # Identify any known coord "type" : latitude/longitude/time/time_period
     # If latitude/longitude, this sets the standard_name of the built AuxCoord
@@ -473,8 +507,27 @@ def action_build_auxiliary_coordinate(engine, auxcoord_fact):
     if coord_type:
         rule_name += f"_{coord_type}"
 
+    # Check if we have a coord_system specified for this coordinate.
+    # (Only possible via extended grid_mapping attribute)
+    coord_systems = engine.cube_parts.get("coordinate_systems")
+    coord_system = None
+
+    cs_mappings = engine.cube_parts.get("coordinate_system_mappings", None)
+    if cs_mappings and coord_systems:
+        if len(coord_systems) == 1 and None in cs_mappings:
+            # Simple grid_mapping - doesn't apply to AuxCoords (we need an explicit mapping)
+            pass
+        else:
+            # Extended grid mapping, e.g.
+            #  `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
+            # We need to search for coord system that references our coordinate.
+            if cs_name := cs_mappings.get(cf_var.cf_name):
+                coord_system = coord_systems.get(cs_name, None)
+
     cf_var = engine.cf_var.cf_group.auxiliary_coordinates[var_name]
-    hh.build_and_add_auxiliary_coordinate(engine, cf_var, coord_name=coord_name)
+    hh.build_and_add_auxiliary_coordinate(
+        engine, cf_var, coord_name=coord_name, coord_system=coord_system
+    )
 
     return rule_name
 
@@ -615,10 +668,9 @@ def run_actions(engine):
     # default (all cubes) action, always runs
     action_default(engine)  # This should run the default rules.
 
-    # deal with grid-mappings
+    # deal with grid-mappings; potentially multiple mappings if extended grid_mapping used.
     grid_mapping_facts = engine.fact_list("grid_mapping")
-    # For now, there should be at most *one* of these.
-    assert len(grid_mapping_facts) in (0, 1)
+
     for grid_mapping_fact in grid_mapping_facts:
         action_provides_grid_mapping(engine, grid_mapping_fact)