Add support for multiple coord_systems

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,42 @@ 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.
+            # Parse the grid_mapping attribute to get coord_system -> coordinate mappings
+            attr_grid_mapping = getattr(engine.cf_var, "grid_mapping", None)
+            if attr_grid_mapping:
+                cs_mappings = hh._parse_extened_grid_mapping(attr_grid_mapping)
+
+                # Simple `grid_mapping = "crs"`
+                # Only one coord_system will be present and cs_grid_mapping will
+                # contain no coordinate references (set to None).
+                if len(coord_systems) == 1 and cs_mappings[0][1] is None:
+                    # Only one grid mapping - apply it.
+                    coord_system = list(coord_systems.values())[0]
+                    cs_name = cs_mappings[0][0]
+
+                # Extended `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
+                # We need to search for coord system that references our coordinate.
+                else:
+                    for name, ref_coords in cs_mappings:
+                        if cf_var.cf_name in ref_coords:
+                            cs_name = name
+                            coord_system = coord_systems[cs_name]
+                            break
+
         # Translate the specific grid-mapping type to a grid-class
         if coord_system is None:
             succeed = True
@@ -352,8 +387,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 +486,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 +514,30 @@ 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
+
+    # get grid_mapping from data variable attribute and parse it
+    grid_mapping_attr = getattr(engine.cf_var, "grid_mapping", None)
+    if grid_mapping_attr:
+        cs_mappings = hh._parse_extened_grid_mapping(grid_mapping_attr)
+
+        if len(coord_systems) == 1 and cs_mappings[0][1] is None:
+            # Simple grid_mapping - doesn't apply to AuxCoords (we need an explicit mapping)
+            pass
+        else:
+            # Extended grid_mapping
+            for crs_name, coords in cs_mappings:
+                if cf_var.cf_name in coords:
+                    coord_system = coord_systems[crs_name]
+                    break
+
     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 +678,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,43 @@
 CF_GRID_MAPPING_OBLIQUE = "oblique_mercator"
 CF_GRID_MAPPING_ROTATED_MERCATOR = "rotated_mercator"
 
+#
+# Regex for parsing grid_mapping (extended format)
+#
+#   (\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 = re.compile(r"(\w+):((?: +(?!\w+:)(\w+))+)+")
+_GRID_MAPPING_PARSE_EXTENDED = re.compile(
+    r"""
+        (\w+):
+        (
+            (?:
+                \s+
+                (?!\w+:)
+                (\w+)
+            )+
+        )+
+    """,
+    re.VERBOSE,
+)
+_GRID_MAPPING_PARSE_SIMPLE = re.compile(r"^\w+$")
+_GRID_MAPPING_VALIDATORS = (
+    (
+        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 +1973,38 @@ def is_grid_mapping(engine, cf_name, grid_mapping):
     return is_valid
 
 
+################################################################################
+def _parse_extened_grid_mapping(grid_mapping):
+    """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
+
+    # TODO(ChrisB): TESTS!!
+
+    # try simple mapping first
+    if _GRID_MAPPING_PARSE_SIMPLE.match(grid_mapping):
+        return [(grid_mapping, None)]  # simple single grid mapping variable
+
+    # 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.IrisError(msg)  # TODO: Better Exception type
+
+    # 2. Parse grid_mapping into list of [cs, (coords, ...)]:
+    mappings = _GRID_MAPPING_PARSE_EXTENDED.findall(grid_mapping)
+    if len(mappings) == 0:
+        msg = f"Failed to parse grid_mapping: {grid_mapping!r}"
+        raise iris.exceptions.IrisError(msg)  # TODO: Better exception type
+
+    # split second match group into list of coordinates:
+    mappings = [(m[0], m[1].split()) for m in mappings]
+    return mappings
+
+
 ################################################################################
 def _is_rotated(engine, cf_name, cf_attr_value):
     """Determine whether the CF coordinate variable is rotated."""

lib/iris/fileformats/cf.py

@@ -24,6 +24,7 @@
 import numpy as np
 import numpy.ma as ma
 
+import iris.fileformats._nc_load_rules.helpers as hh
 from iris.fileformats.netcdf import _thread_safe_nc
 from iris.mesh.components import Connectivity
 import iris.util
@@ -666,17 +667,31 @@ def identify(cls, variables, ignore=None, target=None, warn=True):
             if nc_var_att is not None:
                 name = nc_var_att.strip()
 
-                if name not in ignore:
-                    if name not in variables:
-                        if warn:
-                            message = "Missing CF-netCDF grid mapping variable %r, referenced by netCDF variable %r"
-                            warnings.warn(
-                                message % (name, nc_var_name),
-                                category=iris.warnings.IrisCfMissingVarWarning,
-                            )
-                    else:
-                        result[name] = CFGridMappingVariable(name, variables[name])
+                # parse the grid_mappings
+                mappings = hh._parse_extened_grid_mapping(name)
 
+                for name, coords in mappings:
+                    if name not in ignore:
+                        if name not in variables:
+                            if warn:
+                                message = "Missing CF-netCDF grid mapping variable %r, referenced by netCDF variable %r"
+                                warnings.warn(
+                                    message % (name, nc_var_name),
+                                    category=iris.warnings.IrisCfMissingVarWarning,
+                                )
+                        else:
+                            # For extended grid_mapping, also check coord references exist:
+                            if coords:
+                                for coord_name in coords:
+                                    if coord_name not in variables:
+                                        message = "Missing CF-netCDF coordinate variable %r (associated with grid mapping variable %r), referenced by netCDF variable %r"
+                                        warnings.warn(
+                                            message % (coord_name, name, nc_var_name),
+                                            category=iris.warnings.IrisCfMissingVarWarning,
+                                        )
+                            # TODO: Question: A missing coord reference will not stop the coord_system from
+                            #  being added as a CFGridMappingVariable. Is this ok?
+                            result[name] = CFGridMappingVariable(name, variables[name])
         return result
 
 

lib/iris/fileformats/netcdf/loader.py

@@ -79,6 +79,7 @@ def _assert_case_specific_facts(engine, cf, cf_group):
     engine.cube_parts["coordinates"] = []
     engine.cube_parts["cell_measures"] = []
     engine.cube_parts["ancillary_variables"] = []
+    engine.cube_parts["coordinate_systems"] = {}
 
     # Assert facts for CF coordinates.
     for cf_name in cf_group.coordinates.keys():

lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py

@@ -27,6 +27,7 @@ def _make_testcase_cdl(
         mapping_missingradius=False,
         mapping_type_name=None,
         mapping_scalefactor=None,
+        phenom_grid_mapping="grid",
         yco_values=None,
         xco_name=None,
         yco_name=None,
@@ -226,7 +227,7 @@ def _make_testcase_cdl(
                 double phenom({ydim_name}, {xdim_name}) ;
                     phenom:standard_name = "air_temperature" ;
                     phenom:units = "K" ;
-                    phenom:grid_mapping = "grid" ;
+                    phenom:grid_mapping = "{phenom_grid_mapping}" ;
 {phenom_coords_string}
                 double yco({ydim_name}) ;
                     yco:axis = "Y" ;
@@ -755,6 +756,45 @@ def test_mapping__mismatch__nonll_coords_missing_system(self):
         )
         self.check_result(result, cube_no_cs=True, xco_stdname=False, yco_stdname=False)
 
+    def test_extended_mapping_basic_latlon(self):
+        # A basic reference example with a lat-long grid, but using extended
+        # grid mapping syntax on phenomenon.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(latitude)
+        #     004 : fc_provides_coordinate_(longitude)
+        #     005 : fc_build_coordinate_(latitude)
+        #     006 : fc_build_coordinate_(longitude)
+        # Notes:
+        #     * grid-mapping identified : regular latlon
+        #     * dim-coords identified : lat+lon
+        #     * coords built : standard latlon (with latlon coord-system)
+        result = self.run_testcase(phenom_grid_mapping="grid: yco xco")
+        self.check_result(result)
+
+    def test_extended_mapping_basic_latlon_missing_coords(self):
+        # A basic reference example with a lat-long grid, but using extended
+        # grid mapping syntax on phenomenon.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(latitude)
+        #     004 : fc_provides_coordinate_(longitude)
+        #     005 : fc_build_coordinate_(latitude)
+        #     006 : fc_build_coordinate_(longitude)
+        # Notes:
+        #     * grid-mapping identified : regular latlon
+        #     * dim-coords identified : lat+lon
+        #     * coords built : standard latlon (with latlon coord-system)
+        result = self.run_testcase(
+            phenom_grid_mapping="grid: yco bad_coord",
+            warning_regex="Missing CF-netCDF coordinate variable 'bad_coord'",
+        )
+        self.check_result(result, xco_no_cs=True)
+
 
 class Test__aux_latlons(Mixin__grid_mapping, tests.IrisTest):
     # Testcases for translating auxiliary latitude+longitude variables
@@ -838,6 +878,55 @@ def test_aux_lat_rotated(self):
         )
         self.check_result(result, yco_is_aux=True, yco_no_cs=True)
 
+    def test_extended_grid_mapping_aux_lon(self):
+        # Change the name of xdim, and put xco on the coords list.
+        # Uses extended grid mapping syntax.
+        # In this case, the Aux coord WILL have a coordinate system
+        # as extended grid mapping allows for specification of
+        # explicit coordinate_systems for individual coordinate.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(latitude)
+        #     004 : fc_build_coordinate_(latitude)
+        #     005 : fc_build_auxiliary_coordinate_longitude
+        result = self.run_testcase(
+            xco_is_dim=False, phenom_grid_mapping="grid: yco xco"
+        )
+        self.check_result(result, xco_is_aux=True, xco_no_cs=False)
+
+    def test_extended_grid_mapping_aux_lat(self):
+        # As previous, but with the Y coord.
+        # Uses extended grid mapping syntax
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_provides_coordinate_(longitude)
+        #     004 : fc_build_coordinate_(longitude)
+        #     005 : fc_build_auxiliary_coordinate_latitude
+        result = self.run_testcase(
+            yco_is_dim=False, phenom_grid_mapping="grid: yco xco"
+        )
+        self.check_result(result, yco_is_aux=True, yco_no_cs=False)
+
+    def test_extended_grid_mapping_aux_lat_and_lon(self):
+        # Make *both* X and Y coords into aux-coords.
+        # Uses extended grid mapping syntax; allows coord_system
+        # to be added to an AuxCoord, so cube.coord_system() returns
+        # a valid coordinate_system.
+        #
+        # Rules Triggered:
+        #     001 : fc_default
+        #     002 : fc_provides_grid_mapping_(latitude_longitude)
+        #     003 : fc_build_auxiliary_coordinate_longitude
+        #     004 : fc_build_auxiliary_coordinate_latitude
+        result = self.run_testcase(
+            xco_is_dim=False, yco_is_dim=False, phenom_grid_mapping="grid: yco xco"
+        )
+        self.check_result(result, xco_is_aux=True, yco_is_aux=True, cube_no_cs=False)
+
 
 class Test__nondimcoords(Mixin__grid_mapping, tests.IrisTest):
     @classmethod