diff --git a/docs/src/further_topics/netcdf_io.rst b/docs/src/further_topics/netcdf_io.rst index 682918d5f4..47d85ffeab 100644 --- a/docs/src/further_topics/netcdf_io.rst +++ b/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 +`_ +, 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 + + : [] [: ...] + +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 ; + + + int crsWGS84 ; + crsWGS84:grid_mapping_name = "latitude_longitude" ; + crsWGS84:longitude_of_prime_meridian = 0. ; + + + +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 \ No newline at end of file diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst index 592b99b486..55daabbad2 100644 --- a/docs/src/whatsnew/latest.rst +++ b/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/ diff --git a/lib/iris/cube.py b/lib/iris/cube.py index d8006c1f25..04060160b9 100644 --- a/lib/iris/cube.py +++ b/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) diff --git a/lib/iris/exceptions.py b/lib/iris/exceptions.py index 450afce3a6..870a64901f 100644 --- a/lib/iris/exceptions.py +++ b/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 diff --git a/lib/iris/fileformats/_nc_load_rules/actions.py b/lib/iris/fileformats/_nc_load_rules/actions.py index ee19b54b4b..af84a84fa9 100644 --- a/lib/iris/fileformats/_nc_load_rules/actions.py +++ b/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) diff --git a/lib/iris/fileformats/_nc_load_rules/helpers.py b/lib/iris/fileformats/_nc_load_rules/helpers.py index fc47625943..35c2e96924 100644 --- a/lib/iris/fileformats/_nc_load_rules/helpers.py +++ b/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 ':' 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 followed by colon +# \w+ # Matches a +# )+ # 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+:"), + "`:` identifier followed immediately by another `:` identifier", + ), + ( + re.compile(r"\w+: *$"), + "`:` is empty - missing coordinate list", + ), + ( + re.compile(r"^\w+ +\w+"), + "Multiple coordinates found without `:` 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.""" diff --git a/lib/iris/fileformats/cf.py b/lib/iris/fileformats/cf.py index e69d686c0f..23186072a4 100644 --- a/lib/iris/fileformats/cf.py +++ b/lib/iris/fileformats/cf.py @@ -24,6 +24,8 @@ import numpy as np import numpy.ma as ma +import iris.exceptions +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 @@ -654,7 +656,9 @@ class CFGridMappingVariable(CFVariable): cf_identity = "grid_mapping" @classmethod - def identify(cls, variables, ignore=None, target=None, warn=True): + def identify( + cls, variables, ignore=None, target=None, warn=True, coord_system_mappings=None + ): result = {} ignore, target = cls._identify_common(variables, ignore, target) @@ -664,19 +668,61 @@ def identify(cls, variables, ignore=None, target=None, warn=True): nc_var_att = getattr(nc_var, cls.cf_identity, None) 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]) + # All `grid_mapping` attributes will already have been parsed prior + # to `identify` being called and passed in as an argument. We can + # ignore the attribute here (it's just used to identify that a grid + # mapping exists for this data variable) and get the pre-parsed + # mapping from the `coord_mapping_systems` keyword: + cs_mappings = None + if coord_system_mappings: + cs_mappings = coord_system_mappings.get(nc_var_name, None) + + if not cs_mappings: + # If cs_mappings is None, some parse error must have occurred and the + # user will have already been warned by `_parse_extended_grid_mappings` + continue + # group the cs_mappings by coordinate system, as we want to iterate over coord systems: + uniq_cs = set(cs_mappings.values()) + cs_coord_mappings = { + cs: [ + coord + for coord, coord_cs in cs_mappings.items() + if cs == coord_cs + ] + for cs in uniq_cs + } + + for name, coords in cs_coord_mappings.items(): + 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: + has_a_valid_coord = False + if coords: + for coord_name in coords: + # coord_name could be None if simple grid_mapping is used. + if coord_name is None or ( + coord_name and coord_name in variables + ): + has_a_valid_coord = True + else: + 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, + ) + # Only add as a CFGridMappingVariable if at least one of its referenced coords exists: + if has_a_valid_coord: + result[name] = CFGridMappingVariable( + name, variables[name] + ) return result @@ -1327,6 +1373,9 @@ def __init__(self, file_source, warn=False, monotonic=False): #: Collection of CF-netCDF variables associated with this netCDF file self.cf_group = self.CFGroup() + # Result of parsing "grid_mapping" attribute; mapping of coordinate_system => coordinates + self._coord_system_mappings = {} + # Issue load optimisation warning. if warn and self._dataset.file_format in [ "NETCDF3_CLASSIC", @@ -1395,6 +1444,18 @@ def _translate(self, variables): """Classify the netCDF variables into CF-netCDF variables.""" netcdf_variable_names = list(variables.keys()) + # Parse all instances of "grid_mapping" attributes and store in CFReader + # This avoids re-parsing the grid_mappings each time they are needed. + for nc_var in variables.values(): + if grid_mapping_attr := getattr(nc_var, "grid_mapping", None): + try: + cs_mappings = hh._parse_extended_grid_mapping(grid_mapping_attr) + self._coord_system_mappings[nc_var.name] = cs_mappings + except iris.exceptions.CFParseError as e: + msg = f"Error parsing `grid_mapping` attribute for {nc_var.name}: {str(e)}" + warnings.warn(msg, category=iris.warnings.IrisCfWarning) + continue + # Identify all CF coordinate variables first. This must be done # first as, by CF convention, the definition of a CF auxiliary # coordinate variable may include a scalar CF coordinate variable, @@ -1413,7 +1474,15 @@ def _translate(self, variables): if issubclass(variable_type, CFGridMappingVariable) else coordinate_names ) - self.cf_group.update(variable_type.identify(variables, ignore=ignore)) + kwargs = ( + {"coord_system_mappings": self._coord_system_mappings} + if issubclass(variable_type, CFGridMappingVariable) + else {} + ) + + self.cf_group.update( + variable_type.identify(variables, ignore=ignore, **kwargs) + ) # Identify global netCDF attributes. attr_dict = { @@ -1555,12 +1624,18 @@ def _span_check( # Build CF variable relationships. for variable_type in self._variable_types: ignore = [] + kwargs = {} # Avoid UGridAuxiliaryCoordinateVariables also being # processed as CFAuxiliaryCoordinateVariables. if not is_mesh_var: ignore += ugrid_coord_names # Prevent grid mapping variables being mis-identified as CF coordinate variables. - if not issubclass(variable_type, CFGridMappingVariable): + if issubclass(variable_type, CFGridMappingVariable): + # pass parsed grid_mappings to CFGridMappingVariable types + kwargs.update( + {"coord_system_mappings": self._coord_system_mappings} + ) + else: ignore += coordinate_names match = variable_type.identify( @@ -1568,6 +1643,7 @@ def _span_check( ignore=ignore, target=cf_variable.cf_name, warn=False, + **kwargs, ) # Sanity check dimensionality coverage. for cf_name in match: diff --git a/lib/iris/fileformats/netcdf/loader.py b/lib/iris/fileformats/netcdf/loader.py index b7895c4dda..216df67590 100644 --- a/lib/iris/fileformats/netcdf/loader.py +++ b/lib/iris/fileformats/netcdf/loader.py @@ -79,6 +79,12 @@ 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"] = {} + + # Add the parsed coordinate reference system mappings + engine.cube_parts["coordinate_system_mappings"] = cf._coord_system_mappings.get( + engine.cf_var.cf_name, None + ) # Assert facts for CF coordinates. for cf_name in cf_group.coordinates.keys(): @@ -445,6 +451,13 @@ def fix_attributes_all_elements(role_name): for method in cube.cell_methods ] + # Set extended_grid_mapping property ONLY if extended grid_mapping was used. + # This avoids having an unnecessary `iris_extended_grid_mapping` attribute entry. + if cs_mappings := engine.cube_parts.get("coordinate_system_mappings", None): + # `None` as a mapping key implies simple mapping syntax (single coord system) + if None not in cs_mappings: + cube.extended_grid_mapping = True + if DEBUG: # Show activation statistics for this data-var (i.e. cube). _actions_activation_stats(engine, cf_var.cf_name) diff --git a/lib/iris/fileformats/netcdf/saver.py b/lib/iris/fileformats/netcdf/saver.py index 95bd92d8a1..553bb1e199 100644 --- a/lib/iris/fileformats/netcdf/saver.py +++ b/lib/iris/fileformats/netcdf/saver.py @@ -104,7 +104,12 @@ _CF_GLOBAL_ATTRS = ["conventions", "featureType", "history", "title"] # UKMO specific attributes that should not be global. -_UKMO_DATA_ATTRS = ["STASH", "um_stash_source", "ukmo__process_flags"] +_UKMO_DATA_ATTRS = [ + "STASH", + "um_stash_source", + "ukmo__process_flags", + "iris_extended_grid_mapping", +] # TODO: whenever we advance to CF-1.11 we should then discuss a completion date # for the deprecation of Rotated Mercator in coord_systems.py and @@ -1934,6 +1939,214 @@ def _create_cf_cell_methods(self, cube, dimension_names): return " ".join(cell_methods) + def _add_grid_mapping_to_dataset(self, cs, extended_grid_mapping=False): + """Create a CF-netCDF grid mapping variable and add to the dataset. + + Parameters + ---------- + cs : :class:`iris.coord_system.CoordSystem` + The :class:`iris.coord_system.CoordSystem` used to generate the + new netCDF CF grid mapping variable. + + Returns + ------- + None + """ + cf_var_grid = self._dataset.createVariable(cs.grid_mapping_name, np.int32) + _setncattr(cf_var_grid, "grid_mapping_name", cs.grid_mapping_name) + + def add_ellipsoid(ellipsoid): + cf_var_grid.longitude_of_prime_meridian = ( + ellipsoid.longitude_of_prime_meridian + ) + semi_major = ellipsoid.semi_major_axis + semi_minor = ellipsoid.semi_minor_axis + if semi_minor == semi_major: + cf_var_grid.earth_radius = semi_major + else: + cf_var_grid.semi_major_axis = semi_major + cf_var_grid.semi_minor_axis = semi_minor + if ellipsoid.datum is not None: + cf_var_grid.horizontal_datum_name = ellipsoid.datum + + # latlon + if isinstance(cs, iris.coord_systems.GeogCS): + add_ellipsoid(cs) + + # rotated latlon + elif isinstance(cs, iris.coord_systems.RotatedGeogCS): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.grid_north_pole_latitude = cs.grid_north_pole_latitude + cf_var_grid.grid_north_pole_longitude = cs.grid_north_pole_longitude + cf_var_grid.north_pole_grid_longitude = cs.north_pole_grid_longitude + + # tmerc + elif isinstance(cs, iris.coord_systems.TransverseMercator): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_central_meridian = cs.longitude_of_central_meridian + cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + cf_var_grid.scale_factor_at_central_meridian = ( + cs.scale_factor_at_central_meridian + ) + + # merc + elif isinstance(cs, iris.coord_systems.Mercator): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_projection_origin = ( + cs.longitude_of_projection_origin + ) + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + # Only one of these should be set + if cs.standard_parallel is not None: + cf_var_grid.standard_parallel = cs.standard_parallel + elif cs.scale_factor_at_projection_origin is not None: + cf_var_grid.scale_factor_at_projection_origin = ( + cs.scale_factor_at_projection_origin + ) + + # lcc + elif isinstance(cs, iris.coord_systems.LambertConformal): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.standard_parallel = cs.secant_latitudes + cf_var_grid.latitude_of_projection_origin = cs.central_lat + cf_var_grid.longitude_of_central_meridian = cs.central_lon + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + + # polar stereo (have to do this before Stereographic because it subclasses it) + elif isinstance(cs, iris.coord_systems.PolarStereographic): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.latitude_of_projection_origin = cs.central_lat + cf_var_grid.straight_vertical_longitude_from_pole = cs.central_lon + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + # Only one of these should be set + if cs.true_scale_lat is not None: + cf_var_grid.true_scale_lat = cs.true_scale_lat + elif cs.scale_factor_at_projection_origin is not None: + cf_var_grid.scale_factor_at_projection_origin = ( + cs.scale_factor_at_projection_origin + ) + else: + cf_var_grid.scale_factor_at_projection_origin = 1.0 + + # stereo + elif isinstance(cs, iris.coord_systems.Stereographic): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_projection_origin = cs.central_lon + cf_var_grid.latitude_of_projection_origin = cs.central_lat + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + # Only one of these should be set + if cs.true_scale_lat is not None: + msg = ( + "It is not valid CF to save a true_scale_lat for " + "a Stereographic grid mapping." + ) + raise ValueError(msg) + elif cs.scale_factor_at_projection_origin is not None: + cf_var_grid.scale_factor_at_projection_origin = ( + cs.scale_factor_at_projection_origin + ) + else: + cf_var_grid.scale_factor_at_projection_origin = 1.0 + + # osgb (a specific tmerc) + elif isinstance(cs, iris.coord_systems.OSGB): + warnings.warn( + "OSGB coordinate system not yet handled", + category=iris.warnings.IrisSaveWarning, + ) + + # lambert azimuthal equal area + elif isinstance(cs, iris.coord_systems.LambertAzimuthalEqualArea): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_projection_origin = ( + cs.longitude_of_projection_origin + ) + cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + + # albers conical equal area + elif isinstance(cs, iris.coord_systems.AlbersEqualArea): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_central_meridian = cs.longitude_of_central_meridian + cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + cf_var_grid.standard_parallel = cs.standard_parallels + + # vertical perspective + elif isinstance(cs, iris.coord_systems.VerticalPerspective): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_projection_origin = ( + cs.longitude_of_projection_origin + ) + cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + cf_var_grid.perspective_point_height = cs.perspective_point_height + + # geostationary + elif isinstance(cs, iris.coord_systems.Geostationary): + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.longitude_of_projection_origin = ( + cs.longitude_of_projection_origin + ) + cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + cf_var_grid.perspective_point_height = cs.perspective_point_height + cf_var_grid.sweep_angle_axis = cs.sweep_angle_axis + + # oblique mercator (and rotated variant) + # Use duck-typing over isinstance() - subclasses (i.e. + # RotatedMercator) upset mock tests. + elif getattr(cs, "grid_mapping_name", None) == "oblique_mercator": + # RotatedMercator subclasses ObliqueMercator, and RM + # instances are implicitly saved as OM due to inherited + # properties. This is correct because CF 1.11 is removing + # all mention of RM. + if cs.ellipsoid: + add_ellipsoid(cs.ellipsoid) + cf_var_grid.azimuth_of_central_line = cs.azimuth_of_central_line + cf_var_grid.latitude_of_projection_origin = cs.latitude_of_projection_origin + cf_var_grid.longitude_of_projection_origin = ( + cs.longitude_of_projection_origin + ) + cf_var_grid.false_easting = cs.false_easting + cf_var_grid.false_northing = cs.false_northing + cf_var_grid.scale_factor_at_projection_origin = ( + cs.scale_factor_at_projection_origin + ) + + # other + else: + warnings.warn( + "Unable to represent the horizontal " + "coordinate system. The coordinate system " + "type %r is not yet implemented." % type(cs), + category=iris.warnings.IrisSaveWarning, + ) + + # add WKT string + if extended_grid_mapping: + cf_var_grid.crs_wkt = cs.as_cartopy_crs().to_wkt() + def _create_cf_grid_mapping(self, cube, cf_var_cube): """Create CF-netCDF grid mapping and associated CF-netCDF variable. @@ -1953,227 +2166,108 @@ def _create_cf_grid_mapping(self, cube, cf_var_cube): None """ - cs = cube.coord_system("CoordSystem") - if cs is not None: + coord_systems = [] + if cube.extended_grid_mapping: + # get unique list of all coord_systems on cube coords: + for coord in cube.coords(): + if coord.coord_system and coord.coord_system not in coord_systems: + coord_systems.append(coord.coord_system) + else: + # will only return a single coord system (if one exists): + if cs := cube.coord_system("CoordSystem"): + coord_systems.append(cs) + + grid_mappings = {} + matched_all_coords = True + + for cs in coord_systems: # Grid var not yet created? if cs not in self._coord_systems: + # handle potential duplicate netCDF variable names: while cs.grid_mapping_name in self._dataset.variables: aname = self._increment_name(cs.grid_mapping_name) cs.grid_mapping_name = aname - cf_var_grid = self._dataset.createVariable( - cs.grid_mapping_name, np.int32 + # create grid mapping variable on dataset for this coordinate system: + self._add_grid_mapping_to_dataset( + cs, extended_grid_mapping=cube.extended_grid_mapping ) - _setncattr(cf_var_grid, "grid_mapping_name", cs.grid_mapping_name) - - def add_ellipsoid(ellipsoid): - cf_var_grid.longitude_of_prime_meridian = ( - ellipsoid.longitude_of_prime_meridian - ) - semi_major = ellipsoid.semi_major_axis - semi_minor = ellipsoid.semi_minor_axis - if semi_minor == semi_major: - cf_var_grid.earth_radius = semi_major - else: - cf_var_grid.semi_major_axis = semi_major - cf_var_grid.semi_minor_axis = semi_minor - if ellipsoid.datum is not None: - cf_var_grid.horizontal_datum_name = ellipsoid.datum - - # latlon - if isinstance(cs, iris.coord_systems.GeogCS): - add_ellipsoid(cs) - - # rotated latlon - elif isinstance(cs, iris.coord_systems.RotatedGeogCS): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.grid_north_pole_latitude = cs.grid_north_pole_latitude - cf_var_grid.grid_north_pole_longitude = cs.grid_north_pole_longitude - cf_var_grid.north_pole_grid_longitude = cs.north_pole_grid_longitude - - # tmerc - elif isinstance(cs, iris.coord_systems.TransverseMercator): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_central_meridian = ( - cs.longitude_of_central_meridian - ) - cf_var_grid.latitude_of_projection_origin = ( - cs.latitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - cf_var_grid.scale_factor_at_central_meridian = ( - cs.scale_factor_at_central_meridian - ) - - # merc - elif isinstance(cs, iris.coord_systems.Mercator): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_projection_origin = ( - cs.longitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - # Only one of these should be set - if cs.standard_parallel is not None: - cf_var_grid.standard_parallel = cs.standard_parallel - elif cs.scale_factor_at_projection_origin is not None: - cf_var_grid.scale_factor_at_projection_origin = ( - cs.scale_factor_at_projection_origin - ) + self._coord_systems.append(cs) - # lcc - elif isinstance(cs, iris.coord_systems.LambertConformal): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.standard_parallel = cs.secant_latitudes - cf_var_grid.latitude_of_projection_origin = cs.central_lat - cf_var_grid.longitude_of_central_meridian = cs.central_lon - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - - # polar stereo (have to do this before Stereographic because it subclasses it) - elif isinstance(cs, iris.coord_systems.PolarStereographic): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.latitude_of_projection_origin = cs.central_lat - cf_var_grid.straight_vertical_longitude_from_pole = cs.central_lon - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - # Only one of these should be set - if cs.true_scale_lat is not None: - cf_var_grid.true_scale_lat = cs.true_scale_lat - elif cs.scale_factor_at_projection_origin is not None: - cf_var_grid.scale_factor_at_projection_origin = ( - cs.scale_factor_at_projection_origin - ) - else: - cf_var_grid.scale_factor_at_projection_origin = 1.0 - - # stereo - elif isinstance(cs, iris.coord_systems.Stereographic): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_projection_origin = cs.central_lon - cf_var_grid.latitude_of_projection_origin = cs.central_lat - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - # Only one of these should be set - if cs.true_scale_lat is not None: + # create the `grid_mapping` attribute for the data variable: + if cube.extended_grid_mapping: + # Order the coordinates as per the order in the CRS/WKT string. + # (We should only ever have a coordinate system for horizontal + # spatial coords, so check for east/north directions) + ordered_coords = [] + for ax_info in cs.as_cartopy_crs().axis_info: + try: + match ax_info.direction.lower(): + case "east": + ordered_coords.append( + cube.coord(axis="X", coord_system=cs) + ) + case "north": + ordered_coords.append( + cube.coord(axis="Y", coord_system=cs) + ) + case _: + msg = ( + f"Can't handle axis direction {ax_info.direction!r}" + ) + raise NotImplementedError(msg) + except iris.exceptions.CoordinateNotFoundError as e: msg = ( - "It is not valid CF to save a true_scale_lat for " - "a Stereographic grid mapping." + f"Failed to assign coordinate for {ax_info.name} axis of " + f"coordinate system {cs.grid_mapping_name}: {str(e)}" ) - raise ValueError(msg) - elif cs.scale_factor_at_projection_origin is not None: - cf_var_grid.scale_factor_at_projection_origin = ( - cs.scale_factor_at_projection_origin + warnings.warn(msg, category=iris.warnings.IrisSaveWarning) + + # fall back to simple grid mapping (single crs entry for DimCoord only) + matched_all_coords = False + + # Get variable names, and store them if necessary: + cfvar_names = [] + for coord in ordered_coords: + cfvar = self._name_coord_map.name(coord) + if not cfvar: + # not found - create and store it: + cfvar = self._get_coord_variable_name(cube, coord) + self._name_coord_map.append( + cfvar, self._get_coord_variable_name(cube, coord) ) - else: - cf_var_grid.scale_factor_at_projection_origin = 1.0 - - # osgb (a specific tmerc) - elif isinstance(cs, iris.coord_systems.OSGB): - warnings.warn( - "OSGB coordinate system not yet handled", - category=iris.warnings.IrisSaveWarning, - ) - - # lambert azimuthal equal area - elif isinstance(cs, iris.coord_systems.LambertAzimuthalEqualArea): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_projection_origin = ( - cs.longitude_of_projection_origin - ) - cf_var_grid.latitude_of_projection_origin = ( - cs.latitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - - # albers conical equal area - elif isinstance(cs, iris.coord_systems.AlbersEqualArea): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_central_meridian = ( - cs.longitude_of_central_meridian - ) - cf_var_grid.latitude_of_projection_origin = ( - cs.latitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - cf_var_grid.standard_parallel = cs.standard_parallels - - # vertical perspective - elif isinstance(cs, iris.coord_systems.VerticalPerspective): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_projection_origin = ( - cs.longitude_of_projection_origin - ) - cf_var_grid.latitude_of_projection_origin = ( - cs.latitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - cf_var_grid.perspective_point_height = cs.perspective_point_height - - # geostationary - elif isinstance(cs, iris.coord_systems.Geostationary): - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.longitude_of_projection_origin = ( - cs.longitude_of_projection_origin + cfvar_names.append(cfvar) + + coord_string = " ".join(cfvar_names) + grid_mappings[cs.grid_mapping_name] = coord_string + + # Refer to grid var + if len(coord_systems): + grid_mapping = None + if cube.extended_grid_mapping: + if matched_all_coords: + grid_mapping = " ".join( + f"{cs_name}: {cs_coords}" + for cs_name, cs_coords in grid_mappings.items() ) - cf_var_grid.latitude_of_projection_origin = ( - cs.latitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - cf_var_grid.perspective_point_height = cs.perspective_point_height - cf_var_grid.sweep_angle_axis = cs.sweep_angle_axis - - # oblique mercator (and rotated variant) - # Use duck-typing over isinstance() - subclasses (i.e. - # RotatedMercator) upset mock tests. - elif getattr(cs, "grid_mapping_name", None) == "oblique_mercator": - # RotatedMercator subclasses ObliqueMercator, and RM - # instances are implicitly saved as OM due to inherited - # properties. This is correct because CF 1.11 is removing - # all mention of RM. - if cs.ellipsoid: - add_ellipsoid(cs.ellipsoid) - cf_var_grid.azimuth_of_central_line = cs.azimuth_of_central_line - cf_var_grid.latitude_of_projection_origin = ( - cs.latitude_of_projection_origin - ) - cf_var_grid.longitude_of_projection_origin = ( - cs.longitude_of_projection_origin - ) - cf_var_grid.false_easting = cs.false_easting - cf_var_grid.false_northing = cs.false_northing - cf_var_grid.scale_factor_at_projection_origin = ( - cs.scale_factor_at_projection_origin - ) - - # other else: - warnings.warn( - "Unable to represent the horizontal " - "coordinate system. The coordinate system " - "type %r is not yet implemented." % type(cs), - category=iris.warnings.IrisSaveWarning, - ) - - self._coord_systems.append(cs) + # We didn't match all coords required for grid mapping, default to + # single coordinate system for DimCoords, if set. + for ax in ["x", "y"]: + try: + if dim_cs := cube.coord( + axis=ax, dim_coords=True + ).coord_system: + grid_mapping = dim_cs.grid_mapping_name + break + except iris.exceptions.CoordinateNotFoundError: + pass + else: + # Not using extended_grid_mapping: only ever be on cs in this case: + grid_mapping = coord_systems[0].grid_mapping_name - # Refer to grid var - _setncattr(cf_var_cube, "grid_mapping", cs.grid_mapping_name) + if grid_mapping: + _setncattr(cf_var_cube, "grid_mapping", grid_mapping) def _create_cf_data_variable( self, @@ -2311,8 +2405,8 @@ def set_packing_ncattrs(cfvar): attr_names = set(cube.attributes).intersection(local_keys) for attr_name in sorted(attr_names): - # Do not output 'conventions' attribute. - if attr_name.lower() == "conventions": + # Do not output 'conventions' or extended grid mapping attribute. + if attr_name.lower() in ["conventions", "iris_extended_grid_mapping"]: continue value = cube.attributes[attr_name] diff --git a/lib/iris/tests/integration/netcdf/test_coord_systems.py b/lib/iris/tests/integration/netcdf/test_coord_systems.py index 6c101d9024..9f4f272dd7 100644 --- a/lib/iris/tests/integration/netcdf/test_coord_systems.py +++ b/lib/iris/tests/integration/netcdf/test_coord_systems.py @@ -8,38 +8,23 @@ # importing anything else. import iris.tests as tests # isort:skip -from os.path import join as path_join -import shutil -import tempfile import warnings +import numpy as np import pytest import iris from iris.coords import DimCoord from iris.cube import Cube from iris.tests import stock as stock +from iris.tests._shared_utils import assert_CML from iris.tests.stock.netcdf import ncgen_from_cdl from iris.tests.unit.fileformats.netcdf.loader import test_load_cubes as tlc -@tests.skip_data -class TestCoordSystem(tests.IrisTest): - def setUp(self): - tlc.setUpModule() - - def tearDown(self): - tlc.tearDownModule() - - def test_load_laea_grid(self): - cube = iris.load_cube( - tests.get_data_path( - ("NetCDF", "lambert_azimuthal_equal_area", "euro_air_temp.nc") - ) - ) - self.assertCML(cube, ("netcdf", "netcdf_laea.cml")) - - datum_cf_var_cdl = """ +@pytest.fixture +def datum_cf_var_cdl(): + return """ netcdf output { dimensions: y = 4 ; @@ -84,7 +69,10 @@ def test_load_laea_grid(self): } """ - datum_wkt_cdl = """ + +@pytest.fixture +def datum_wkt_cdl(): + return """ netcdf output5 { dimensions: y = 4 ; @@ -93,7 +81,7 @@ def test_load_laea_grid(self): float data(y, x) ; data :standard_name = "toa_brightness_temperature" ; data :units = "K" ; - data :grid_mapping = "mercator" ; + data :grid_mapping = "mercator: x y" ; int mercator ; mercator:grid_mapping_name = "mercator" ; mercator:longitude_of_prime_meridian = 0. ; @@ -131,25 +119,94 @@ def test_load_laea_grid(self): } """ - def test_load_datum_wkt(self): + +@pytest.fixture +def multi_cs_osgb_wkt(): + return """ +netcdf osgb { +dimensions: + y = 5 ; + x = 4 ; +variables: + double x(x) ; + x:standard_name = "projection_x_coordinate" ; + x:long_name = "Easting" ; + x:units = "m" ; + double y(y) ; + y:standard_name = "projection_y_coordinate" ; + y:long_name = "Northing" ; + 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" ; + float temp(y, x) ; + temp:standard_name = "air_temperature" ; + temp:units = "K" ; + temp:coordinates = "lat lon" ; + temp:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ; + int crsOSGB ; + crsOSGB:grid_mapping_name = "transverse_mercator" ; + crsOSGB:semi_major_axis = 6377563.396 ; + crsOSGB:inverse_flattening = 299.3249646 ; + crsOSGB:longitude_of_prime_meridian = 0. ; + crsOSGB:latitude_of_projection_origin = 49. ; + crsOSGB:longitude_of_central_meridian = -2. ; + crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ; + crsOSGB:false_easting = 400000. ; + crsOSGB:false_northing = -100000. ; + crsOSGB:unit = "metre" ; + crsOSGB:crs_wkt = "PROJCRS[\\"unknown\\",BASEGEOGCRS[\\"unknown\\",DATUM[\\"Unknown based on Airy 1830 ellipsoid\\",ELLIPSOID[\\"Airy 1830\\",6377563.396,299.324964600004,LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]],PRIMEM[\\"Greenwich\\",0,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8901]]],CONVERSION[\\"unknown\\",METHOD[\\"Transverse Mercator\\",ID[\\"EPSG\\",9807]],PARAMETER[\\"Latitude of natural origin\\",49,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8801]],PARAMETER[\\"Longitude of natural origin\\",-2,ANGLEUNIT[\\"degree\\",0.0174532925199433],ID[\\"EPSG\\",8802]],PARAMETER[\\"Scale factor at natural origin\\",0.9996012717,SCALEUNIT[\\"unity\\",1],ID[\\"EPSG\\",8805]],PARAMETER[\\"False easting\\",400000,LENGTHUNIT[\\"metre\\",1],ID[\\"EPSG\\",8806]],PARAMETER[\\"False northing\\",-100000,LENGTHUNIT[\\"metre\\",1],ID[\\"EPSG\\",8807]]],CS[Cartesian,2],AXIS[\\"(E)\\",east,ORDER[1],LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]],AXIS[\\"(N)\\",north,ORDER[2],LENGTHUNIT[\\"metre\\",1,ID[\\"EPSG\\",9001]]]]" ; + int crsWGS84 ; + crsWGS84:grid_mapping_name = "latitude_longitude" ; + crsWGS84:longitude_of_prime_meridian = 0. ; + crsWGS84:semi_major_axis = 6378137. ; + crsWGS84:inverse_flattening = 298.257223563 ; + crsWGS84: crs_wkt = "GEOGCRS[\\"unknown\\",DATUM[\\"Unknown based on WGS 84 ellipsoid\\",ELLIPSOID[\\"WGS 84\\",6378137,298.257223562997,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]]]]" ; +data: + x = 1,2,3,4 ; + y = 1,2,3,4,5 ; +} + """ + + +@tests.skip_data +class TestCoordSystem: + @pytest.fixture(autouse=True) + def _setup(self): + tlc.setUpModule() + yield + tlc.tearDownModule() + + def test_load_laea_grid(self, request): + cube = iris.load_cube( + tests.get_data_path( + ("NetCDF", "lambert_azimuthal_equal_area", "euro_air_temp.nc") + ) + ) + assert_CML(request, cube, ("netcdf", "netcdf_laea.cml")) + + def test_load_datum_wkt(self, datum_wkt_cdl): expected = "OSGB 1936" - nc_path = tlc.cdl_to_nc(self.datum_wkt_cdl) + nc_path = tlc.cdl_to_nc(datum_wkt_cdl) with iris.FUTURE.context(datum_support=True): cube = iris.load_cube(nc_path) test_crs = cube.coord("projection_y_coordinate").coord_system actual = str(test_crs.as_cartopy_crs().datum) assert actual == expected - def test_no_load_datum_wkt(self): - nc_path = tlc.cdl_to_nc(self.datum_wkt_cdl) + def test_no_load_datum_wkt(self, datum_wkt_cdl): + nc_path = tlc.cdl_to_nc(datum_wkt_cdl) with pytest.warns(FutureWarning, match="iris.FUTURE.datum_support"): cube = iris.load_cube(nc_path) test_crs = cube.coord("projection_y_coordinate").coord_system actual = str(test_crs.as_cartopy_crs().datum) assert actual == "unknown" - def test_no_datum_no_warn(self): - new_cdl = self.datum_wkt_cdl.splitlines() + def test_no_datum_no_warn(self, datum_wkt_cdl): + new_cdl = datum_wkt_cdl.splitlines() new_cdl = [line for line in new_cdl if "DATUM" not in line] new_cdl = "\n".join(new_cdl) nc_path = tlc.cdl_to_nc(new_cdl) @@ -158,24 +215,40 @@ def test_no_datum_no_warn(self): warnings.simplefilter("error", FutureWarning) _ = iris.load_cube(nc_path) - def test_load_datum_cf_var(self): + def test_load_datum_cf_var(self, datum_cf_var_cdl): expected = "OSGB 1936" - nc_path = tlc.cdl_to_nc(self.datum_cf_var_cdl) + nc_path = tlc.cdl_to_nc(datum_cf_var_cdl) with iris.FUTURE.context(datum_support=True): cube = iris.load_cube(nc_path) test_crs = cube.coord("projection_y_coordinate").coord_system actual = str(test_crs.as_cartopy_crs().datum) assert actual == expected - def test_no_load_datum_cf_var(self): - nc_path = tlc.cdl_to_nc(self.datum_cf_var_cdl) + def test_no_load_datum_cf_var(self, datum_cf_var_cdl): + nc_path = tlc.cdl_to_nc(datum_cf_var_cdl) with pytest.warns(FutureWarning, match="iris.FUTURE.datum_support"): cube = iris.load_cube(nc_path) test_crs = cube.coord("projection_y_coordinate").coord_system actual = str(test_crs.as_cartopy_crs().datum) assert actual == "unknown" - def test_save_datum(self): + def test_load_multi_cs_wkt(self, multi_cs_osgb_wkt): + nc_path = tlc.cdl_to_nc(multi_cs_osgb_wkt) + with iris.FUTURE.context(datum_support=True): + cube = iris.load_cube(nc_path) + + assert len(cube.coord_systems()) == 2 + for name in ["projection_y_coordinate", "projection_y_coordinate"]: + assert ( + cube.coord(name).coord_system.grid_mapping_name == "transverse_mercator" + ) + for name in ["latitude", "longitude"]: + assert ( + cube.coord(name).coord_system.grid_mapping_name == "latitude_longitude" + ) + assert cube.extended_grid_mapping is True + + def test_save_datum(self, tmp_path): expected = "OSGB 1936" saved_crs = iris.coord_systems.Mercator( ellipsoid=iris.coord_systems.GeogCS.from_datum("OSGB36") @@ -204,24 +277,73 @@ def test_save_datum(self): (test_lon_coord, 2), ), ) - - with self.temp_filename(suffix=".nc") as filename: - iris.save(test_cube, filename) - with iris.FUTURE.context(datum_support=True): - cube = iris.load_cube(filename) + filename = tmp_path / "output.nc" + iris.save(test_cube, filename) + with iris.FUTURE.context(datum_support=True): + cube = iris.load_cube(filename) test_crs = cube.coord("projection_y_coordinate").coord_system actual = str(test_crs.as_cartopy_crs().datum) assert actual == expected + def test_save_multi_cs_wkt(self, tmp_path): + crsOSGB = iris.coord_systems.OSGB() + crsLatLon = iris.coord_systems.GeogCS(6e6) -class TestLoadMinimalGeostationary(tests.IrisTest): - """Check we can load data with a geostationary grid-mapping, even when the - 'false-easting' and 'false_northing' properties are missing. + dimx_coord = iris.coords.DimCoord( + np.arange(4), "projection_x_coordinate", coord_system=crsOSGB + ) + dimy_coord = iris.coords.DimCoord( + np.arange(5), "projection_y_coordinate", coord_system=crsOSGB + ) - """ + auxlon_coord = iris.coords.AuxCoord( + np.arange(20).reshape((5, 4)), + standard_name="longitude", + coord_system=crsLatLon, + ) + auxlat_coord = iris.coords.AuxCoord( + np.arange(20).reshape((5, 4)), + standard_name="latitude", + coord_system=crsLatLon, + ) + + test_cube = Cube( + np.ones(20).reshape((5, 4)), + standard_name="air_pressure", + units="Pa", + dim_coords_and_dims=( + (dimy_coord, 0), + (dimx_coord, 1), + ), + aux_coords_and_dims=( + (auxlat_coord, (0, 1)), + (auxlon_coord, (0, 1)), + ), + ) - _geostationary_problem_cdl = """ + test_cube.extended_grid_mapping = True + + filename = tmp_path / "output.nc" + iris.save(test_cube, filename) + with iris.FUTURE.context(datum_support=True): + cube = iris.load_cube(filename) + + assert len(cube.coord_systems()) == 2 + for name in ["projection_y_coordinate", "projection_y_coordinate"]: + assert ( + cube.coord(name).coord_system.grid_mapping_name == "transverse_mercator" + ) + for name in ["latitude", "longitude"]: + assert ( + cube.coord(name).coord_system.grid_mapping_name == "latitude_longitude" + ) + assert cube.extended_grid_mapping is True + + +@pytest.fixture(scope="module") +def geostationary_problem_cdl(): + return """ netcdf geostationary_problem_case { dimensions: y = 2 ; @@ -258,35 +380,34 @@ class TestLoadMinimalGeostationary(tests.IrisTest): x = 0, 1, 2 ; } -""" - - @classmethod - def setUpClass(cls): - # Create a temp directory for transient test files. - cls.temp_dir = tempfile.mkdtemp() - cls.path_test_cdl = path_join(cls.temp_dir, "geos_problem.cdl") - cls.path_test_nc = path_join(cls.temp_dir, "geos_problem.nc") - # Create reference CDL and netcdf files from the CDL text. + """ + + +class TestLoadMinimalGeostationary: + """Check we can load data with a geostationary grid-mapping, even when the + 'false-easting' and 'false_northing' properties are missing. + + """ + + @pytest.fixture(scope="class") + def geostationary_problem_ncfile(self, tmp_path_factory, geostationary_problem_cdl): + tmp_path = tmp_path_factory.mktemp("geos") + cdl_path = tmp_path / "geos_problem.cdl" + nc_path = tmp_path / "geos_problem.nc" ncgen_from_cdl( - cdl_str=cls._geostationary_problem_cdl, - cdl_path=cls.path_test_cdl, - nc_path=cls.path_test_nc, + cdl_str=geostationary_problem_cdl, + cdl_path=cdl_path, + nc_path=nc_path, ) + return nc_path - @classmethod - def tearDownClass(cls): - # Destroy the temp directory. - shutil.rmtree(cls.temp_dir) - - def test_geostationary_no_false_offsets(self): + def test_geostationary_no_false_offsets( + self, tmp_path, geostationary_problem_ncfile + ): # Check we can load the test data and coordinate system properties are correct. - cube = iris.load_cube(self.path_test_nc) + cube = iris.load_cube(geostationary_problem_ncfile) # Check the coordinate system properties has the correct default properties. cs = cube.coord_system() assert isinstance(cs, iris.coord_systems.Geostationary) assert cs.false_easting == 0.0 assert cs.false_northing == 0.0 - - -if __name__ == "__main__": - tests.main() diff --git a/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py b/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py index b89477bfe9..a2d6e9c754 100644 --- a/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py +++ b/lib/iris/tests/integration/test_netcdf__loadsaveattrs.py @@ -69,6 +69,9 @@ def global_attr(request): iris.fileformats.netcdf.saver._CF_DATA_ATTRS + iris.fileformats.netcdf.saver._UKMO_DATA_ATTRS ) +# Don't test iris_extended_grid_mapping, as it is a special attribute that +# is not expected to always roundtrip. +_LOCAL_TEST_ATTRS = [a for a in _LOCAL_TEST_ATTRS if a != "iris_extended_grid_mapping"] # Define a fixture to parametrise over the 'local-style' test attributes. @@ -513,6 +516,7 @@ def fetch_results( "standard_error_multiplier", "STASH", "um_stash_source", + "iris_extended_grid_mapping", ] _MATRIX_ATTRNAMES = [attr for attr in _MATRIX_ATTRNAMES if attr not in _SPECIAL_ATTRS] diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs.cdl new file mode 100644 index 0000000000..547813fdf1 --- /dev/null +++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs.cdl @@ -0,0 +1,40 @@ +dimensions: + projection_x_coordinate = 4 ; + projection_y_coordinate = 3 ; +variables: + int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ; + air_pressure_anomaly:standard_name = "air_pressure_anomaly" ; + air_pressure_anomaly:grid_mapping = "transverse_mercator: projection_x_coordinate projection_y_coordinate latitude_longitude: longitude latitude" ; + air_pressure_anomaly:coordinates = "latitude longitude" ; + int transverse_mercator ; + transverse_mercator:grid_mapping_name = "transverse_mercator" ; + transverse_mercator:longitude_of_prime_meridian = 0. ; + transverse_mercator:semi_major_axis = 6377563.396 ; + transverse_mercator:semi_minor_axis = 6356256.909 ; + transverse_mercator:longitude_of_central_meridian = -2. ; + transverse_mercator:latitude_of_projection_origin = 49. ; + transverse_mercator:false_easting = -400000. ; + transverse_mercator:false_northing = 100000. ; + transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ; + transverse_mercator:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",49,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-2,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996012717,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",-400000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",100000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ; + int latitude_longitude ; + latitude_longitude:grid_mapping_name = "latitude_longitude" ; + latitude_longitude:longitude_of_prime_meridian = 0. ; + latitude_longitude:semi_major_axis = 6377563.396 ; + latitude_longitude:semi_minor_axis = 6356256.909 ; + latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,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]]]]" ; + int64 projection_y_coordinate(projection_y_coordinate) ; + projection_y_coordinate:axis = "Y" ; + projection_y_coordinate:units = "m" ; + projection_y_coordinate:standard_name = "projection_y_coordinate" ; + int64 projection_x_coordinate(projection_x_coordinate) ; + projection_x_coordinate:axis = "X" ; + projection_x_coordinate:units = "m" ; + projection_x_coordinate:standard_name = "projection_x_coordinate" ; + int64 latitude(projection_y_coordinate, projection_x_coordinate) ; + latitude:units = "degrees_north" ; + latitude:standard_name = "latitude" ; + int64 longitude(projection_y_coordinate, projection_x_coordinate) ; + longitude:units = "degrees_east" ; + longitude:standard_name = "longitude" ; +} diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs_missing_coord.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs_missing_coord.cdl new file mode 100644 index 0000000000..34c7061c9b --- /dev/null +++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/multi_cs_missing_coord.cdl @@ -0,0 +1,37 @@ +dimensions: + projection_x_coordinate = 4 ; + projection_y_coordinate = 3 ; +variables: + int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ; + air_pressure_anomaly:standard_name = "air_pressure_anomaly" ; + air_pressure_anomaly:grid_mapping = "transverse_mercator" ; + air_pressure_anomaly:coordinates = "longitude" ; + int transverse_mercator ; + transverse_mercator:grid_mapping_name = "transverse_mercator" ; + transverse_mercator:longitude_of_prime_meridian = 0. ; + transverse_mercator:semi_major_axis = 6377563.396 ; + transverse_mercator:semi_minor_axis = 6356256.909 ; + transverse_mercator:longitude_of_central_meridian = -2. ; + transverse_mercator:latitude_of_projection_origin = 49. ; + transverse_mercator:false_easting = -400000. ; + transverse_mercator:false_northing = 100000. ; + transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ; + transverse_mercator:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",49,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-2,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996012717,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",-400000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",100000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ; + int latitude_longitude ; + latitude_longitude:grid_mapping_name = "latitude_longitude" ; + latitude_longitude:longitude_of_prime_meridian = 0. ; + latitude_longitude:semi_major_axis = 6377563.396 ; + latitude_longitude:semi_minor_axis = 6356256.909 ; + latitude_longitude:crs_wkt = "GEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,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]]]]" ; + int64 projection_y_coordinate(projection_y_coordinate) ; + projection_y_coordinate:axis = "Y" ; + projection_y_coordinate:units = "m" ; + projection_y_coordinate:standard_name = "projection_y_coordinate" ; + int64 projection_x_coordinate(projection_x_coordinate) ; + projection_x_coordinate:axis = "X" ; + projection_x_coordinate:units = "m" ; + projection_x_coordinate:standard_name = "projection_x_coordinate" ; + int64 longitude(projection_y_coordinate, projection_x_coordinate) ; + longitude:units = "degrees_east" ; + longitude:standard_name = "longitude" ; +} diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_aux_cs.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_aux_cs.cdl new file mode 100644 index 0000000000..9b5aeb2dfb --- /dev/null +++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_aux_cs.cdl @@ -0,0 +1,34 @@ +dimensions: + projection_x_coordinate = 4 ; + projection_y_coordinate = 3 ; +variables: + int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ; + air_pressure_anomaly:standard_name = "air_pressure_anomaly" ; + air_pressure_anomaly:grid_mapping = "transverse_mercator: projection_x_coordinate projection_y_coordinate" ; + air_pressure_anomaly:coordinates = "latitude longitude" ; + int transverse_mercator ; + transverse_mercator:grid_mapping_name = "transverse_mercator" ; + transverse_mercator:longitude_of_prime_meridian = 0. ; + transverse_mercator:semi_major_axis = 6377563.396 ; + transverse_mercator:semi_minor_axis = 6356256.909 ; + transverse_mercator:longitude_of_central_meridian = -2. ; + transverse_mercator:latitude_of_projection_origin = 49. ; + transverse_mercator:false_easting = -400000. ; + transverse_mercator:false_northing = 100000. ; + transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ; + transverse_mercator:crs_wkt = "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown\",ELLIPSOID[\"unknown\",6377563.396,299.324961266495,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",49,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-2,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996012717,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",-400000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",100000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" ; + int64 projection_y_coordinate(projection_y_coordinate) ; + projection_y_coordinate:axis = "Y" ; + projection_y_coordinate:units = "m" ; + projection_y_coordinate:standard_name = "projection_y_coordinate" ; + int64 projection_x_coordinate(projection_x_coordinate) ; + projection_x_coordinate:axis = "X" ; + projection_x_coordinate:units = "m" ; + projection_x_coordinate:standard_name = "projection_x_coordinate" ; + int64 latitude(projection_y_coordinate, projection_x_coordinate) ; + latitude:units = "degrees_north" ; + latitude:standard_name = "latitude" ; + int64 longitude(projection_y_coordinate, projection_x_coordinate) ; + longitude:units = "degrees_east" ; + longitude:standard_name = "longitude" ; +} diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_cs.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_cs.cdl new file mode 100644 index 0000000000..ce8825b16d --- /dev/null +++ b/lib/iris/tests/results/unit/fileformats/netcdf/saver/Saver/write_extended_grid_mapping/no_cs.cdl @@ -0,0 +1,22 @@ +dimensions: + projection_x_coordinate = 4 ; + projection_y_coordinate = 3 ; +variables: + int64 air_pressure_anomaly(projection_y_coordinate, projection_x_coordinate) ; + air_pressure_anomaly:standard_name = "air_pressure_anomaly" ; + air_pressure_anomaly:coordinates = "latitude longitude" ; + int64 projection_y_coordinate(projection_y_coordinate) ; + projection_y_coordinate:axis = "Y" ; + projection_y_coordinate:units = "m" ; + projection_y_coordinate:standard_name = "projection_y_coordinate" ; + int64 projection_x_coordinate(projection_x_coordinate) ; + projection_x_coordinate:axis = "X" ; + projection_x_coordinate:units = "m" ; + projection_x_coordinate:standard_name = "projection_x_coordinate" ; + int64 latitude(projection_y_coordinate, projection_x_coordinate) ; + latitude:units = "degrees_north" ; + latitude:standard_name = "latitude" ; + int64 longitude(projection_y_coordinate, projection_x_coordinate) ; + longitude:units = "degrees_east" ; + longitude:standard_name = "longitude" ; +} diff --git a/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py b/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py index d80b4311c2..95069ba378 100644 --- a/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py +++ b/lib/iris/tests/unit/coord_systems/test_ObliqueMercator.py @@ -5,7 +5,6 @@ """Unit tests for the :class:`iris.coord_systems.ObliqueMercator` class.""" from typing import List, NamedTuple -from unittest.mock import Mock from cartopy import crs as ccrs import pytest @@ -130,18 +129,20 @@ def make_instance(self) -> ObliqueMercator: def instance(self): return self.make_instance() + @pytest.fixture() + def mock_ccrs(self, mocker): + return mocker.patch("cartopy.crs.ObliqueMercator", autospec=True) + def test_instantiate(self): _ = self.make_instance() - def test_cartopy_crs(self, instance): - ccrs.ObliqueMercator = Mock() + def test_cartopy_crs(self, instance, mock_ccrs): instance.as_cartopy_crs() - ccrs.ObliqueMercator.assert_called_with(**self.cartopy_kwargs_expected) + mock_ccrs.assert_called_with(**self.cartopy_kwargs_expected) - def test_cartopy_projection(self, instance): - ccrs.ObliqueMercator = Mock() + def test_cartopy_projection(self, instance, mock_ccrs): instance.as_cartopy_projection() - ccrs.ObliqueMercator.assert_called_with(**self.cartopy_kwargs_expected) + mock_ccrs.assert_called_with(**self.cartopy_kwargs_expected) @pytest.fixture() def label_class(self, instance): diff --git a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py index 7c79740023..296765f853 100644 --- a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py +++ b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/__init__.py @@ -65,7 +65,7 @@ def tearDownClass(cls): # Destroy a temp directory for temp files. shutil.rmtree(cls.temp_dirpath) - def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path): + def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path, mocker=None): """Load the 'phenom' data variable in a CDL testcase, as a cube. Using ncgen, CFReader and the _load_cube call. @@ -85,8 +85,12 @@ def load_cube_from_cdl(self, cdl_string, cdl_path, nc_path): # If debug enabled, switch on the activation summary debug output. # Use 'patch' so it is restored after the test. - self.patch("iris.fileformats.netcdf.loader.DEBUG", self.debug_info) - + if mocker: + # If pytest mocker provided, use that to patch: + mocker.patch("iris.fileformats.netcdf.loader.DEBUG", self.debug_info) + elif hasattr(self, "patch"): + # keep old UnitTest patch working + self.patch("iris.fileformats.netcdf.loader.DEBUG", self.debug_info) with warnings.catch_warnings(): warnings.filterwarnings( "ignore", diff --git a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py index 8731049ea8..202fb0fa16 100644 --- a/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py +++ b/lib/iris/tests/unit/fileformats/nc_load_rules/actions/test__grid_mappings.py @@ -9,8 +9,12 @@ """ +import iris.coord_systems + import iris.tests as tests # isort: skip +import pytest + import iris.coord_systems as ics import iris.fileformats._nc_load_rules.helpers as hh from iris.loading import LOAD_PROBLEMS @@ -27,6 +31,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 +231,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" ; @@ -325,6 +330,7 @@ def check_result( self.assertIsNone(yco_cs) self.assertIsNone(xco_cs) else: + self.assertIsNotNone(cube_cs) if cube_cstype is not None: self.assertIsInstance(cube_cs, cube_cstype) if xco_no_cs: @@ -755,6 +761,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 +883,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 @@ -869,5 +963,279 @@ def test_nondim_lats(self): self.check_result(result, yco_is_aux=True, load_problems_regex=error) +@pytest.fixture +def latlon_cs(): + return """ + int crsWGS84 ; + crsWGS84:grid_mapping_name = "latitude_longitude" ; + crsWGS84:longitude_of_prime_meridian = 0. ; + crsWGS84:semi_major_axis = 6378137. ; + crsWGS84:inverse_flattening = 298.257223563 ; + double lat(x,y) ; + lat:standard_name = "latitude" ; + lat:units = "degrees_north" ; + double lon(x,y) ; + lon:standard_name = "longitude" ; + lon:units = "degrees_east" ; + """ + + +@pytest.fixture +def latlon_cs_missing_coord(latlon_cs): + # Trin last 3 lines of latlon_cs to remove longitude coord + return "\n".join(latlon_cs.splitlines()[:-4]) + + +@pytest.fixture +def osgb_cs(): + return """ + int crsOSGB ; + crsOSGB:grid_mapping_name = "transverse_mercator" ; + crsOSGB:semi_major_axis = 6377563.396 ; + crsOSGB:inverse_flattening = 299.3249646 ; + crsOSGB:longitude_of_prime_meridian = 0. ; + crsOSGB:latitude_of_projection_origin = 49. ; + crsOSGB:longitude_of_central_meridian = -2. ; + crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ; + crsOSGB:false_easting = 400000. ; + crsOSGB:false_northing = -100000. ; + crsOSGB:unit = "metre" ; + double x(x) ; + x:standard_name = "projection_x_coordinate" ; + x:long_name = "Easting" ; + x:units = "m" ; + double y(y) ; + y:standard_name = "projection_y_coordinate" ; + y:long_name = "Northing" ; + y:units = "m" ; + """ + + +class Test_multi_coordinate_system_grid_mapping(Mixin__nc_load_actions): + def test_two_coord_systems(self, osgb_cs, latlon_cs, mocker, tmp_path): + """Test load a well described multi coordinate system variable.""" + cdl = f""" + netcdf tmp {{ + dimensions: + x = 4 ; + y = 3 ; + variables: + float phenom(y, x) ; + phenom:standard_name = "air_pressure" ; + phenom:units = "Pa" ; + phenom:coordinates = "lat lon" ; + phenom:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ; + {osgb_cs} + {latlon_cs} + }} + """ + nc_path = str(tmp_path / "tmp.nc") + cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker) + + assert len(cube.coord_systems()) == 2 + + for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]: + assert isinstance( + cube.coord(coord_name).coord_system, + iris.coord_systems.TransverseMercator, + ) + for coord_name in ["longitude", "latitude"]: + assert isinstance( + cube.coord(coord_name).coord_system, iris.coord_systems.GeogCS + ) + + # Loading multiple coord systems or using extended grid mapping implies ordered axes: + assert cube.extended_grid_mapping is True + + def test_two_coord_systems_missing_coord( + self, osgb_cs, latlon_cs_missing_coord, mocker, tmp_path + ): + """Test missing coord in grid_mapping raises warning.""" + cdl = f""" + netcdf tmp {{ + dimensions: + x = 4 ; + y = 3 ; + variables: + float phenom(y, x) ; + phenom:standard_name = "air_pressure" ; + phenom:units = "Pa" ; + phenom:coordinates = "lat" ; + phenom:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ; + {osgb_cs} + {latlon_cs_missing_coord} + }} + """ + nc_path = str(tmp_path / "tmp.nc") + + # loader will warn that it can't find the longitude coord + with pytest.warns( + iris.warnings.IrisCfWarning, + match="Missing CF-netCDF coordinate variable 'lon'", + ): + cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker) + + assert len(cube.coords()) == 3 + assert len(cube.coord_systems()) == 2 + + for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]: + assert isinstance( + cube.coord(coord_name).coord_system, + iris.coord_systems.TransverseMercator, + ) + for coord_name in ["latitude"]: + assert isinstance( + cube.coord(coord_name).coord_system, iris.coord_systems.GeogCS + ) + + assert cube.extended_grid_mapping is True + + def test_two_coord_systems_missing_aux_crs( + self, osgb_cs, latlon_cs, mocker, tmp_path + ): + """Test invalid coordinate system mapping for Aux coords.""" + cdl = f""" + netcdf tmp {{ + dimensions: + x = 4 ; + y = 3 ; + variables: + float phenom(y, x) ; + phenom:standard_name = "air_pressure" ; + phenom:units = "Pa" ; + phenom:coordinates = "lat lon" ; + phenom:grid_mapping = "crsOSGB: x y non_existent_crs: lat lon" ; + {osgb_cs} + {latlon_cs} + }} + """ + nc_path = str(tmp_path / "tmp.nc") + + # loader will warn that it can't find the longitude coord + with pytest.warns( + iris.warnings.IrisCfWarning, + match="Missing CF-netCDF grid mapping variable 'non_existent_crs'", + ): + cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker) + + assert len(cube.coords()) == 4 + assert len(cube.coord_systems()) == 1 + + for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]: + assert isinstance( + cube.coord(coord_name).coord_system, + iris.coord_systems.TransverseMercator, + ) + for coord_name in ["latitude", "longitude"]: + assert cube.coord(coord_name).coord_system is None + + assert cube.extended_grid_mapping is True + + def test_two_coord_systems_missing_dim_crs( + self, osgb_cs, latlon_cs, mocker, tmp_path + ): + """Test invalid coordinate system mapping for Dim coords.""" + cdl = f""" + netcdf tmp {{ + dimensions: + x = 4 ; + y = 3 ; + variables: + float phenom(y, x) ; + phenom:standard_name = "air_pressure" ; + phenom:units = "Pa" ; + phenom:coordinates = "lat lon" ; + phenom:grid_mapping = "non_existent_crs: x y crsWGS84: lat lon" ; + {osgb_cs} + {latlon_cs} + }} + """ + nc_path = str(tmp_path / "tmp.nc") + + # loader will warn that it can't find the longitude coord + with pytest.warns( + iris.warnings.IrisCfWarning, + match="Missing CF-netCDF grid mapping variable 'non_existent_crs'", + ): + cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker) + + assert len(cube.coords()) == 2 # DimCoords won't be found + assert len(cube.coord_systems()) == 1 + + for coord_name in ["latitude", "longitude"]: + assert isinstance( + cube.coord(coord_name).coord_system, iris.coord_systems.GeogCS + ) + + assert cube.extended_grid_mapping is True + + def test_two_coord_systems_invalid_grid_mapping( + self, osgb_cs, latlon_cs, mocker, tmp_path + ): + """Test invalid grid mapping doesn't load any coord systems and warns.""" + cdl = f""" + netcdf tmp {{ + dimensions: + x = 4 ; + y = 3 ; + variables: + float phenom(y, x) ; + phenom:standard_name = "air_pressure" ; + phenom:units = "Pa" ; + phenom:coordinates = "lat lon" ; + phenom:grid_mapping = "crsOSGB: crsWGS84:" ; + {osgb_cs} + {latlon_cs} + }} + """ + nc_path = str(tmp_path / "tmp.nc") + + # loader will warn that grid_mapping is invalid + with pytest.warns( + iris.warnings.IrisCfWarning, match="Error parsing `grid_mapping` attribute" + ): + cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker) + + assert len(cube.coords()) == 2 + assert len(cube.coord_systems()) == 0 + for coord in cube.coords(): + assert coord.coord_system is None + + assert cube.extended_grid_mapping is False + + def test_one_coord_system_simple(self, osgb_cs, latlon_cs, mocker, tmp_path): + """Make sure the simple coord system syntax still works.""" + cdl = f""" + netcdf tmp {{ + dimensions: + x = 4 ; + y = 3 ; + variables: + float phenom(y, x) ; + phenom:standard_name = "air_pressure" ; + phenom:units = "Pa" ; + phenom:coordinates = "lat lon" ; + phenom:grid_mapping = "crsOSGB" ; + {osgb_cs} + {latlon_cs} + }} + """ + nc_path = str(tmp_path / "tmp.nc") + cube = self.load_cube_from_cdl(cdl, None, nc_path, mocker) + + assert len(cube.coord_systems()) == 1 + + for coord_name in ["projection_x_coordinate", "projection_x_coordinate"]: + assert isinstance( + cube.coord(coord_name).coord_system, + iris.coord_systems.TransverseMercator, + ) + for coord_name in ["longitude", "latitude"]: + assert cube.coord(coord_name).coord_system is None + + # Loading multiple coord systems or using extended grid mapping implies ordered axes: + assert cube.extended_grid_mapping is False + + if __name__ == "__main__": tests.main() diff --git a/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py b/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py index 292a21b18b..0905c3d2a9 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py +++ b/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver.py @@ -16,6 +16,7 @@ import numpy as np from numpy import ma +import pytest import iris from iris.coord_systems import ( @@ -35,6 +36,7 @@ from iris.coords import AncillaryVariable, AuxCoord, DimCoord from iris.cube import Cube from iris.fileformats.netcdf import Saver, _thread_safe_nc +from iris.tests._shared_utils import assert_CDL import iris.tests.stock as stock @@ -906,12 +908,122 @@ def test_passthrough_units(self): ) -class Test__create_cf_grid_mapping(tests.IrisTest): +@pytest.fixture +def transverse_mercator_cube_multi_cs(): + """A transverse mercator cube with an auxiliary GeogGS coordinate system.""" + data = np.arange(12).reshape(3, 4) + cube = Cube(data, "air_pressure_anomaly") + cube.extended_grid_mapping = True + + geog_cs = GeogCS(6377563.396, 6356256.909) + trans_merc = TransverseMercator( + 49.0, -2.0, -400000.0, 100000.0, 0.9996012717, geog_cs + ) + coord = DimCoord( + np.arange(3), + "projection_y_coordinate", + units="m", + coord_system=trans_merc, + ) + cube.add_dim_coord(coord, 0) + coord = DimCoord( + np.arange(4), + "projection_x_coordinate", + units="m", + coord_system=trans_merc, + ) + cube.add_dim_coord(coord, 1) + + # Add auxiliary lat/lon coords with a GeogCS coord system + coord = AuxCoord( + np.arange(3 * 4).reshape((3, 4)), + "longitude", + units="degrees", + coord_system=geog_cs, + ) + cube.add_aux_coord(coord, (0, 1)) + + coord = AuxCoord( + np.arange(3 * 4).reshape((3, 4)), + "latitude", + units="degrees", + coord_system=geog_cs, + ) + cube.add_aux_coord(coord, (0, 1)) + + return cube + + +class Test_write_extended_grid_mapping: + def test_multi_cs(self, transverse_mercator_cube_multi_cs, tmp_path, request): + """Test writing a cube with multiple coordinate systems. + Should generate a grid mapping using extended syntax that references + both coordinate systems and the coords. + """ + cube = transverse_mercator_cube_multi_cs + nc_path = tmp_path / "tmp.nc" + with Saver(nc_path, "NETCDF4") as saver: + saver.write(cube) + assert_CDL(request, nc_path) + + def test_no_aux_cs(self, transverse_mercator_cube_multi_cs, tmp_path, request): + """Test when DimCoords have coord system, but AuxCoords do not. + Should write extended grid mapping for just DimCoords. + """ + cube = transverse_mercator_cube_multi_cs + cube.coord("latitude").coord_system = None + cube.coord("longitude").coord_system = None + + nc_path = tmp_path / "tmp.nc" + with Saver(nc_path, "NETCDF4") as saver: + saver.write(cube) + assert_CDL(request, nc_path) + + def test_multi_cs_missing_coord( + self, transverse_mercator_cube_multi_cs, tmp_path, request + ): + """Test when we have a missing coordinate. + Grid mapping will fall back to simple mapping to DimCoord CS (no coords referenced). + """ + cube = transverse_mercator_cube_multi_cs + cube.remove_coord("latitude") + nc_path = tmp_path / "tmp.nc" + with Saver(nc_path, "NETCDF4") as saver: + saver.write(cube) + assert_CDL(request, nc_path) + + def test_no_cs(self, transverse_mercator_cube_multi_cs, tmp_path, request): + """Test when no coordinate systems associated with cube coords. + Grid mapping will not be generated at all. + """ + cube = transverse_mercator_cube_multi_cs + for coord in cube.coords(): + coord.coord_system = None + + nc_path = tmp_path / "tmp.nc" + with Saver(nc_path, "NETCDF4") as saver: + saver.write(cube) + assert_CDL(request, nc_path) + + +class Test_create_cf_grid_mapping: + """Tests correct generation of CF grid_mapping variable attributes. + + Note: The first 3 tests are run with the "extended grid" mapping + both enabled (the default for all these tests) and disabled. This + controls the output of the WKT attribute. + """ + + @pytest.fixture(autouse=True) + def _setup(self): + self._extended_grid_mapping = True # forces WKT strings to be written + def _cube_with_cs(self, coord_system): """Return a simple 2D cube that uses the given coordinate system.""" cube = stock.lat_lon_cube() x, y = cube.coord("longitude"), cube.coord("latitude") x.coord_system = y.coord_system = coord_system + cube.extended_grid_mapping = self._extended_grid_mapping return cube def _grid_mapping_variable(self, coord_system): @@ -930,14 +1042,15 @@ def setncattr(self, name, attr): grid_variable = NCMock(name="NetCDFVariable") create_var_fn = mock.Mock(side_effect=[grid_variable]) dataset = mock.Mock(variables=[], createVariable=create_var_fn) - saver = mock.Mock(spec=Saver, _coord_systems=[], _dataset=dataset) variable = NCMock() - # This is the method we're actually testing! - Saver._create_cf_grid_mapping(saver, cube, variable) + saver = Saver(dataset, "NETCDF4", compute=False) + + # The method we want to test: + saver._create_cf_grid_mapping(cube, variable) - self.assertEqual(create_var_fn.call_count, 1) - self.assertEqual(variable.grid_mapping, grid_variable.grid_mapping_name) + assert create_var_fn.call_count == 1 + assert variable.grid_mapping, grid_variable.grid_mapping_name return grid_variable def _variable_attributes(self, coord_system): @@ -957,12 +1070,14 @@ def _test(self, coord_system, expected): actual = self._variable_attributes(coord_system) # To see obvious differences, check that they keys are the same. - self.assertEqual(sorted(actual.keys()), sorted(expected.keys())) + assert sorted(actual.keys()) == sorted(expected.keys()) # Now check that the values are equivalent. - self.assertEqual(actual, expected) + assert actual == expected - def test_rotated_geog_cs(self): + def test_rotated_geog_cs(self, extended_grid_mapping): + self._extended_grid_mapping = extended_grid_mapping coord_system = RotatedGeogCS(37.5, 177.5, ellipsoid=GeogCS(6371229.0)) + expected = { "grid_mapping_name": b"rotated_latitude_longitude", "north_pole_grid_longitude": 0.0, @@ -971,18 +1086,47 @@ def test_rotated_geog_cs(self): "longitude_of_prime_meridian": 0.0, "earth_radius": 6371229.0, } + if extended_grid_mapping: + expected["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]]]]' + ) + self._test(coord_system, expected) - def test_spherical_geog_cs(self): + def test_spherical_geog_cs(self, extended_grid_mapping): + self._extended_grid_mapping = extended_grid_mapping coord_system = GeogCS(6371229.0) expected = { "grid_mapping_name": b"latitude_longitude", "longitude_of_prime_meridian": 0.0, "earth_radius": 6371229.0, } + if extended_grid_mapping: + expected["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]]]]' + ) self._test(coord_system, expected) - def test_elliptic_geog_cs(self): + def test_elliptic_geog_cs(self, extended_grid_mapping): + self._extended_grid_mapping = extended_grid_mapping coord_system = GeogCS(637, 600) expected = { "grid_mapping_name": b"latitude_longitude", @@ -990,6 +1134,17 @@ def test_elliptic_geog_cs(self): "semi_minor_axis": 600.0, "semi_major_axis": 637.0, } + if extended_grid_mapping: + expected["crs_wkt"] = ( + 'GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",637,' + '17.2162162162162,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],' + 'PRIMEM["Reference meridian",0,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]]]]" + ) self._test(coord_system, expected) def test_lambert_conformal(self): @@ -1002,6 +1157,25 @@ def test_lambert_conformal(self): ellipsoid=GeogCS(6371000), ) expected = { + "crs_wkt": ( + 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",' + 'ELLIPSOID["unknown",6371000,0,LENGTHUNIT["metre",1,ID["EPSG"' + ',9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",' + '0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",' + 'METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],' + 'PARAMETER["Latitude of false origin",44,ANGLEUNIT["degree",' + '0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of ' + 'false origin",2,ANGLEUNIT["degree",0.0174532925199433],ID[' + '"EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",' + '38,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],' + 'PARAMETER["Latitude of 2nd standard parallel",50,ANGLEUNIT' + '["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER[' + '"Easting at false origin",-2,LENGTHUNIT["metre",1],ID["EPSG",' + '8826]],PARAMETER["Northing at false origin",-5,LENGTHUNIT[' + '"metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,' + 'ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",' + 'north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]' + ), "grid_mapping_name": b"lambert_conformal_conic", "latitude_of_projection_origin": 44, "longitude_of_central_meridian": 2, @@ -1022,6 +1196,21 @@ def test_laea_cs(self): ellipsoid=GeogCS(6377563.396, 6356256.909), ) expected = { + "crs_wkt": ( + 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP' + 'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre' + '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree' + '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",' + 'METHOD["Lambert Azimuthal Equal Area",ID["EPSG",9820]],PARAME' + 'TER["Latitude of natural origin",52,ANGLEUNIT["degree",0.0174' + '532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natura' + 'l origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG"' + ',8802]],PARAMETER["False easting",100,LENGTHUNIT["metre",1],I' + 'D["EPSG",8806]],PARAMETER["False northing",200,LENGTHUNIT["me' + 'tre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORD' + 'ER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north' + ',ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]' + ), "grid_mapping_name": b"lambert_azimuthal_equal_area", "latitude_of_projection_origin": 52, "longitude_of_projection_origin": 10, @@ -1043,6 +1232,25 @@ def test_aea_cs(self): ellipsoid=GeogCS(6377563.396, 6356256.909), ) expected = { + "crs_wkt": ( + 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP' + 'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre' + '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree' + '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",' + 'METHOD["Albers Equal Area",ID["EPSG",9822]],PARAMETER["Latitu' + 'de of false origin",52,ANGLEUNIT["degree",0.0174532925199433]' + ',ID["EPSG",8821]],PARAMETER["Longitude of false origin",10,AN' + 'GLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMET' + 'ER["Latitude of 1st standard parallel",38,ANGLEUNIT["degree",' + '0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2' + 'nd standard parallel",50,ANGLEUNIT["degree",0.017453292519943' + '3],ID["EPSG",8824]],PARAMETER["Easting at false origin",100,L' + 'ENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at ' + 'false origin",200,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[' + 'Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID' + '["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",' + '1,ID["EPSG",9001]]]]' + ), "grid_mapping_name": b"albers_conical_equal_area", "latitude_of_projection_origin": 52, "longitude_of_central_meridian": 10, @@ -1075,6 +1283,24 @@ def test_vp_cs(self): ellipsoid=ellipsoid, ) expected = { + "crs_wkt": ( + 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP' + 'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre' + '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree' + '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",' + 'METHOD["Vertical Perspective",ID["EPSG",9838]],PARAMETER["Lat' + 'itude of topocentric origin",1,ANGLEUNIT["degree",0.017453292' + '5199433],ID["EPSG",8834]],PARAMETER["Longitude of topocentric' + ' origin",2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8' + '835]],PARAMETER["Ellipsoidal height of topocentric origin",0,' + 'LENGTHUNIT["metre",1],ID["EPSG",8836]],PARAMETER["Viewpoint h' + 'eight",2000000,LENGTHUNIT["metre",1],ID["EPSG",8840]],PARAMET' + 'ER["False easting",100,LENGTHUNIT["metre",1],ID["EPSG",8806]]' + ',PARAMETER["False northing",200,LENGTHUNIT["metre",1],ID["EPS' + 'G",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNI' + 'T["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGT' + 'HUNIT["metre",1,ID["EPSG",9001]]]]' + ), "grid_mapping_name": b"vertical_perspective", "latitude_of_projection_origin": latitude_of_projection_origin, "longitude_of_projection_origin": longitude_of_projection_origin, @@ -1109,6 +1335,23 @@ def test_geo_cs(self): ellipsoid=ellipsoid, ) expected = { + "crs_wkt": ( + 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP' + 'SOID["unknown",6377563.396,299.324961266495,LENGTHUNIT["metre' + '",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree' + '",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",' + 'METHOD["Geostationary Satellite (Sweep X)"],PARAMETER["Longit' + 'ude of natural origin",2,ANGLEUNIT["degree",0.017453292519943' + '3],ID["EPSG",8802]],PARAMETER["Satellite Height",2000000,LENG' + 'THUNIT["metre",1,ID["EPSG",9001]]],PARAMETER["False easting",' + '100,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False n' + 'orthing",200,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Carte' + 'sian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPS' + 'G",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID[' + '"EPSG",9001]]],REMARK["PROJ CRS string: +proj=geos +a=6377563' + ".396 +b=6356256.909 +lon_0=2.0 +lat_0=0.0 +h=2000000.0 +x_0=1" + '00.0 +y_0=200.0 +units=m +sweep=x +no_defs"]]' + ), "grid_mapping_name": b"geostationary", "latitude_of_projection_origin": latitude_of_projection_origin, "longitude_of_projection_origin": longitude_of_projection_origin, @@ -1126,8 +1369,30 @@ def test_oblique_cs(self): # Some none-default settings to confirm all parameters are being # handled. + wkt_template = ( + 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIP' + 'SOID["unknown",1,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PR' + 'IMEM["Reference meridian",0,ANGLEUNIT["degree",0.017453292519' + '9433,ID["EPSG",9122]]]],CONVERSION["unknown",METHOD["Hotine O' + 'blique Mercator (variant B)",ID["EPSG",9815]],PARAMETER["Lati' + 'tude of projection centre",89.9,ANGLEUNIT["degree",0.01745329' + '25199433],ID["EPSG",8811]],PARAMETER["Longitude of projection' + ' centre",45,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",' + '8812]],PARAMETER["Azimuth at projection centre",{angle},ANGLEUNIT[' + '"degree",0.0174532925199433],ID["EPSG",8813]],PARAMETER["Angl' + 'e from Rectified to Skew Grid",{angle},ANGLEUNIT["degree",0.017453' + '2925199433],ID["EPSG",8814]],PARAMETER["Scale factor at proje' + 'ction centre",0.939692620786,SCALEUNIT["unity",1],ID["EPSG",8' + '815]],PARAMETER["Easting at projection centre",1000000,LENGTH' + 'UNIT["metre",1],ID["EPSG",8816]],PARAMETER["Northing at proje' + 'ction centre",-2000000,LENGTHUNIT["metre",1],ID["EPSG",8817]]' + '],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre"' + ',1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["me' + 'tre",1,ID["EPSG",9001]]]]' + ) + kwargs_rotated = dict( - latitude_of_projection_origin=90.0, + latitude_of_projection_origin=89.9, longitude_of_projection_origin=45.0, false_easting=1000000.0, false_northing=-2000000.0, @@ -1142,8 +1407,9 @@ def test_oblique_cs(self): expected_rotated = dict( # Automatically converted to oblique_mercator in line with CF 1.11 . grid_mapping_name=b"oblique_mercator", - # Azimuth should be automatically populated. + # Azimuth and crs_wkt should be automatically populated. azimuth_of_central_line=90.0, + crs_wkt=wkt_template.format(angle="89.999"), **kwargs_rotated, ) # Convert the ellipsoid @@ -1156,6 +1422,7 @@ def test_oblique_cs(self): # Same as rotated, but different azimuth. expected_oblique = dict(expected_rotated, **oblique_azimuth) + expected_oblique["crs_wkt"] = wkt_template.format(angle="45") oblique = ObliqueMercator(**kwargs_oblique) rotated = RotatedMercator(**kwargs_rotated) @@ -1167,5 +1434,16 @@ def test_oblique_cs(self): self._test(coord_system, expected) +@pytest.fixture( + params=[ + pytest.param(True, id="extended_grid_mapping"), + pytest.param(False, id="no_extended_grid_mapping"), + ] +) +def extended_grid_mapping(request): + """Fixture for enabling/disabling extended grid mapping.""" + return request.param + + if __name__ == "__main__": tests.main()