Add support for multiple coord_systems

docs/src/further_topics/netcdf_io.rst

@@ -235,3 +235,70 @@ Worked example:
     >>> my_coord.ignore_axis = True
     >>> print(guess_coord_axis(my_coord))
     None
+
+Multiple Coordinate Systems and Ordered Axes
+--------------------------------------------
+
+Traditionally, the coordinate system for coordinate of a data variable
+is specified by variable a single name in the ``grid_mapping`` attribute,
+e.g.:
+
+::
+
+    grid_mapping = 'latitude_longitude`
+
+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 the concept of **ordered axes** in the definition
+of a coordinate system which allows for more explicit specification of the
+associated coordinate variables.
+
+This is achieved via use of an extended syntax in the ``grid_mapping``
+attribute of a data variable:
+
+::
+
+    <gridMappingVariable>: <coordinatesVariable> [<coordinatesVariable> …] [<gridMappingVariable>: <coordinatesVariable>…]
+
+where each ``gridMappingVariable`` identifies a grid mapping variable
+followed by the list of associated coordinate variables. Note that with
+this syntax it is possible to specify multiple coordinate systems for a
+data variable.
+
+The order of the axes in the extended grid mapping specification is
+significant, but only when used in conjunction with a "well known text"
+(WKT) representation of the coordinate system _REF_ 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.coord.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.ordered_axes` 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. 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.ordered_axes` 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.ordered_axes=True``.

docs/src/whatsnew/latest.rst

@@ -52,6 +52,18 @@ This document explains the changes made to Iris for this release
    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 **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.ordered_coords = True`` will
+    enable extended grid mapping syntax when saving a NetCDF file and also
+    enerated an associated well known text attribute (``crs_wks``; as
+    described in :issue:`3796`).
+    See `CRS Grid Mappings and Projections`_ for more information.
+    (:issue:3388:, :pull:6536:)
+
 
 🐛 Bugs Fixed
 =============
@@ -146,4 +158,5 @@ This document explains the changes made to Iris for this release
     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

lib/iris/cube.py

@@ -1291,6 +1291,9 @@ def __init__(
         ] = []
         self._aux_factories: list[AuxCoordFactory] = []
 
+        # Default ordered_axes property to False; requires explicit opt-in
+        self._ordered_axes = False
+
         # Cell Measures
         self._cell_measures_and_dims: list[tuple[CellMeasure, tuple[int, ...]]] = []
 
@@ -2451,6 +2454,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 ordered_axes(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._ordered_axes
+
+    @ordered_axes.setter
+    def ordered_axes(self, ordered: bool) -> None:
+        """Set to True to enable extended grid mapping syntax."""
+        self._ordered_axes = 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)
 

lib/iris/fileformats/_nc_load_rules/helpers.py

@@ -143,6 +143,46 @@
 CF_GRID_MAPPING_OBLIQUE = "oblique_mercator"
 CF_GRID_MAPPING_ROTATED_MERCATOR = "rotated_mercator"
 
+#
+# Regex for parsing grid_mapping (extended format)
+# Link to online regex101 playground: https://regex101.com/r/NcKzkQ/1
+#
+#   (\w+):                # Matches '<word>:' and stores in CAPTURE GROUP 1
+#   (                     # CAPTURE GROUP 2 for capturing multiple coords
+#       (?:               #  Non-capturing group for composing match
+#           \s+           #   Matches one or more blank characters
+#           (?!\w+:)      #   Negative look-ahead: don't match <word> followed by colon
+#           \w+           #   Matches a <word>
+#       )+                #  Repeats non-capturing group at least once.
+#   )                     # End of CAPTURE GROUP 2
+_GRID_MAPPING_PARSE_EXTENDED = re.compile(
+    r"""
+        (\w+):
+        (
+            (?:
+                \s+
+                (?!\w+:)
+                \w+
+            )+
+        )+
+    """,
+    re.VERBOSE,
+)
+_GRID_MAPPING_PARSE_SIMPLE = re.compile(r"^\w+$")
+_GRID_MAPPING_VALIDATORS = (  # fmt: skip
+    (
+        re.compile(r"\w+: +\w+:"),
+        "`<coord_system>:` identifier followed immediately by another `<coord_system>:` identifier",
+    ),
+    (
+        re.compile(r"\w+: *$"),
+        "`<coord_system>:` is empty - missing coordinate list",
+    ),
+    (
+        re.compile(r"^\w+ +\w+"),
+        "Multiple coordinates found without `<coord_system>:` identifier",
+    ),
+)
 #
 # CF Attribute Names.
 #
@@ -1936,6 +1976,41 @@ def is_grid_mapping(engine, cf_name, grid_mapping):
     return is_valid
 
 
+################################################################################
+def _parse_extended_grid_mapping(grid_mapping: str) -> dict[None | str, str]:
+    """Parse `grid_mapping` attribute and return list of coordinate system variables and associated coords."""
+    # Handles extended grid_mapping too. Possibilities:
+    #  grid_mapping = "crs"  : simple mapping; a single variable name with no coords
+    #  grid_mapping = "crs: lat lon"  : extended mapping; a variable name and list of coords
+    #  grid_mapping = "crs: lat lon other: var1 var2"  : multiple extended mappings
+
+    mappings: dict[None | str, str]
+
+    # try simple mapping first
+    if _GRID_MAPPING_PARSE_SIMPLE.match(grid_mapping):
+        mappings = {None: grid_mapping}  # simple single grid mapping variable
+    else:
+        # Try extended mapping:
+        # 1. Run validators to check for invalid expressions:
+        for v_re, v_msg in _GRID_MAPPING_VALIDATORS:
+            if len(match := v_re.findall(grid_mapping)):
+                msg = f"Invalid syntax in extended grid_mapping: {grid_mapping!r}\n{v_msg} : {match}"
+                raise iris.exceptions.CFParseError(msg)
+
+        # 2. Parse grid_mapping into list of [cs, (coords, ...)]:
+        result = _GRID_MAPPING_PARSE_EXTENDED.findall(grid_mapping)
+        if len(result) == 0:
+            msg = f"Failed to parse grid_mapping: {grid_mapping!r}"
+            raise iris.exceptions.CFParseError(msg)
+
+        # split second match group into list of coordinates:
+        mappings = {}
+        for cs, coords in result:
+            mappings.update({coord: cs for coord in coords.split()})
+
+    return mappings
+
+
 ################################################################################
 def _is_rotated(engine, cf_name, cf_attr_value):
     """Determine whether the CF coordinate variable is rotated."""