diff --git a/lib/iris/fileformats/_nc_load_rules/__init__.py b/lib/iris/fileformats/_nc_load_rules/__init__.py new file mode 100644 index 0000000000..baea3cf555 --- /dev/null +++ b/lib/iris/fileformats/_nc_load_rules/__init__.py @@ -0,0 +1,16 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Support for cube-specific CF-to-Iris translation operations. + +Interprets CF concepts identified by :mod:`iris.fileformats.cf` to add +components into loaded cubes. + +For now : the API which mimics :class:`pyke.knowledge_engine.engine`. +As this is aiming to replace the old Pyke-based logic rules. +TODO: simplify once the parallel operation with Pyke is no longer required. + +""" diff --git a/lib/iris/fileformats/_nc_load_rules/actions.py b/lib/iris/fileformats/_nc_load_rules/actions.py new file mode 100644 index 0000000000..4676389cb7 --- /dev/null +++ b/lib/iris/fileformats/_nc_load_rules/actions.py @@ -0,0 +1,562 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Replacement code for the Pyke rules. + +For now, we are still emulating various aspects of how our original Pyke-based +code used the Pyke 'engine' to hold translation data, both Pyke-specific and +not : +1) basic details from the iris.fileformats.cf analysis of the file are + recorded before translating each output cube, using + "engine.assert_case_specific_fact(name, args)". + +2) this is also used to store intermediate info passed between rules, which + used to be done with a "facts_cf.provides" statement in rule actions. + +3) Iris-specific info is (still) stored in additional properties created on + the engine object : + engine.cf_var, .cube, .cube_parts, .requires, .rule_triggered, .filename + +Our "rules" are just action routines. +The top-level 'run_actions' routine decides which actions to call, based on the +info recorded when processing each cube output. It does this in a simple +explicit way, which doesn't use any clever chaining, "trigger conditions" or +other rule-type logic. + +Each 'action' function can replace several similar 'rules'. +E.G. 'action_provides_grid_mapping' replaces all 'fc_provides_grid_mapping_'. +To aid debug, each returns a 'rule_name' string, indicating which original rule +this particular action call is emulating : In some cases, this may include a +textual note that this rule 'failed', aka "did not trigger", which would not be +recorded in the original implementation. + +The top-level 'run_actions' ensures that the individual rules actions are +called, with various arguments, as appropriate to ensure the whole cube is +built as it was by the original rules implementation. + +TODO: remove the use of intermediate "facts" to carry information between +actions. This mimics older behaviour, so is still useful while we are still +comparing behaviour with the old Pyke rules (debugging). But once that is no +longer useful, this can be considerably simplified. + +""" + +from . import helpers as hh + +from functools import wraps + +import iris.fileformats.pp as pp +import iris.fileformats.cf + + +def _default_rulenamesfunc(func_name): + # A simple default function to deduce the rules-name from an action-name. + funcname_prefix = "action_" + rulename_prefix = "fc_" # To match existing behaviours + rule_name = func_name + if rule_name.startswith(funcname_prefix): + rule_name = rule_name[len(funcname_prefix) :] + if not rule_name.startswith(rulename_prefix): + rule_name = rulename_prefix + rule_name + return rule_name + + +def action_function(func): + # Wrap an action function with some standard behaviour. + # Notably : engages with the rules logging process. + @wraps(func) + def inner(engine, *args, **kwargs): + # Call the original rules-func + rule_name = func(engine, *args, **kwargs) + if rule_name is None: + # Work out the corresponding rule name, and log it. + # Note: an action returns a name string, which identifies it, + # but also may vary depending on whether it successfully + # triggered, and if so what it matched. + rule_name = _default_rulenamesfunc(func.__name__) + engine.rule_triggered.add(rule_name) + + func._rulenames_func = _default_rulenamesfunc + return inner + + +@action_function +def action_default(engine): + """Standard operations for every cube.""" + hh.build_cube_metadata(engine) + + +# Lookup table used by 'action_provides_grid_mapping'. +# Maps each supported CF grid-mapping-name to a pair of handling ("helper") +# routines: +# (@0) a validity-checker (or None) +# (@1) a coord-system builder function. +_GRIDTYPE_CHECKER_AND_BUILDER = { + hh.CF_GRID_MAPPING_LAT_LON: (None, hh.build_coordinate_system), + hh.CF_GRID_MAPPING_ROTATED_LAT_LON: ( + None, + hh.build_rotated_coordinate_system, + ), + hh.CF_GRID_MAPPING_MERCATOR: ( + hh.has_supported_mercator_parameters, + hh.build_mercator_coordinate_system, + ), + hh.CF_GRID_MAPPING_TRANSVERSE: ( + None, + hh.build_transverse_mercator_coordinate_system, + ), + hh.CF_GRID_MAPPING_STEREO: ( + hh.has_supported_stereographic_parameters, + hh.build_stereographic_coordinate_system, + ), + hh.CF_GRID_MAPPING_LAMBERT_CONFORMAL: ( + None, + hh.build_lambert_conformal_coordinate_system, + ), + hh.CF_GRID_MAPPING_LAMBERT_AZIMUTHAL: ( + None, + hh.build_lambert_azimuthal_equal_area_coordinate_system, + ), + hh.CF_GRID_MAPPING_ALBERS: ( + None, + hh.build_albers_equal_area_coordinate_system, + ), + hh.CF_GRID_MAPPING_VERTICAL: ( + None, + hh.build_vertical_perspective_coordinate_system, + ), + hh.CF_GRID_MAPPING_GEOSTATIONARY: ( + None, + hh.build_geostationary_coordinate_system, + ), +} + + +@action_function +def action_provides_grid_mapping(engine, gridmapping_fact): + """Convert a CFGridMappingVariable into a cube coord-system.""" + (var_name,) = gridmapping_fact + rule_name = "fc_provides_grid_mapping" + cf_var = engine.cf_var.cf_group[var_name] + grid_mapping_type = getattr(cf_var, hh.CF_ATTR_GRID_MAPPING_NAME, None) + + succeed = True + if grid_mapping_type is None: + succeed = False + rule_name += " --FAILED(no grid-mapping attr)" + else: + grid_mapping_type = grid_mapping_type.lower() + + if succeed: + if grid_mapping_type in _GRIDTYPE_CHECKER_AND_BUILDER: + checker, builder = _GRIDTYPE_CHECKER_AND_BUILDER[grid_mapping_type] + rule_name += f"_({grid_mapping_type})" + else: + succeed = False + rule_name += f" --FAILED(unhandled type {grid_mapping_type})" + + if succeed: + if checker is not None and not checker(engine, var_name): + succeed = False + rule_name += f" --(FAILED check {checker.__name__})" + + if succeed: + coordinate_system = builder(engine, cf_var) + # Check there is not an existing one. + old_gridtype_fact = engine.fact_list("grid-type") + if old_gridtype_fact: + (old_gridtype,) = old_gridtype_fact + succeed = False + rule_name += ( + f" --(FAILED overwrite coord-sytem " + f"{old_gridtype} with {grid_mapping_type})" + ) + + if succeed: + engine.cube_parts["coordinate_system"] = coordinate_system + engine.add_fact("grid-type", (grid_mapping_type,)) + + return rule_name + + +@action_function +def action_provides_coordinate(engine, dimcoord_fact): + """Identify the coordinate 'type' of a CFCoordinateVariable.""" + (var_name,) = dimcoord_fact + + # Identify the "type" of a coordinate variable + coord_type = None + # NOTE: must test for rotated cases *first*, as 'is_longitude' and + # 'is_latitude' functions also accept rotated cases. + if hh.is_rotated_latitude(engine, var_name): + coord_type = "rotated_latitude" + elif hh.is_rotated_longitude(engine, var_name): + coord_type = "rotated_longitude" + elif hh.is_latitude(engine, var_name): + coord_type = "latitude" + elif hh.is_longitude(engine, var_name): + coord_type = "longitude" + elif hh.is_time(engine, var_name): + coord_type = "time" + elif hh.is_time_period(engine, var_name): + coord_type = "time_period" + elif hh.is_projection_x_coordinate(engine, var_name): + coord_type = "projection_x" + elif hh.is_projection_y_coordinate(engine, var_name): + coord_type = "projection_y" + + if coord_type is None: + # Not identified as a specific known coord_type. + # N.B. in the original rules, this does *not* trigger separate + # 'provides' and 'build' phases : there is just a single + # 'fc_default_coordinate' rule. + # Rationalise this for now by making it more like the others. + # FOR NOW: ~matching old code, but they could *all* be simplified. + # TODO: combine 2 operation into 1 for ALL of these. + coord_type = "miscellaneous" + rule_name = "fc_default_coordinate_(provide-phase)" + else: + rule_name = f"fc_provides_coordinate_({coord_type})" + + engine.add_fact("provides-coordinate-(oftype)", (coord_type, var_name)) + return rule_name + + +# Lookup table used by 'action_build_dimension_coordinate'. +# Maps each supported coordinate-type name (a rules-internal concept) to a pair +# of information values : +# (@0) A grid "type", one of latlon/rotated/projected (or None) +# If set, the cube should have a coord-system, which is set on the +# resulting coordinate. If None, the coord has no coord_system. +# (@1) an (optional) fixed standard-name for the coordinate, or None +# If None, the coordinate name is copied from the source variable +_COORDTYPE_GRIDTYPES_AND_COORDNAMES = { + "latitude": ("latlon", hh.CF_VALUE_STD_NAME_LAT), + "longitude": ("latlon", hh.CF_VALUE_STD_NAME_LON), + "rotated_latitude": ( + "rotated", + hh.CF_VALUE_STD_NAME_GRID_LAT, + ), + "rotated_longitude": ( + "rotated", + hh.CF_VALUE_STD_NAME_GRID_LON, + ), + "projection_x": ("projected", hh.CF_VALUE_STD_NAME_PROJ_X), + "projection_y": ("projected", hh.CF_VALUE_STD_NAME_PROJ_Y), + "time": (None, None), + "time_period": (None, None), + "miscellaneous": (None, None), +} + + +@action_function +def action_build_dimension_coordinate(engine, providescoord_fact): + """Convert a CFCoordinateVariable into a cube dim-coord.""" + coord_type, var_name = providescoord_fact + cf_var = engine.cf_var.cf_group[var_name] + rule_name = f"fc_build_coordinate_({coord_type})" + coord_grid_class, coord_name = _COORDTYPE_GRIDTYPES_AND_COORDNAMES[ + coord_type + ] + if coord_grid_class is None: + # Coordinates not identified with a specific grid-type class (latlon, + # rotated or projected) are always built, but can have no coord-system. + coord_system = None # no coord-system can be used + succeed = True + else: + grid_classes = ("latlon", "rotated", "projected") + assert coord_grid_class in grid_classes + # If a coord is of a type identified with a grid, we may have a + # coordinate system (i.e. a valid grid-mapping). + # N.B. this requires each grid-type identification to validate the + # coord var (e.g. "is_longitude"). + # 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") + # Translate the specific grid-mapping type to a grid-class + if coord_system is None: + succeed = True + cs_gridclass = None + else: + # 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 + if cs_gridtype == "latitude_longitude": + cs_gridclass = "latlon" + elif cs_gridtype == "rotated_latitude_longitude": + cs_gridclass = "rotated" + else: + # Other specific projections + assert cs_gridtype is not None + cs_gridclass = "projected" + + assert cs_gridclass in grid_classes + (None,) + + if coord_grid_class == "latlon": + if cs_gridclass == "latlon": + succeed = True + elif cs_gridclass is None: + succeed = True + rule_name += "(no-cs)" + elif cs_gridclass == "rotated": + # We disallow this case + succeed = False + rule_name += "(FAILED : latlon coord with rotated cs)" + else: + assert cs_gridclass == "projected" + # succeed, no error, but discards the coord-system + # TODO: could issue a warning in this case ? + succeed = True + coord_system = None + rule_name += "(no-cs : discarded projected cs)" + elif coord_grid_class == "rotated": + if cs_gridclass == "rotated": + succeed = True + rule_name += "(rotated)" + elif cs_gridclass is None: + succeed = True + rule_name += "(rotated no-cs)" + elif cs_gridclass == "latlon": + # We disallow this case + succeed = False + rule_name += "(FAILED rotated coord with latlon cs)" + else: + assert cs_gridclass == "projected" + succeed = True + coord_system = None + rule_name += "(rotated no-cs : discarded projected cs)" + elif coord_grid_class == "projected": + # In this case, can *only* build a coord at all if there is a + # coord-system of the correct class (i.e. 'projected'). + succeed = cs_gridclass == "projected" + if not succeed: + rule_name += "(FAILED projected coord with non-projected cs)" + else: + msg = ( + f'Unexpected coord grid-class "{coord_grid_class}" ' + f"for coord {var_name}." + ) + raise ValueError(msg) + + if succeed: + hh.build_dimension_coordinate( + engine, cf_var, coord_name=coord_name, coord_system=coord_system + ) + return rule_name + + +@action_function +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" + + # Identify any known coord "type" : latitude/longitude/time/time_period + # If latitude/longitude, this sets the standard_name of the built AuxCoord + coord_type = "" # unidentified : can be OK + coord_name = None + if hh.is_time(engine, var_name): + coord_type = "time" + elif hh.is_time_period(engine, var_name): + coord_type = "time_period" + elif hh.is_longitude(engine, var_name): + coord_type = "longitude" + if hh.is_rotated_longitude(engine, var_name): + coord_type += "_rotated" + coord_name = hh.CF_VALUE_STD_NAME_GRID_LON + else: + coord_name = hh.CF_VALUE_STD_NAME_LON + elif hh.is_latitude(engine, var_name): + coord_type = "latitude" + if hh.is_rotated_latitude(engine, var_name): + coord_type += "_rotated" + coord_name = hh.CF_VALUE_STD_NAME_GRID_LAT + else: + coord_name = hh.CF_VALUE_STD_NAME_LAT + + if coord_type: + rule_name += f"_{coord_type}" + + cf_var = engine.cf_var.cf_group.auxiliary_coordinates[var_name] + hh.build_auxiliary_coordinate(engine, cf_var, coord_name=coord_name) + + return rule_name + + +@action_function +def action_ukmo_stash(engine): + """Convert 'ukmo stash' cf property into a cube attribute.""" + rule_name = "fc_attribute_ukmo__um_stash_source" + var = engine.cf_var + attr_name = "ukmo__um_stash_source" + attr_value = getattr(var, attr_name, None) + if attr_value is None: + attr_altname = "um_stash_source" # legacy form + attr_value = getattr(var, attr_altname, None) + if attr_value is None: + rule_name += "(NOT-TRIGGERED)" + else: + # No helper routine : just do it + engine.cube.attributes["STASH"] = pp.STASH.from_msi(attr_value) + + return rule_name + + +@action_function +def action_ukmo_processflags(engine): + """Convert 'ukmo process flags' cf property into a cube attribute.""" + rule_name = "fc_attribute_ukmo__process_flags" + var = engine.cf_var + attr_name = "ukmo__process_flags" + attr_value = getattr(var, attr_name, None) + if attr_value is None: + rule_name += "(NOT-TRIGGERED)" + else: + # No helper routine : just do it + flags = [x.replace("_", " ") for x in attr_value.split(" ")] + engine.cube.attributes["ukmo__process_flags"] = tuple(flags) + + return rule_name + + +@action_function +def action_build_cell_measure(engine, cellm_fact): + """Convert a CFCellMeasureVariable into a cube cell-measure.""" + (var_name,) = cellm_fact + var = engine.cf_var.cf_group.cell_measures[var_name] + hh.build_cell_measures(engine, var) + + +@action_function +def action_build_ancil_var(engine, ancil_fact): + """Convert a CFAncillaryVariable into a cube ancil-var.""" + (var_name,) = ancil_fact + var = engine.cf_var.cf_group.ancillary_variables[var_name] + hh.build_ancil_var(engine, var) + + +@action_function +def action_build_label_coordinate(engine, label_fact): + """Convert a CFLabelVariable into a cube string-type aux-coord.""" + (var_name,) = label_fact + var = engine.cf_var.cf_group.labels[var_name] + hh.build_auxiliary_coordinate(engine, var) + + +@action_function +def action_formula_type(engine, formula_root_fact): + """Register a CFVariable as a formula root.""" + rule_name = "fc_formula_type" + (var_name,) = formula_root_fact + cf_var = engine.cf_var.cf_group[var_name] + # var.standard_name is a formula type (or we should never get here). + formula_type = getattr(cf_var, "standard_name", None) + succeed = True + if formula_type not in iris.fileformats.cf.reference_terms: + succeed = False + rule_name += f"(FAILED - unrecognised formula type = {formula_type!r})" + if succeed: + # Check we don't already have one. + existing_type = engine.requires.get("formula_type") + if existing_type: + succeed = False + rule_name += ( + f"(FAILED - new formula type ={formula_type!r} " + f"collided with existing one ={existing_type!r}.)" + ) + if succeed: + rule_name += f"_{formula_type}" + # Set 'requires' info for iris.fileformats.netcdf._load_aux_factory. + engine.requires["formula_type"] = formula_type + + return rule_name + + +@action_function +def action_formula_term(engine, formula_term_fact): + """Register a CFVariable as a formula term.""" + # Must run AFTER formula root identification. + (termvar_name, rootvar_name, term_name) = formula_term_fact + # The rootname is implicit : have only one per cube + # TODO: change when we adopt cf-1.7 advanced grid-mping syntax + engine.requires.setdefault("formula_terms", {})[term_name] = termvar_name + rule_name = f"fc_formula_term({term_name})" + return rule_name + + +def run_actions(engine): + """ + Run all actions for a cube. + + This is the top-level "activation" function which runs all the appropriate + rules actions to translate facts and build all the cube elements. + + The specific cube being translated is "engine.cube". + + """ + + # default (all cubes) action, always runs + action_default(engine) # This should run the default rules. + + # deal with grid-mappings + 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) + + # identify + record aka "PROVIDE" specific named coordinates + # N.B. cf.py has identified that these are dim-coords, NOT aux-coords + # (which are recorded separately). + # TODO: can probably remove this step ?? + dimcoord_facts = engine.fact_list("coordinate") + for dimcoord_fact in dimcoord_facts: + action_provides_coordinate(engine, dimcoord_fact) + + # build (dimension) coordinates + # The 'provides' step and the grid-mapping must have already been done. + providescoord_facts = engine.fact_list("provides-coordinate-(oftype)") + for providescoord_fact in providescoord_facts: + action_build_dimension_coordinate(engine, providescoord_fact) + + # build aux-coords + auxcoord_facts = engine.fact_list("auxiliary_coordinate") + for auxcoord_fact in auxcoord_facts: + action_build_auxiliary_coordinate(engine, auxcoord_fact) + + # Detect + process and special 'ukmo' attributes + # Run on every cube : they choose themselves whether to trigger. + action_ukmo_stash(engine) + action_ukmo_processflags(engine) + + # cell measures + cellm_facts = engine.fact_list("cell_measure") + for cellm_fact in cellm_facts: + action_build_cell_measure(engine, cellm_fact) + + # ancillary variables + ancil_facts = engine.fact_list("ancillary_variable") + for ancil_fact in ancil_facts: + action_build_ancil_var(engine, ancil_fact) + + # label coords + label_facts = engine.fact_list("label") + for label_fact in label_facts: + action_build_label_coordinate(engine, label_fact) + + # formula root variables + formula_root_facts = engine.fact_list("formula_root") + for root_fact in formula_root_facts: + action_formula_type(engine, root_fact) + + # formula terms + # The 'formula_root's must have already been done. + formula_term_facts = engine.fact_list("formula_term") + for term_fact in formula_term_facts: + action_formula_term(engine, term_fact) diff --git a/lib/iris/fileformats/_nc_load_rules/engine.py b/lib/iris/fileformats/_nc_load_rules/engine.py new file mode 100644 index 0000000000..60f956d4d1 --- /dev/null +++ b/lib/iris/fileformats/_nc_load_rules/engine.py @@ -0,0 +1,131 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +A simple mimic of the Pyke 'knowledge_engine', for interfacing to the routines +in 'iris.fileformats.netcdf' with minimal changes to that code. + +This allows us to replace the Pyke rules operation with the simpler pure-Python +translation operations in :mod:`iris.fileformats._nc_load_rules.actions`. + +The core of this is the 'Engine' class, which mimics the Pyke engine operations, +as used by our code to translate each data cube. + +engine.get_kb() also returns a FactEntity object, which mimics *just enough* +API of a Pyke.knowlege_base, so that we can list its case-specific facts, as +used in :meth:`iris.fileformats.netcdf.pyke_stats`. + +""" +from .actions import run_actions + + +class FactList: + def __init__(self): + self.case_specific_facts = [] + + +class FactEntity: + # To support: + """ + kb_facts = engine.get_kb(_PYKE_FACT_BASE) + + for key in kb_facts.entity_lists.keys(): + for arg in kb_facts.entity_lists[key].case_specific_facts: + print("\t%s%s" % (key, arg)) + + """ + + def __init__(self): + self.entity_lists = {} + + def add_fact(self, fact_name, args): + if fact_name not in self.entity_lists: + self.entity_lists[fact_name] = FactList() + fact_list = self.entity_lists[fact_name] + fact_list.case_specific_facts.append(tuple(args)) + + def sect_facts(self, entity_name): + if entity_name in self.entity_lists: + facts = self.entity_lists.get(entity_name).case_specific_facts + else: + facts = [] + return facts + + +class Engine: + """ + A minimal mimic of a Pyke.engine. + + Provides just enough API so that the existing code in + :mod:`iris.fileformats.netcdf` can interface with our new rules functions. + + """ + + def __init__(self): + """Init new engine.""" + self.reset() + + def reset(self): + """Reset the engine = remove all facts.""" + self.facts = FactEntity() + + def activate(self, rules_base_str=None): + """ + Run all the translation rules to produce a single output cube. + + This implicitly references the output variable for this operation, + set by engine.cf_var (the variable name). + + The rules operation itself is coded elsewhere, + in :mod:`iris.fileformats.netcdf._nc_load_rules.rules`. + + """ + run_actions(self) + + def print_stats(self): + """No-op, called by :meth:`iris.fileformats.netcdf.pyke_stats`.""" + pass + + def add_case_specific_fact(self, kb_name, fact_name, fact_arglist): + """ + Record a fact about the current output operation. + + Roughly, self.facts.entity_lists[fact_name].append(fact_arglist). + + """ + self.facts.add_fact(fact_name, fact_arglist) + + def get_kb(self, fact_base_str=None): + """ + Get a FactEntity, which mimic (bits of) a knowledge-base. + + Just allowing + :meth:`iris.fileformats.netcdf.pyke_stats` to list the facts. + + """ + return self.facts + + def fact_list(self, fact_name): + """ + Return the facts (arg-lists) for one fact name. + + A shorthand form used only by the new rules routines. + + AKA 'case-specific-facts', in the original. + Roughly "return self.facts.entity_lists[fact_name]". + + """ + return self.facts.sect_facts(fact_name) + + def add_fact(self, fact_name, fact_arglist): + """ + Add a new fact. + + A shorthand form used only by the new rules routines. + + """ + self.add_case_specific_fact( + kb_name="", fact_name=fact_name, fact_arglist=fact_arglist + ) diff --git a/lib/iris/fileformats/_nc_load_rules/helpers.py b/lib/iris/fileformats/_nc_load_rules/helpers.py new file mode 100644 index 0000000000..ce7a194b35 --- /dev/null +++ b/lib/iris/fileformats/_nc_load_rules/helpers.py @@ -0,0 +1,1308 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +All the pure-Python 'helper' functions which were previously included in the +Pyke rules database 'fc_rules_cf.krb'. + +Initially these haven't changed. +The new rules approach is still calling most of them. + +""" + +import warnings + +import cf_units +import numpy as np +import numpy.ma as ma + +import iris.aux_factory +from iris.common.mixin import _get_valid_standard_name +import iris.coords +import iris.coord_systems +import iris.fileformats.cf as cf +import iris.fileformats.netcdf +from iris.fileformats.netcdf import ( + _get_cf_var_data, + parse_cell_methods, + UnknownCellMethodWarning, +) +import iris.exceptions +import iris.std_names +import iris.util + + +# +# UD Units Constants (based on Unidata udunits.dat definition file) +# +UD_UNITS_LAT = [ + "degrees_north", + "degree_north", + "degree_n", + "degrees_n", + "degreen", + "degreesn", + "degrees", + "degrees north", + "degree north", + "degree n", + "degrees n", +] +UD_UNITS_LON = [ + "degrees_east", + "degree_east", + "degree_e", + "degrees_e", + "degreee", + "degreese", + "degrees", + "degrees east", + "degree east", + "degree e", + "degrees e", +] +UNKNOWN_UNIT_STRING = "?" +NO_UNIT_STRING = "-" + +# +# CF Dimensionless Vertical Coordinates +# +CF_COORD_VERTICAL = { + "atmosphere_ln_pressure_coordinate": ["p0", "lev"], + "atmosphere_sigma_coordinate": ["sigma", "ps", "ptop"], + "atmosphere_hybrid_sigma_pressure_coordinate": ["a", "b", "ps", "p0"], + "atmosphere_hybrid_height_coordinate": ["a", "b", "orog"], + "atmosphere_sleve_coordinate": [ + "a", + "b1", + "b2", + "ztop", + "zsurf1", + "zsurf2", + ], + "ocean_sigma_coordinate": ["sigma", "eta", "depth"], + "ocean_s_coordinate": ["s", "eta", "depth", "a", "b", "depth_c"], + "ocean_sigma_z_coordinate": [ + "sigma", + "eta", + "depth", + "depth_c", + "nsigma", + "zlev", + ], + "ocean_double_sigma_coordinate": [ + "sigma", + "depth", + "z1", + "z2", + "a", + "href", + "k_c", + ], + "ocean_s_coordinate_g1": ["s", "eta", "depth", "depth_c", "C"], + "ocean_s_coordinate_g2": ["s", "eta", "depth", "depth_c", "C"], +} + +# +# CF Grid Mappings +# +CF_GRID_MAPPING_ALBERS = "albers_conical_equal_area" +CF_GRID_MAPPING_AZIMUTHAL = "azimuthal_equidistant" +CF_GRID_MAPPING_LAMBERT_AZIMUTHAL = "lambert_azimuthal_equal_area" +CF_GRID_MAPPING_LAMBERT_CONFORMAL = "lambert_conformal_conic" +CF_GRID_MAPPING_LAMBERT_CYLINDRICAL = "lambert_cylindrical_equal_area" +CF_GRID_MAPPING_LAT_LON = "latitude_longitude" +CF_GRID_MAPPING_MERCATOR = "mercator" +CF_GRID_MAPPING_ORTHO = "orthographic" +CF_GRID_MAPPING_POLAR = "polar_stereographic" +CF_GRID_MAPPING_ROTATED_LAT_LON = "rotated_latitude_longitude" +CF_GRID_MAPPING_STEREO = "stereographic" +CF_GRID_MAPPING_TRANSVERSE = "transverse_mercator" +CF_GRID_MAPPING_VERTICAL = "vertical_perspective" +CF_GRID_MAPPING_GEOSTATIONARY = "geostationary" + +# +# CF Attribute Names. +# +CF_ATTR_AXIS = "axis" +CF_ATTR_BOUNDS = "bounds" +CF_ATTR_CALENDAR = "calendar" +CF_ATTR_CLIMATOLOGY = "climatology" +CF_ATTR_GRID_INVERSE_FLATTENING = "inverse_flattening" +CF_ATTR_GRID_EARTH_RADIUS = "earth_radius" +CF_ATTR_GRID_MAPPING_NAME = "grid_mapping_name" +CF_ATTR_GRID_NORTH_POLE_LAT = "grid_north_pole_latitude" +CF_ATTR_GRID_NORTH_POLE_LON = "grid_north_pole_longitude" +CF_ATTR_GRID_NORTH_POLE_GRID_LON = "north_pole_grid_longitude" +CF_ATTR_GRID_SEMI_MAJOR_AXIS = "semi_major_axis" +CF_ATTR_GRID_SEMI_MINOR_AXIS = "semi_minor_axis" +CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN = "latitude_of_projection_origin" +CF_ATTR_GRID_LON_OF_PROJ_ORIGIN = "longitude_of_projection_origin" +CF_ATTR_GRID_STANDARD_PARALLEL = "standard_parallel" +CF_ATTR_GRID_FALSE_EASTING = "false_easting" +CF_ATTR_GRID_FALSE_NORTHING = "false_northing" +CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN = "scale_factor_at_projection_origin" +CF_ATTR_GRID_SCALE_FACTOR_AT_CENT_MERIDIAN = "scale_factor_at_central_meridian" +CF_ATTR_GRID_LON_OF_CENT_MERIDIAN = "longitude_of_central_meridian" +CF_ATTR_GRID_STANDARD_PARALLEL = "standard_parallel" +CF_ATTR_GRID_PERSPECTIVE_HEIGHT = "perspective_point_height" +CF_ATTR_GRID_SWEEP_ANGLE_AXIS = "sweep_angle_axis" +CF_ATTR_POSITIVE = "positive" +CF_ATTR_STD_NAME = "standard_name" +CF_ATTR_LONG_NAME = "long_name" +CF_ATTR_UNITS = "units" +CF_ATTR_CELL_METHODS = "cell_methods" + +# +# CF Attribute Value Constants. +# +# Attribute - axis. +CF_VALUE_AXIS_X = "x" +CF_VALUE_AXIS_Y = "y" +CF_VALUE_AXIS_T = "t" +CF_VALUE_AXIS_Z = "z" + + +# Attribute - positive. +CF_VALUE_POSITIVE = ["down", "up"] + +# Attribute - standard_name. +CF_VALUE_STD_NAME_LAT = "latitude" +CF_VALUE_STD_NAME_LON = "longitude" +CF_VALUE_STD_NAME_GRID_LAT = "grid_latitude" +CF_VALUE_STD_NAME_GRID_LON = "grid_longitude" +CF_VALUE_STD_NAME_PROJ_X = "projection_x_coordinate" +CF_VALUE_STD_NAME_PROJ_Y = "projection_y_coordinate" + + +################################################################################ +def build_cube_metadata(engine): + """Add the standard meta data to the cube.""" + + cf_var = engine.cf_var + cube = engine.cube + + # Determine the cube's name attributes + cube.var_name = cf_var.cf_name + standard_name = getattr(cf_var, CF_ATTR_STD_NAME, None) + long_name = getattr(cf_var, CF_ATTR_LONG_NAME, None) + cube.long_name = long_name + + if standard_name is not None: + try: + cube.standard_name = _get_valid_standard_name(standard_name) + except ValueError: + if cube.long_name is not None: + cube.attributes["invalid_standard_name"] = standard_name + else: + cube.long_name = standard_name + + # Determine the cube units. + attr_units = get_attr_units(cf_var, cube.attributes) + cube.units = attr_units + + # Incorporate cell methods + nc_att_cell_methods = getattr(cf_var, CF_ATTR_CELL_METHODS, None) + with warnings.catch_warnings(record=True) as warning_records: + cube.cell_methods = parse_cell_methods(nc_att_cell_methods) + # Filter to get the warning we are interested in. + warning_records = [ + record + for record in warning_records + if issubclass(record.category, UnknownCellMethodWarning) + ] + if len(warning_records) > 0: + # Output an enhanced warning message. + warn_record = warning_records[0] + name = "{}".format(cf_var.cf_name) + msg = warn_record.message.args[0] + msg = msg.replace("variable", "variable {!r}".format(name)) + warnings.warn(message=msg, category=UnknownCellMethodWarning) + + # Set the cube global attributes. + for attr_name, attr_value in cf_var.cf_group.global_attributes.items(): + try: + cube.attributes[str(attr_name)] = attr_value + except ValueError as e: + msg = "Skipping global attribute {!r}: {}" + warnings.warn(msg.format(attr_name, str(e))) + + +################################################################################ +def _get_ellipsoid(cf_grid_var): + """Return the ellipsoid definition.""" + major = getattr(cf_grid_var, CF_ATTR_GRID_SEMI_MAJOR_AXIS, None) + minor = getattr(cf_grid_var, CF_ATTR_GRID_SEMI_MINOR_AXIS, None) + inverse_flattening = getattr( + cf_grid_var, CF_ATTR_GRID_INVERSE_FLATTENING, None + ) + + # Avoid over-specification exception. + if major is not None and minor is not None: + inverse_flattening = None + + # Check for a default spherical earth. + if major is None and minor is None and inverse_flattening is None: + major = getattr(cf_grid_var, CF_ATTR_GRID_EARTH_RADIUS, None) + + return major, minor, inverse_flattening + + +################################################################################ +def build_coordinate_system(engine, cf_grid_var): + """Create a coordinate system from the CF-netCDF grid mapping variable.""" + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + return iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + +################################################################################ +def build_rotated_coordinate_system(engine, cf_grid_var): + """Create a rotated coordinate system from the CF-netCDF grid mapping variable.""" + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + north_pole_latitude = getattr( + cf_grid_var, CF_ATTR_GRID_NORTH_POLE_LAT, 90.0 + ) + north_pole_longitude = getattr( + cf_grid_var, CF_ATTR_GRID_NORTH_POLE_LON, 0.0 + ) + if north_pole_latitude is None or north_pole_longitude is None: + warnings.warn("Rotated pole position is not fully specified") + + north_pole_grid_lon = getattr( + cf_grid_var, CF_ATTR_GRID_NORTH_POLE_GRID_LON, 0.0 + ) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + rcs = iris.coord_systems.RotatedGeogCS( + north_pole_latitude, + north_pole_longitude, + north_pole_grid_lon, + ellipsoid, + ) + + return rcs + + +################################################################################ +def build_transverse_mercator_coordinate_system(engine, cf_grid_var): + """ + Create a transverse Mercator coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_central_meridian = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_CENT_MERIDIAN, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + scale_factor_at_central_meridian = getattr( + cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_CENT_MERIDIAN, None + ) + + # The following accounts for the inconsistancy in the transverse + # mercator description within the CF spec. + if longitude_of_central_meridian is None: + longitude_of_central_meridian = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None + ) + if scale_factor_at_central_meridian is None: + scale_factor_at_central_meridian = getattr( + cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None + ) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.TransverseMercator( + latitude_of_projection_origin, + longitude_of_central_meridian, + false_easting, + false_northing, + scale_factor_at_central_meridian, + ellipsoid, + ) + + return cs + + +################################################################################ +def build_lambert_conformal_coordinate_system(engine, cf_grid_var): + """ + Create a Lambert conformal conic coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_central_meridian = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_CENT_MERIDIAN, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + standard_parallel = getattr( + cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None + ) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.LambertConformal( + latitude_of_projection_origin, + longitude_of_central_meridian, + false_easting, + false_northing, + standard_parallel, + ellipsoid, + ) + + return cs + + +################################################################################ +def build_stereographic_coordinate_system(engine, cf_grid_var): + """ + Create a stereographic coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + # Iris currently only supports Stereographic projections with a scale + # factor of 1.0. This is checked elsewhere. + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.Stereographic( + latitude_of_projection_origin, + longitude_of_projection_origin, + false_easting, + false_northing, + true_scale_lat=None, + ellipsoid=ellipsoid, + ) + + return cs + + +################################################################################ +def build_mercator_coordinate_system(engine, cf_grid_var): + """ + Create a Mercator coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + longitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None + ) + # Iris currently only supports Mercator projections with specific + # values for false_easting, false_northing, + # scale_factor_at_projection_origin and standard_parallel. These are + # checked elsewhere. + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.Mercator( + longitude_of_projection_origin, ellipsoid=ellipsoid + ) + + return cs + + +################################################################################ +def build_lambert_azimuthal_equal_area_coordinate_system(engine, cf_grid_var): + """ + Create a lambert azimuthal equal area coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.LambertAzimuthalEqualArea( + latitude_of_projection_origin, + longitude_of_projection_origin, + false_easting, + false_northing, + ellipsoid, + ) + + return cs + + +################################################################################ +def build_albers_equal_area_coordinate_system(engine, cf_grid_var): + """ + Create a albers conical equal area coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_central_meridian = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_CENT_MERIDIAN, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + standard_parallels = getattr( + cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None + ) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.AlbersEqualArea( + latitude_of_projection_origin, + longitude_of_central_meridian, + false_easting, + false_northing, + standard_parallels, + ellipsoid, + ) + + return cs + + +################################################################################ +def build_vertical_perspective_coordinate_system(engine, cf_grid_var): + """ + Create a vertical perspective coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None + ) + perspective_point_height = getattr( + cf_grid_var, CF_ATTR_GRID_PERSPECTIVE_HEIGHT, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.VerticalPerspective( + latitude_of_projection_origin, + longitude_of_projection_origin, + perspective_point_height, + false_easting, + false_northing, + ellipsoid, + ) + + return cs + + +################################################################################ +def build_geostationary_coordinate_system(engine, cf_grid_var): + """ + Create a geostationary coordinate system from the CF-netCDF + grid mapping variable. + + """ + major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var) + + latitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None + ) + longitude_of_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_LON_OF_PROJ_ORIGIN, None + ) + perspective_point_height = getattr( + cf_grid_var, CF_ATTR_GRID_PERSPECTIVE_HEIGHT, None + ) + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + sweep_angle_axis = getattr( + cf_grid_var, CF_ATTR_GRID_SWEEP_ANGLE_AXIS, None + ) + + ellipsoid = None + if ( + major is not None + or minor is not None + or inverse_flattening is not None + ): + ellipsoid = iris.coord_systems.GeogCS(major, minor, inverse_flattening) + + cs = iris.coord_systems.Geostationary( + latitude_of_projection_origin, + longitude_of_projection_origin, + perspective_point_height, + sweep_angle_axis, + false_easting, + false_northing, + ellipsoid, + ) + + return cs + + +################################################################################ +def get_attr_units(cf_var, attributes): + attr_units = getattr(cf_var, CF_ATTR_UNITS, UNKNOWN_UNIT_STRING) + if not attr_units: + attr_units = UNKNOWN_UNIT_STRING + + # Sanitise lat/lon units. + if attr_units in UD_UNITS_LAT or attr_units in UD_UNITS_LON: + attr_units = "degrees" + + # Graceful loading of invalid units. + try: + cf_units.as_unit(attr_units) + except ValueError: + # Using converted unicode message. Can be reverted with Python 3. + msg = "Ignoring netCDF variable {!r} invalid units {!r}".format( + cf_var.cf_name, attr_units + ) + warnings.warn(msg) + attributes["invalid_units"] = attr_units + attr_units = UNKNOWN_UNIT_STRING + + if np.issubdtype(cf_var.dtype, np.str_): + attr_units = NO_UNIT_STRING + + if any( + hasattr(cf_var.cf_data, name) + for name in ("flag_values", "flag_masks", "flag_meanings") + ): + attr_units = cf_units._NO_UNIT_STRING + + # Get any assoicated calendar for a time reference coordinate. + if cf_units.as_unit(attr_units).is_time_reference(): + attr_calendar = getattr(cf_var, CF_ATTR_CALENDAR, None) + + if attr_calendar: + attr_units = cf_units.Unit(attr_units, calendar=attr_calendar) + + return attr_units + + +################################################################################ +def get_names(cf_coord_var, coord_name, attributes): + """Determine the standard_name, long_name and var_name attributes.""" + + standard_name = getattr(cf_coord_var, CF_ATTR_STD_NAME, None) + long_name = getattr(cf_coord_var, CF_ATTR_LONG_NAME, None) + cf_name = str(cf_coord_var.cf_name) + + if standard_name is not None: + try: + standard_name = _get_valid_standard_name(standard_name) + except ValueError: + if long_name is not None: + attributes["invalid_standard_name"] = standard_name + if coord_name is not None: + standard_name = coord_name + else: + standard_name = None + else: + if coord_name is not None: + attributes["invalid_standard_name"] = standard_name + standard_name = coord_name + else: + standard_name = None + + else: + if coord_name is not None: + standard_name = coord_name + + # Last attempt to set the standard name to something meaningful. + if standard_name is None: + if cf_name in iris.std_names.STD_NAMES: + standard_name = cf_name + + return (standard_name, long_name, cf_name) + + +################################################################################ +def get_cf_bounds_var(cf_coord_var): + """ + Return the CF variable representing the bounds of a coordinate + variable. + + """ + attr_bounds = getattr(cf_coord_var, CF_ATTR_BOUNDS, None) + attr_climatology = getattr(cf_coord_var, CF_ATTR_CLIMATOLOGY, None) + + # Determine bounds, prefering standard bounds over climatology. + # NB. No need to raise a warning if the bounds/climatology + # variable is missing, as that will already have been done by + # iris.fileformats.cf. + cf_bounds_var = None + climatological = False + if attr_bounds is not None: + bounds_vars = cf_coord_var.cf_group.bounds + if attr_bounds in bounds_vars: + cf_bounds_var = bounds_vars[attr_bounds] + elif attr_climatology is not None: + climatology_vars = cf_coord_var.cf_group.climatology + if attr_climatology in climatology_vars: + cf_bounds_var = climatology_vars[attr_climatology] + climatological = True + + if attr_bounds is not None and attr_climatology is not None: + warnings.warn( + "Ignoring climatology in favour of bounds attribute " + "on NetCDF variable {!r}.".format(cf_coord_var.cf_name) + ) + + return cf_bounds_var, climatological + + +################################################################################ +def reorder_bounds_data(bounds_data, cf_bounds_var, cf_coord_var): + """ + Return a bounds_data array with the vertex dimension as the most + rapidly varying. + + .. note:: + + This function assumes the dimension names of the coordinate + variable match those of the bounds variable in order to determine + which is the vertex dimension. + + + """ + vertex_dim_names = set(cf_bounds_var.dimensions).difference( + cf_coord_var.dimensions + ) + if len(vertex_dim_names) != 1: + msg = ( + "Too many dimension names differ between coordinate " + "variable {!r} and the bounds variable {!r}. " + "Expected 1, got {}." + ) + raise ValueError( + msg.format( + str(cf_coord_var.cf_name), + str(cf_bounds_var.cf_name), + len(vertex_dim_names), + ) + ) + vertex_dim = cf_bounds_var.dimensions.index(*vertex_dim_names) + bounds_data = np.rollaxis( + bounds_data.view(), vertex_dim, len(bounds_data.shape) + ) + return bounds_data + + +################################################################################ +def build_dimension_coordinate( + engine, cf_coord_var, coord_name=None, coord_system=None +): + """Create a dimension coordinate (DimCoord) and add it to the cube.""" + + cf_var = engine.cf_var + cube = engine.cube + attributes = {} + + attr_units = get_attr_units(cf_coord_var, attributes) + points_data = cf_coord_var[:] + # Gracefully fill points masked array. + if ma.is_masked(points_data): + points_data = ma.filled(points_data) + msg = "Gracefully filling {!r} dimension coordinate masked points" + warnings.warn(msg.format(str(cf_coord_var.cf_name))) + + # Get any coordinate bounds. + cf_bounds_var, climatological = get_cf_bounds_var(cf_coord_var) + if cf_bounds_var is not None: + bounds_data = cf_bounds_var[:] + # Gracefully fill bounds masked array. + if ma.is_masked(bounds_data): + bounds_data = ma.filled(bounds_data) + msg = "Gracefully filling {!r} dimension coordinate masked bounds" + warnings.warn(msg.format(str(cf_coord_var.cf_name))) + # Handle transposed bounds where the vertex dimension is not + # the last one. Test based on shape to support different + # dimension names. + if cf_bounds_var.shape[:-1] != cf_coord_var.shape: + bounds_data = reorder_bounds_data( + bounds_data, cf_bounds_var, cf_coord_var + ) + else: + bounds_data = None + + # Determine whether the coordinate is circular. + circular = False + if ( + points_data.ndim == 1 + and coord_name in [CF_VALUE_STD_NAME_LON, CF_VALUE_STD_NAME_GRID_LON] + and cf_units.Unit(attr_units) + in [cf_units.Unit("radians"), cf_units.Unit("degrees")] + ): + modulus_value = cf_units.Unit(attr_units).modulus + circular = iris.util._is_circular( + points_data, modulus_value, bounds=bounds_data + ) + + # Determine the name of the dimension/s shared between the CF-netCDF data variable + # and the coordinate being built. + common_dims = [ + dim for dim in cf_coord_var.dimensions if dim in cf_var.dimensions + ] + data_dims = None + if common_dims: + # Calculate the offset of each common dimension. + data_dims = [cf_var.dimensions.index(dim) for dim in common_dims] + + # Determine the standard_name, long_name and var_name + standard_name, long_name, var_name = get_names( + cf_coord_var, coord_name, attributes + ) + + # Create the coordinate. + try: + coord = iris.coords.DimCoord( + points_data, + standard_name=standard_name, + long_name=long_name, + var_name=var_name, + units=attr_units, + bounds=bounds_data, + attributes=attributes, + coord_system=coord_system, + circular=circular, + climatological=climatological, + ) + except ValueError as e_msg: + # Attempt graceful loading. + coord = iris.coords.AuxCoord( + points_data, + standard_name=standard_name, + long_name=long_name, + var_name=var_name, + units=attr_units, + bounds=bounds_data, + attributes=attributes, + coord_system=coord_system, + climatological=climatological, + ) + cube.add_aux_coord(coord, data_dims) + msg = ( + "Failed to create {name!r} dimension coordinate: {error}\n" + "Gracefully creating {name!r} auxiliary coordinate instead." + ) + warnings.warn(msg.format(name=str(cf_coord_var.cf_name), error=e_msg)) + else: + # Add the dimension coordinate to the cube. + if data_dims: + cube.add_dim_coord(coord, data_dims) + else: + # Scalar coords are placed in the aux_coords container. + cube.add_aux_coord(coord, data_dims) + + # Update the coordinate to CF-netCDF variable mapping. + engine.cube_parts["coordinates"].append((coord, cf_coord_var.cf_name)) + + +################################################################################ +def build_auxiliary_coordinate( + engine, cf_coord_var, coord_name=None, coord_system=None +): + """Create an auxiliary coordinate (AuxCoord) and add it to the cube.""" + + cf_var = engine.cf_var + cube = engine.cube + attributes = {} + + # Get units + attr_units = get_attr_units(cf_coord_var, attributes) + + # Get any coordinate point data. + if isinstance(cf_coord_var, cf.CFLabelVariable): + points_data = cf_coord_var.cf_label_data(cf_var) + else: + points_data = _get_cf_var_data(cf_coord_var, engine.filename) + + # Get any coordinate bounds. + cf_bounds_var, climatological = get_cf_bounds_var(cf_coord_var) + if cf_bounds_var is not None: + bounds_data = _get_cf_var_data(cf_bounds_var, engine.filename) + + # Handle transposed bounds where the vertex dimension is not + # the last one. Test based on shape to support different + # dimension names. + if cf_bounds_var.shape[:-1] != cf_coord_var.shape: + # Resolving the data to a numpy array (i.e. *not* masked) for + # compatibility with array creators (i.e. dask) + bounds_data = np.asarray(bounds_data) + bounds_data = reorder_bounds_data( + bounds_data, cf_bounds_var, cf_coord_var + ) + else: + bounds_data = None + + # Determine the name of the dimension/s shared between the CF-netCDF data variable + # and the coordinate being built. + common_dims = [ + dim for dim in cf_coord_var.dimensions if dim in cf_var.dimensions + ] + data_dims = None + if common_dims: + # Calculate the offset of each common dimension. + data_dims = [cf_var.dimensions.index(dim) for dim in common_dims] + + # Determine the standard_name, long_name and var_name + standard_name, long_name, var_name = get_names( + cf_coord_var, coord_name, attributes + ) + + # Create the coordinate + coord = iris.coords.AuxCoord( + points_data, + standard_name=standard_name, + long_name=long_name, + var_name=var_name, + units=attr_units, + bounds=bounds_data, + attributes=attributes, + coord_system=coord_system, + climatological=climatological, + ) + + # Add it to the cube + cube.add_aux_coord(coord, data_dims) + + # Make a list with names, stored on the engine, so we can find them all later. + engine.cube_parts["coordinates"].append((coord, cf_coord_var.cf_name)) + + +################################################################################ +def build_cell_measures(engine, cf_cm_var): + """Create a CellMeasure instance and add it to the cube.""" + cf_var = engine.cf_var + cube = engine.cube + attributes = {} + + # Get units + attr_units = get_attr_units(cf_cm_var, attributes) + + # Get (lazy) content array + data = _get_cf_var_data(cf_cm_var, engine.filename) + + # Determine the name of the dimension/s shared between the CF-netCDF data variable + # and the coordinate being built. + common_dims = [ + dim for dim in cf_cm_var.dimensions if dim in cf_var.dimensions + ] + data_dims = None + if common_dims: + # Calculate the offset of each common dimension. + data_dims = [cf_var.dimensions.index(dim) for dim in common_dims] + + # Determine the standard_name, long_name and var_name + standard_name, long_name, var_name = get_names(cf_cm_var, None, attributes) + + # Obtain the cf_measure. + measure = cf_cm_var.cf_measure + + # Create the CellMeasure + cell_measure = iris.coords.CellMeasure( + data, + standard_name=standard_name, + long_name=long_name, + var_name=var_name, + units=attr_units, + attributes=attributes, + measure=measure, + ) + + # Add it to the cube + cube.add_cell_measure(cell_measure, data_dims) + + # Make a list with names, stored on the engine, so we can find them all later. + engine.cube_parts["cell_measures"].append( + (cell_measure, cf_cm_var.cf_name) + ) + + +################################################################################ +def build_ancil_var(engine, cf_av_var): + """Create an AncillaryVariable instance and add it to the cube.""" + cf_var = engine.cf_var + cube = engine.cube + attributes = {} + + # Get units + attr_units = get_attr_units(cf_av_var, attributes) + + # Get (lazy) content array + data = _get_cf_var_data(cf_av_var, engine.filename) + + # Determine the name of the dimension/s shared between the CF-netCDF data variable + # and the AV being built. + common_dims = [ + dim for dim in cf_av_var.dimensions if dim in cf_var.dimensions + ] + data_dims = None + if common_dims: + # Calculate the offset of each common dimension. + data_dims = [cf_var.dimensions.index(dim) for dim in common_dims] + + # Determine the standard_name, long_name and var_name + standard_name, long_name, var_name = get_names(cf_av_var, None, attributes) + + # Create the AncillaryVariable + av = iris.coords.AncillaryVariable( + data, + standard_name=standard_name, + long_name=long_name, + var_name=var_name, + units=attr_units, + attributes=attributes, + ) + + # Add it to the cube + cube.add_ancillary_variable(av, data_dims) + + # Make a list with names, stored on the engine, so we can find them all later. + engine.cube_parts["ancillary_variables"].append((av, cf_av_var.cf_name)) + + +################################################################################ +def _is_lat_lon( + cf_var, ud_units, std_name, std_name_grid, axis_name, prefixes +): + """ + Determine whether the CF coordinate variable is a latitude/longitude variable. + + Ref: [CF] Section 4.1 Latitude Coordinate. + [CF] Section 4.2 Longitude Coordinate. + + """ + is_valid = False + attr_units = getattr(cf_var, CF_ATTR_UNITS, None) + + if attr_units is not None: + attr_units = attr_units.lower() + is_valid = attr_units in ud_units + + # Special case - Check for rotated pole. + if attr_units == "degrees": + attr_std_name = getattr(cf_var, CF_ATTR_STD_NAME, None) + if attr_std_name is not None: + is_valid = attr_std_name.lower() == std_name_grid + else: + is_valid = False + # TODO: check that this interpretation of axis is correct. + attr_axis = getattr(cf_var, CF_ATTR_AXIS, None) + if attr_axis is not None: + is_valid = attr_axis.lower() == axis_name + else: + # Alternative is to check standard_name or axis. + attr_std_name = getattr(cf_var, CF_ATTR_STD_NAME, None) + + if attr_std_name is not None: + attr_std_name = attr_std_name.lower() + is_valid = attr_std_name in [std_name, std_name_grid] + if not is_valid: + is_valid = any( + [attr_std_name.startswith(prefix) for prefix in prefixes] + ) + else: + attr_axis = getattr(cf_var, CF_ATTR_AXIS, None) + + if attr_axis is not None: + is_valid = attr_axis.lower() == axis_name + + return is_valid + + +################################################################################ +def is_latitude(engine, cf_name): + """Determine whether the CF coordinate variable is a latitude variable.""" + cf_var = engine.cf_var.cf_group[cf_name] + return _is_lat_lon( + cf_var, + UD_UNITS_LAT, + CF_VALUE_STD_NAME_LAT, + CF_VALUE_STD_NAME_GRID_LAT, + CF_VALUE_AXIS_Y, + ["lat", "rlat"], + ) + + +################################################################################ +def is_longitude(engine, cf_name): + """Determine whether the CF coordinate variable is a longitude variable.""" + cf_var = engine.cf_var.cf_group[cf_name] + return _is_lat_lon( + cf_var, + UD_UNITS_LON, + CF_VALUE_STD_NAME_LON, + CF_VALUE_STD_NAME_GRID_LON, + CF_VALUE_AXIS_X, + ["lon", "rlon"], + ) + + +################################################################################ +def is_projection_x_coordinate(engine, cf_name): + """ + Determine whether the CF coordinate variable is a + projection_x_coordinate variable. + + """ + cf_var = engine.cf_var.cf_group[cf_name] + attr_name = getattr(cf_var, CF_ATTR_STD_NAME, None) or getattr( + cf_var, CF_ATTR_LONG_NAME, None + ) + return attr_name == CF_VALUE_STD_NAME_PROJ_X + + +################################################################################ +def is_projection_y_coordinate(engine, cf_name): + """ + Determine whether the CF coordinate variable is a + projection_y_coordinate variable. + + """ + cf_var = engine.cf_var.cf_group[cf_name] + attr_name = getattr(cf_var, CF_ATTR_STD_NAME, None) or getattr( + cf_var, CF_ATTR_LONG_NAME, None + ) + return attr_name == CF_VALUE_STD_NAME_PROJ_Y + + +################################################################################ +def is_time(engine, cf_name): + """ + Determine whether the CF coordinate variable is a time variable. + + Ref: [CF] Section 4.4 Time Coordinate. + + """ + cf_var = engine.cf_var.cf_group[cf_name] + attr_units = getattr(cf_var, CF_ATTR_UNITS, None) + + attr_std_name = getattr(cf_var, CF_ATTR_STD_NAME, None) + attr_axis = getattr(cf_var, CF_ATTR_AXIS, "") + try: + is_time_reference = cf_units.Unit(attr_units or 1).is_time_reference() + except ValueError: + is_time_reference = False + + return is_time_reference and ( + attr_std_name == "time" or attr_axis.lower() == CF_VALUE_AXIS_T + ) + + +################################################################################ +def is_time_period(engine, cf_name): + """Determine whether the CF coordinate variable represents a time period.""" + is_valid = False + cf_var = engine.cf_var.cf_group[cf_name] + attr_units = getattr(cf_var, CF_ATTR_UNITS, None) + + if attr_units is not None: + try: + is_valid = cf_units.is_time(attr_units) + except ValueError: + is_valid = False + + return is_valid + + +################################################################################ +def is_grid_mapping(engine, cf_name, grid_mapping): + """Determine whether the CF grid mapping variable is of the appropriate type.""" + + is_valid = False + cf_var = engine.cf_var.cf_group[cf_name] + attr_mapping_name = getattr(cf_var, CF_ATTR_GRID_MAPPING_NAME, None) + + if attr_mapping_name is not None: + is_valid = attr_mapping_name.lower() == grid_mapping + + return is_valid + + +################################################################################ +def _is_rotated(engine, cf_name, cf_attr_value): + """Determine whether the CF coordinate variable is rotated.""" + + is_valid = False + cf_var = engine.cf_var.cf_group[cf_name] + attr_std_name = getattr(cf_var, CF_ATTR_STD_NAME, None) + + if attr_std_name is not None: + is_valid = attr_std_name.lower() == cf_attr_value + else: + attr_units = getattr(cf_var, CF_ATTR_UNITS, None) + if attr_units is not None: + is_valid = attr_units.lower() == "degrees" + + return is_valid + + +################################################################################ +def is_rotated_latitude(engine, cf_name): + """Determine whether the CF coodinate variable is rotated latitude.""" + return _is_rotated(engine, cf_name, CF_VALUE_STD_NAME_GRID_LAT) + + +############################################################################### +def is_rotated_longitude(engine, cf_name): + """Determine whether the CF coordinate variable is rotated longitude.""" + return _is_rotated(engine, cf_name, CF_VALUE_STD_NAME_GRID_LON) + + +################################################################################ +def has_supported_mercator_parameters(engine, cf_name): + """Determine whether the CF grid mapping variable has the supported + values for the parameters of the Mercator projection.""" + + is_valid = True + cf_grid_var = engine.cf_var.cf_group[cf_name] + + false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None) + false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None) + scale_factor_at_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None + ) + standard_parallel = getattr( + cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None + ) + + if false_easting is not None and false_easting != 0: + warnings.warn( + "False eastings other than 0.0 not yet supported " + "for Mercator projections" + ) + is_valid = False + if false_northing is not None and false_northing != 0: + warnings.warn( + "False northings other than 0.0 not yet supported " + "for Mercator projections" + ) + is_valid = False + if ( + scale_factor_at_projection_origin is not None + and scale_factor_at_projection_origin != 1 + ): + warnings.warn( + "Scale factors other than 1.0 not yet supported for " + "Mercator projections" + ) + is_valid = False + if standard_parallel is not None and standard_parallel != 0: + warnings.warn( + "Standard parallels other than 0.0 not yet " + "supported for Mercator projections" + ) + is_valid = False + + return is_valid + + +################################################################################ +def has_supported_stereographic_parameters(engine, cf_name): + """Determine whether the CF grid mapping variable has a value of 1.0 + for the scale_factor_at_projection_origin attribute.""" + + is_valid = True + cf_grid_var = engine.cf_var.cf_group[cf_name] + + scale_factor_at_projection_origin = getattr( + cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None + ) + + if ( + scale_factor_at_projection_origin is not None + and scale_factor_at_projection_origin != 1 + ): + warnings.warn( + "Scale factors other than 1.0 not yet supported for " + "stereographic projections" + ) + is_valid = False + + return is_valid diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 206f7526c6..37e6f39b41 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -384,7 +384,7 @@ def coord(self, name): return result -def _pyke_kb_engine(): +def _pyke_kb_engine_real(): """Return the PyKE knowledge engine for CF->cube conversion.""" pyke_dir = os.path.join(os.path.dirname(__file__), "_pyke_rules") @@ -415,6 +415,22 @@ def _pyke_kb_engine(): return engine +LOAD_PYKE = False +CHECK_BOTH = True + + +def _pyke_kb_engine(use_pyke=True): + """Return a knowledge engine, or replacement object.""" + if use_pyke: + engine = _pyke_kb_engine_real() + else: + # Deferred import to avoid circularity. + import iris.fileformats._nc_load_rules.engine as nonpyke_engine + + engine = nonpyke_engine.Engine() + return engine + + class NetCDFDataProxy: """A reference to the data payload of a single NetCDF file variable.""" @@ -584,6 +600,17 @@ def _get_cf_var_data(cf_var, filename): return as_lazy_data(proxy, chunks=chunks) +class OrderedAddableList(list): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self._n_add = 0 + + def add(self, msg): + self._n_add += 1 + n_add = self._n_add + self.append(f"#{n_add:03d} : {msg}") + + def _load_cube(engine, cf, cf_var, filename): """Create the cube associated with the CF-netCDF data variable.""" data = _get_cf_var_data(cf_var, filename) @@ -597,7 +624,7 @@ def _load_cube(engine, cf, cf_var, filename): engine.cube = cube engine.cube_parts = {} engine.requires = {} - engine.rule_triggered = set() + engine.rule_triggered = OrderedAddableList() # set() engine.filename = filename # Assert any case-specific facts. @@ -790,7 +817,10 @@ def load_cubes(filenames, callback=None): """ # Initialise the pyke inference engine. - engine = _pyke_kb_engine() + if LOAD_PYKE or CHECK_BOTH: + pyke_engine = _pyke_kb_engine(use_pyke=True) + if not LOAD_PYKE or CHECK_BOTH: + nonpyke_engine = _pyke_kb_engine(use_pyke=False) if isinstance(filenames, str): filenames = [filenames] @@ -804,14 +834,59 @@ def load_cubes(filenames, callback=None): cf.cf_group.promoted.values() ) for cf_var in data_variables: - cube = _load_cube(engine, cf, cf_var, filename) - # Process any associated formula terms and attach - # the corresponding AuxCoordFactory. - try: - _load_aux_factory(engine, cube) - except ValueError as e: - warnings.warn("{}".format(e)) + def one_cube(engine, show_warnings=True): + cube = _load_cube(engine, cf, cf_var, filename) + + # Process any associated formula terms and attach + # the corresponding AuxCoordFactory. + try: + _load_aux_factory(engine, cube) + except ValueError as e: + if show_warnings: + warnings.warn("{}".format(e)) + + return cube + + if LOAD_PYKE or CHECK_BOTH: + # Make pyke cube if needed, show warnings if "primary" cube + cube_pyke = one_cube(pyke_engine, show_warnings=LOAD_PYKE) + if not LOAD_PYKE or CHECK_BOTH: + # Make nonpyke cube if needed, show warnings if "primary" cube + cube_nonpyke = one_cube( + nonpyke_engine, show_warnings=not LOAD_PYKE + ) + + if not LOAD_PYKE or CHECK_BOTH: + msgs = [] + if not LOAD_PYKE: + msgs.append("NONPYKE-LOAD") + warnings.warn("((NONPYKE-LOAD))") + if CHECK_BOTH: + msgs.append("PYKE-SAMECHECK") + msg = f"(({': '.join(msgs)}))" + print(msg) + + if CHECK_BOTH: + + def unmask_cube(cube): + cube = cube.copy() + if isinstance(cube.data, np.ma.MaskedArray): + cube.data = cube.data.filled(0) + return cube + + if unmask_cube(cube_nonpyke) != unmask_cube(cube_pyke): + warnings.warn("PYKE-SAMECHECK : unmasked cubes == failed.") + full_test = cube_nonpyke.copy() == cube_pyke.copy() + if full_test not in (True, np.ma.masked): + # NOTE: this one can fail: non-bool answer if masked ? + warnings.warn("PYKE-SAMECHECK : copied cubes == failed.") + + # Select "primary" cube to return. + if LOAD_PYKE: + cube = cube_pyke + else: + cube = cube_nonpyke # Perform any user registered callback function. cube = iris.io.run_callback(callback, cube, cf_var, filename) diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/__init__.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/__init__.py new file mode 100644 index 0000000000..8bc429a906 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/__init__.py @@ -0,0 +1,6 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +"""Unit tests for the `iris.fileformats.netcdf._load_cube` function.""" diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/__init__.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/__init__.py new file mode 100644 index 0000000000..dd9f843783 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/__init__.py @@ -0,0 +1,253 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for the engine.activate() call within the +`iris.fileformats.netcdf._load_cube` function. + +For now, these tests are designed to function with **either** the "old" +Pyke-rules implementation in :mod:`iris.fileformats._pyke_rules`, **or** the +"new" :mod:`iris.fileformats._nc_load_rules`. +Both of those supply an "engine" with an "activate" method + -- at least for now : may be simplified in future. + +""" +from pathlib import Path +import shutil +import subprocess +import tempfile + +import numpy as np + +from iris.fileformats.cf import CFReader +import iris.fileformats.netcdf +from iris.fileformats.netcdf import _load_cube +import iris.fileformats._nc_load_rules.engine + +""" +Notes on testing method. + +IN cf : "def _load_cube(engine, cf, cf_var, filename)" +WHERE: + - engine is a :class:`pyke.knowledge_engine.engine` + -- **OR** :class:`iris.fileformats._nc_load_rules.engine.Engine` + - cf is a CFReader + - cf_var is a CFDAtaVariable + +As it's hard to construct a suitable CFReader from scratch, it would seem +simpler (for now) to use an ACTUAL FILE. +Likewise, the easiest approach to that is with CDL and "ncgen". +To do this, we need a test "fixture" that can create suitable test files in a +temporary directory. + +""" + + +class Mixin__nc_load_actions: + """ + Class to make testcases for rules or actions code, and check results. + + Defines standard setUpClass/tearDownClass methods, to create a temporary + directory for intermediate files. + NOTE: owing to peculiarities of unittest, these must be explicitly called + from a setUpClass/tearDownClass within the 'final' inheritor, i.e. the + actual Test_XXX class which also inherits unittest.TestCase. + + Testcases are manufactured by the '_make_testcase_cdl' method. + These are based on a 'standard simple latlon grid' example. + Various kwargs control variations on this. + + The 'run_testcase' method takes the '_make_testcase_cdl' kwargs and makes + a result cube (by: producing cdl, converting to netcdf, and loading). + + The 'check_result' method performs various checks on the result, with + kwargs controlling the expected properties to be tested against. + This usage is *also* based on the 'standard simple latlon grid' example, + the kwargs specify expected differences from that. + + Can also test with either the Pyke(rules) or non-Pyke (actions) + implementations (for now). + + """ + + # + # "global" test settings + # + + # whether to test 'rules' or 'actions' implementations + # TODO: remove when Pyke is gone + use_pyke = True + use_pyke_override = False + debug_confirm_which = True + + # whether to output various debug info + # TODO: ?possibly? remove when development is complete + debug = False + + # whether to perform action in both ways and compare results. + compare_pyke_nonpyke = True + + @classmethod + def setUpClass(cls): + # # Control which testing method we are applying. + # Create a temp directory for temp files. + cls.temp_dirpath = Path(tempfile.mkdtemp()) + + @classmethod + 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): + """ + Load the 'phenom' data variable in a CDL testcase, as a cube. + + Using ncgen, CFReader and the _load_cube call. + Can use a genuine Pyke engine, or the actions mimic engine, + selected by `self.use_pyke`. + + """ + # Write the CDL to a file. + with open(cdl_path, "w") as f_out: + f_out.write(cdl_string) + + # Create a netCDF file from the CDL file. + command = "ncgen -o {} {}".format(nc_path, cdl_path) + subprocess.check_call(command, shell=True) + + # Simulate the inner part of the file reading process. + cf = CFReader(nc_path) + # Grab a data variable : FOR NOW, should be only 1 + cf_var = list(cf.cf_group.data_variables.values())[0] + cf_var = cf.cf_group.data_variables["phenom"] + + use_pyke = self.use_pyke + if self.use_pyke_override is not None: + use_pyke = self.use_pyke_override + do_pyke = use_pyke or self.compare_pyke_nonpyke + do_nonpyke = not use_pyke or self.compare_pyke_nonpyke + if do_pyke: + pyke_engine = iris.fileformats.netcdf._pyke_kb_engine_real() + if do_nonpyke: + nonpyke_engine = iris.fileformats._nc_load_rules.engine.Engine() + + iris.fileformats.netcdf.DEBUG = self.debug + + # Call the main translation function to load a single cube. + def load_single_cube(engine): + # _load_cube establishes per-cube facts, activates rules and + # produces an actual cube. + cube = _load_cube(engine, cf, cf_var, nc_path) + + # Also Record, on the cubes, which hybrid coord elements were identified + # by the rules operation. + # Unlike the other translations, _load_cube does *not* convert this + # information into actual cube elements. That is instead done by + # `iris.fileformats.netcdf._load_aux_factory`. + # For rules testing, it is anyway more convenient to deal with the raw + # data, as each factory type has different validity requirements to + # build it, and none of that is relevant to the rules operation. + cube._formula_type_name = engine.requires.get("formula_type") + cube._formula_terms_byname = engine.requires.get("formula_terms") + + return cube + + if do_pyke: + pyke_cube = load_single_cube(pyke_engine) + if do_nonpyke: + nonpyke_cube = load_single_cube(nonpyke_engine) + + # If requested, directly compare the pyke and non-pyke outputs. + if self.compare_pyke_nonpyke: + # Compare the loaded cubes from both engines. + # print("\nPYKE-NONPYKE COMPARE") + + # First zap cube-data, as masked data does not compare well. + def unmask_cube(cube): + # preserve the original, we're going to realise.. + cube = cube.copy() + if isinstance(cube.data, np.ma.MaskedArray): + cube.data = cube.data.filled(0) + return cube + + pyke_cube_copy = unmask_cube(pyke_cube) + nonpyke_cube_copy = unmask_cube(nonpyke_cube) + if self.debug: + if nonpyke_cube_copy != pyke_cube_copy: + + def show_cube(cube): + result = str(cube) + result += "\n--coords--" + for coord in cube.coords(): + result += "\n " + str(coord) + result += "\n--attributes--" + if not cube.attributes: + result += "\n (none)" + else: + for key, value in cube.attributes.items(): + result += f"\n {key}: {value}" + return result + + print("\nPyke/nonpyke mismatch.") + print("Pyke cube:\n----") + print(show_cube(pyke_cube)) + print() + print("NONPyke cube:\n----") + print(show_cube(nonpyke_cube)) + print("") + else: + self.assertEqual(pyke_cube_copy, nonpyke_cube_copy) + + # Return the right thing, whether we did 'both' or not + if use_pyke: + if self.debug_confirm_which: + print("(PYKE_TESTED)") + result_cube = pyke_cube + else: + if self.debug_confirm_which: + print("(NONPYKE_TESTED)") + result_cube = nonpyke_cube + + # Always returns a single cube. + return result_cube + + def run_testcase(self, warning=None, **testcase_kwargs): + """ + Run a testcase with chosen options, returning a test cube. + + The kwargs apply to the '_make_testcase_cdl' method. + + """ + cdl_path = str(self.temp_dirpath / "test.cdl") + nc_path = cdl_path.replace(".cdl", ".nc") + + cdl_string = self._make_testcase_cdl(**testcase_kwargs) + if self.debug: + print("CDL file content:") + print(cdl_string) + print("------\n") + + if warning is None: + context = self.assertNoWarningsRegexp() + else: + context = self.assertWarnsRegexp(warning) + with context: + cube = self.load_cube_from_cdl(cdl_string, cdl_path, nc_path) + + if self.debug: + print("\nCube:") + print(cube) + print("") + return cube + + def _make_testcase_cdl(self, **kwargs): + """Make a testcase CDL string.""" + # Override for specific uses... + raise NotImplementedError() + + def check_result(self, cube, **kwargs): + """Test a result cube.""" + # Override for specific uses... + raise NotImplementedError() diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__additional.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__additional.py new file mode 100644 index 0000000000..957c736501 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__additional.py @@ -0,0 +1,57 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for the engine.activate() call within the +`iris.fileformats.netcdf._load_cube` function. + +For now, these tests are designed to function with **either** the "old" +Pyke-rules implementation in :mod:`iris.fileformats._pyke_rules`, **or** the +"new" :mod:`iris.fileformats._nc_load_rules`. +Both of those supply an "engine" with an "activate" method + -- at least for now : may be simplified in future. + +""" +import iris.tests as tests + + +from iris.tests.unit.fileformats.netcdf.load_cube.load_cube__activate.test__grid_mappings import ( + Mixin__grid_mapping, +) + + +class Test__additional(Mixin__grid_mapping, tests.IrisTest): + # Run grid-mapping tests with non-Pyke (actions) + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + def test_nondim_lats(self): + # Check what happens when values don't allow a coord to be dim-coord. + # + # 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: + # in terms of rule triggers, this is not distinct from a normal case + # - but the latitude is now an aux-coord. + warning = "must be.* monotonic" + result = self.run_testcase(warning=warning, yco_values=[0.0, 0.0]) + self.check_result(result, yco_is_aux=True) + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__grid_mappings.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__grid_mappings.py new file mode 100644 index 0000000000..5a11a2cc88 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__grid_mappings.py @@ -0,0 +1,869 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for the engine.activate() call within the +`iris.fileformats.netcdf._load_cube` function. + +Here, *specifically* testcases relating to grid-mappings and dim-coords. + +""" +import iris.tests as tests + +import iris.coord_systems as ics +import iris.fileformats._nc_load_rules.helpers as hh + +from iris.tests.unit.fileformats.netcdf.load_cube.load_cube__activate import ( + Mixin__nc_load_actions, +) + + +class Mixin__grid_mapping(Mixin__nc_load_actions): + # Testcase support routines for testing translation of grid-mappings + def _make_testcase_cdl( + self, + latitude_units=None, + gridmapvar_name=None, + gridmapvar_mappropertyname=None, + mapping_missingradius=False, + mapping_type_name=None, + mapping_scalefactor=None, + yco_values=None, + xco_name=None, + yco_name=None, + xco_units=None, + yco_units=None, + xco_is_dim=True, + yco_is_dim=True, + ): + """ + Create a CDL string for a testcase. + + This is the "master" routine for creating all our testcases. + Kwarg options modify a simple default testcase with a latlon grid. + The routine handles the various testcase options and their possible + interactions. This includes knowing what extra changes are required + to support different grid-mapping types (for example). + + """ + # The grid-mapping options are standard-latlon, rotated, or non-latlon. + # This affects names+units of the X and Y coords. + # We don't have an option to *not* include a grid-mapping variable, but + # we can mimic a missing grid-mapping by changing the varname from that + # which the data-variable refers to, with "gridmapvar_name=xxx". + # Likewise, an invalid (unrecognised) grid-mapping can be mimicked by + # selecting an unkown 'grid_mapping_name' property, with + # "gridmapvar_mappropertyname=xxx". + if mapping_type_name is None: + # Default grid-mapping and coords are standard lat-lon. + mapping_type_name = hh.CF_GRID_MAPPING_LAT_LON + xco_name_default = hh.CF_VALUE_STD_NAME_LON + yco_name_default = hh.CF_VALUE_STD_NAME_LAT + xco_units_default = "degrees_east" + # Special kwarg overrides some of the values. + if latitude_units is None: + yco_units_default = "degrees_north" + else: + # Override the latitude units (to invalidate). + yco_units_default = latitude_units + + elif mapping_type_name == hh.CF_GRID_MAPPING_ROTATED_LAT_LON: + # Rotated lat-lon coordinates. + xco_name_default = hh.CF_VALUE_STD_NAME_GRID_LON + yco_name_default = hh.CF_VALUE_STD_NAME_GRID_LAT + xco_units_default = "degrees" + yco_units_default = "degrees" + + else: + # General non-latlon coordinates + # Exactly which depends on the grid_mapping name. + xco_name_default = hh.CF_VALUE_STD_NAME_PROJ_X + yco_name_default = hh.CF_VALUE_STD_NAME_PROJ_Y + xco_units_default = "m" + yco_units_default = "m" + + # Options can override coord (standard) names and units. + if xco_name is None: + xco_name = xco_name_default + if yco_name is None: + yco_name = yco_name_default + if xco_units is None: + xco_units = xco_units_default + if yco_units is None: + yco_units = yco_units_default + + phenom_auxcoord_names = [] + if xco_is_dim: + # xdim has same name as xco, making xco a dim-coord + xdim_name = "xco" + else: + # use alternate dim-name, and put xco on the 'coords' list + # This makes the X coord an aux-coord + xdim_name = "xdim_altname" + phenom_auxcoord_names.append("xco") + if yco_is_dim: + # ydim has same name as yco, making yco a dim-coord + ydim_name = "yco" # This makes the Y coord a dim-coord + else: + # use alternate dim-name, and put yco on the 'coords' list + # This makes the Y coord an aux-coord + ydim_name = "ydim_altname" + phenom_auxcoord_names.append("yco") + + # Build a 'phenom:coords' string if needed. + if phenom_auxcoord_names: + phenom_coords_string = " ".join(phenom_auxcoord_names) + phenom_coords_string = f""" + phenom:coordinates = "{phenom_coords_string}" ; +""" + else: + phenom_coords_string = "" + + grid_mapping_name = "grid" + # Options can override the gridvar name and properties. + g_varname = gridmapvar_name + g_mapname = gridmapvar_mappropertyname + if g_varname is None: + g_varname = grid_mapping_name + if g_mapname is None: + # If you change this, it is no longer a valid grid-mapping var. + g_mapname = "grid_mapping_name" + + # Omit the earth radius, if requested. + if mapping_missingradius: + g_radius_string = "" + else: + g_radius_string = f"{g_varname}:earth_radius = 6.e6 ;" + g_string = f""" + int {g_varname} ; + {g_varname}:{g_mapname} = "{mapping_type_name}"; + {g_radius_string} + """ + + # Add a specified scale-factor, if requested. + if mapping_scalefactor is not None: + # Add a specific scale-factor term to the grid mapping. + # (Non-unity scale is not supported for Mercator/Stereographic). + sfapo_name = hh.CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN + g_string += f""" + {g_varname}:{sfapo_name} = {mapping_scalefactor} ; + """ + + # + # Add various additional (minimal) required properties for different + # grid mapping types. + # + + # Those which require 'latitude of projection origin' + if mapping_type_name in ( + hh.CF_GRID_MAPPING_TRANSVERSE, + hh.CF_GRID_MAPPING_STEREO, + hh.CF_GRID_MAPPING_GEOSTATIONARY, + hh.CF_GRID_MAPPING_VERTICAL, + ): + latpo_name = hh.CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN + g_string += f""" + {g_varname}:{latpo_name} = 0.0 ; + """ + # Those which require 'longitude of projection origin' + if mapping_type_name in ( + hh.CF_GRID_MAPPING_STEREO, + hh.CF_GRID_MAPPING_GEOSTATIONARY, + hh.CF_GRID_MAPPING_VERTICAL, + ): + lonpo_name = hh.CF_ATTR_GRID_LON_OF_PROJ_ORIGIN + g_string += f""" + {g_varname}:{lonpo_name} = 0.0 ; + """ + # Those which require 'longitude of central meridian' + if mapping_type_name in (hh.CF_GRID_MAPPING_TRANSVERSE,): + latcm_name = hh.CF_ATTR_GRID_LON_OF_CENT_MERIDIAN + g_string += f""" + {g_varname}:{latcm_name} = 0.0 ; + """ + # Those which require 'perspective point height' + if mapping_type_name in ( + hh.CF_GRID_MAPPING_VERTICAL, + hh.CF_GRID_MAPPING_GEOSTATIONARY, + ): + pph_name = hh.CF_ATTR_GRID_PERSPECTIVE_HEIGHT + g_string += f""" + {g_varname}:{pph_name} = 600000.0 ; + """ + # Those which require 'sweep angle axis' + if mapping_type_name in (hh.CF_GRID_MAPPING_GEOSTATIONARY,): + saa_name = hh.CF_ATTR_GRID_SWEEP_ANGLE_AXIS + g_string += f""" + {g_varname}:{saa_name} = "y" ; + """ + + # y-coord values + if yco_values is None: + yco_values = [10.0, 20.0] + yco_value_strings = [str(val) for val in yco_values] + yco_values_string = ", ".join(yco_value_strings) + + # Construct the total CDL string + cdl_string = f""" + netcdf test {{ + dimensions: + {ydim_name} = 2 ; + {xdim_name} = 3 ; + variables: + double phenom({ydim_name}, {xdim_name}) ; + phenom:standard_name = "air_temperature" ; + phenom:units = "K" ; + phenom:grid_mapping = "grid" ; +{phenom_coords_string} + double yco({ydim_name}) ; + yco:axis = "Y" ; + yco:units = "{yco_units}" ; + yco:standard_name = "{yco_name}" ; + double xco({xdim_name}) ; + xco:axis = "X" ; + xco:units = "{xco_units}" ; + xco:standard_name = "{xco_name}" ; + {g_string} + data: + yco = {yco_values_string} ; + xco = 100., 110., 120. ; + }} + """ + return cdl_string + + def check_result( + self, + cube, + cube_cstype=None, + cube_no_cs=False, + cube_no_xycoords=False, + xco_no_cs=False, # N.B. no effect if cube_no_cs is True + yco_no_cs=False, # N.B. no effect if cube_no_cs is True + xco_is_aux=False, + yco_is_aux=False, + xco_stdname=True, + yco_stdname=True, + ): + """ + Check key properties of a result cube. + + Various options control the expected things which are tested. + """ + self.assertEqual(cube.standard_name, "air_temperature") + self.assertEqual(cube.var_name, "phenom") + + x_coords = cube.coords(dimensions=(1,)) + y_coords = cube.coords(dimensions=(0,)) + expected_dim_coords = [] + expected_aux_coords = [] + if yco_is_aux: + expected_aux_coords += y_coords + else: + expected_dim_coords += y_coords + if xco_is_aux: + expected_aux_coords += x_coords + else: + expected_dim_coords += x_coords + + self.assertEqual( + set(expected_dim_coords), set(cube.coords(dim_coords=True)) + ) + if cube_no_xycoords: + self.assertEqual(expected_dim_coords, []) + x_coord = None + y_coord = None + else: + self.assertEqual(len(x_coords), 1) + (x_coord,) = x_coords + self.assertEqual(len(y_coords), 1) + (y_coord,) = y_coords + + self.assertEqual( + set(expected_aux_coords), set(cube.coords(dim_coords=False)) + ) + + if x_coord: + if xco_stdname is None: + # no check + pass + elif xco_stdname is True: + self.assertIsNotNone(x_coord.standard_name) + elif xco_stdname is False: + self.assertIsNone(x_coord.standard_name) + else: + self.assertEqual(x_coord.standard_name, xco_stdname) + + if y_coord: + if yco_stdname is None: + # no check + pass + if yco_stdname is True: + self.assertIsNotNone(y_coord.standard_name) + elif yco_stdname is False: + self.assertIsNone(y_coord.standard_name) + else: + self.assertEqual(y_coord.standard_name, yco_stdname) + + cube_cs = cube.coord_system() + if cube_no_xycoords: + yco_cs = None + xco_cs = None + else: + yco_cs = y_coord.coord_system + xco_cs = x_coord.coord_system + if cube_no_cs: + self.assertIsNone(cube_cs) + self.assertIsNone(yco_cs) + self.assertIsNone(xco_cs) + else: + if cube_cstype is not None: + self.assertIsInstance(cube_cs, cube_cstype) + if xco_no_cs: + self.assertIsNone(xco_cs) + else: + self.assertEqual(xco_cs, cube_cs) + if yco_no_cs: + self.assertIsNone(yco_cs) + else: + self.assertEqual(yco_cs, cube_cs) + + +class Mixin__grid_mapping__tests(Mixin__grid_mapping): + # Various testcases for translation of grid-mappings + + def test_basic_latlon(self): + # A basic reference example with a lat-long grid. + # 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: regular latlon + # dim-coords: lat+lon + # coords-build: standard latlon coords (with latlon coord-system) + result = self.run_testcase() + self.check_result(result) + + def test_missing_latlon_radius(self): + # Lat-long with a missing earth-radius causes an error. + # One of very few cases where activation may encounter an error. + # N.B. doesn't really test rules-activation, but maybe worth doing. + # (no rules trigger) + with self.assertRaisesRegex(ValueError, "No ellipsoid"): + self.run_testcase(mapping_missingradius=True) + + def test_bad_gridmapping_nameproperty(self): + # Fix the 'grid' var so it does not register as a grid-mapping. + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_latitude + # 003 : fc_provides_coordinate_longitude + # 004 : fc_build_coordinate_latitude_nocs + # 005 : fc_build_coordinate_longitude_nocs + # Notes: + # grid-mapping: NONE + # dim-coords: lat+lon + # coords-build: latlon coords NO coord-system + result = self.run_testcase(gridmapvar_mappropertyname="mappy") + self.check_result(result, cube_no_cs=True) + + def test_latlon_bad_gridmapping_varname(self): + # rename the grid-mapping variable so it is effectively 'missing' + # (I.E. the var named in "data-variable:grid_mapping" does not exist). + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_latitude + # 003 : fc_provides_coordinate_longitude + # 004 : fc_build_coordinate_latitude_nocs + # 005 : fc_build_coordinate_longitude_nocs + # Notes: + # no coord-system + # all the same as test_bad_gridmapping_nameproperty + warning = "Missing.*grid mapping variable 'grid'" + result = self.run_testcase(warning=warning, gridmapvar_name="grid_2") + self.check_result(result, cube_no_cs=True) + + def test_latlon_bad_latlon_unit(self): + # Check with bad latitude units : 'degrees' in place of 'degrees_north'. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_latitude_longitude + # 003 : fc_provides_coordinate_longitude + # 004 : fc_build_coordinate_longitude + # 005 : fc_default_coordinate + # Notes: + # grid-mapping: regular latlon + # dim-coords: + # x is regular longitude dim-coord + # y is 'default' coord ==> builds as an 'extra' dim-coord + # coords-build: + # x(lon) is regular latlon with coord-system + # y(lat) is a dim-coord, but NO coord-system + # = "fc_provides_coordinate_latitude" does not trigger, because it is + # not a valid latitude coordinate. + result = self.run_testcase(latitude_units="degrees") + self.check_result(result, yco_no_cs=True) + + def test_mapping_rotated(self): + # Test with rotated-latlon grid-mapping + # Distinct from both regular-latlon and non-latlon cases, as the + # coordinate standard names and units are different. + # (run_testcase/_make_testcase_cdl know how to handle that). + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_rotated_latitude_longitude + # 003 : fc_provides_coordinate_latitude + # 004 : fc_provides_coordinate_longitude + # 005 : fc_build_coordinate_latitude_rotated + # 006 : fc_build_coordinate_longitude_rotated + # Notes: + # grid-mapping: rotated lat-lon + # dim-coords: lat+lon + # coords-build: lat+lon coords ROTATED, with coord-system + # (rotated means different name + units) + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ROTATED_LAT_LON + ) + self.check_result(result, cube_cstype=ics.RotatedGeogCS) + + # + # All non-latlon coordinate systems ... + # These all have projection-x/y coordinates with units of metres. + # They all work the same way, except that Mercator/Stereographic have + # parameter checking routines that can fail. + # NOTE: various mapping types *require* certain addtional properties + # - without which an error will occur during translation. + # - run_testcase/_make_testcase_cdl know how to provide these + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_ + # 003 : fc_provides_projection_x_coordinate + # 004 : fc_provides_projection_y_coordinate + # 005 : fc_build_coordinate_projection_x_ + # 006 : fc_build_coordinate_projection_y_ + # Notes: + # grid-mapping: + # dim-coords: proj-x and -y + # coords-build: proj-x/-y_, with coord-system + + def test_mapping_albers(self): + result = self.run_testcase(mapping_type_name=hh.CF_GRID_MAPPING_ALBERS) + self.check_result(result, cube_cstype=ics.AlbersEqualArea) + + def test_mapping_geostationary(self): + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_GEOSTATIONARY + ) + self.check_result(result, cube_cstype=ics.Geostationary) + + def test_mapping_lambert_azimuthal(self): + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_LAMBERT_AZIMUTHAL + ) + self.check_result(result, cube_cstype=ics.LambertAzimuthalEqualArea) + + def test_mapping_lambert_conformal(self): + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_LAMBERT_CONFORMAL + ) + self.check_result(result, cube_cstype=ics.LambertConformal) + + def test_mapping_mercator(self): + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_MERCATOR + ) + self.check_result(result, cube_cstype=ics.Mercator) + + def test_mapping_mercator__fail_unsupported(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_projection_x_coordinate + # 003 : fc_provides_projection_y_coordinate + # Notes: + # grid-mapping: NONE + # dim-coords: proj-x and -y + # coords-build: NONE + # = NO coord-system + # = NO dim-coords built (cube has no coords) + # Set a non-unity scale factor, which mercator cannot handle. + warning = "not yet supported for Mercator" + result = self.run_testcase( + warning=warning, + mapping_type_name=hh.CF_GRID_MAPPING_MERCATOR, + mapping_scalefactor=2.0, + ) + self.check_result(result, cube_no_cs=True, cube_no_xycoords=True) + + def test_mapping_stereographic(self): + result = self.run_testcase(mapping_type_name=hh.CF_GRID_MAPPING_STEREO) + self.check_result(result, cube_cstype=ics.Stereographic) + + def test_mapping_stereographic__fail_unsupported(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_projection_x_coordinate + # 003 : fc_provides_projection_y_coordinate + # Notes: + # as for 'mercator__fail_unsupported', above + # = NO dim-coords built (cube has no coords) + # + # Set a non-unity scale factor, which stereo cannot handle. + warning = "not yet supported for stereographic" + result = self.run_testcase( + warning=warning, + mapping_type_name=hh.CF_GRID_MAPPING_STEREO, + mapping_scalefactor=2.0, + ) + self.check_result(result, cube_no_cs=True, cube_no_xycoords=True) + + def test_mapping_transverse_mercator(self): + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_TRANSVERSE + ) + self.check_result(result, cube_cstype=ics.TransverseMercator) + + def test_mapping_vertical_perspective(self): + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_VERTICAL + ) + self.check_result(result, cube_cstype=ics.VerticalPerspective) + + def test_mapping_unsupported(self): + # Use azimuthal, which is a real thing but we don't yet support it. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_projection_x_coordinate + # 003 : fc_provides_projection_y_coordinate + # NOTES: + # - there is no warning for this. + # TODO: perhaps there should be ? + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_AZIMUTHAL + ) + self.check_result(result, cube_no_cs=True, cube_no_xycoords=True) + + def test_mapping_undefined(self): + # Use a random, unknown "mapping type". + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_projection_x_coordinate + # 003 : fc_provides_projection_y_coordinate + # NOTES: + # - there is no warning for this. + # TODO: perhaps there should be ? + result = self.run_testcase(mapping_type_name="unknown") + self.check_result(result, cube_no_cs=True, cube_no_xycoords=True) + + # + # Cases where names(+units) of coords don't match the grid-mapping type + # Effectively, there are 9 possibilities for (latlon/rotated/projected) + # coords against (latlon/rotated/projected/missing) coord-systems. + # N.B. the results are not all the same ... + # + + def test_mapping__mismatch__latlon_coords_rotated_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_rotated_latitude_longitude + # 003 : fc_provides_coordinate_latitude + # 004 : fc_provides_coordinate_longitude + # NOTES: + # no build_coord triggers, as it requires the correct mapping type + # so no dim-coords at all in this case + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ROTATED_LAT_LON, + xco_name="longitude", + xco_units="degrees_east", + yco_name="latitude", + yco_units="degrees_north", + ) + self.check_result(result, cube_no_cs=True, cube_no_xycoords=True) + + def test_mapping__mismatch__latlon_coords_nonll_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_albers_equal_area + # 003 : fc_provides_coordinate_latitude + # 004 : fc_provides_coordinate_longitude + # 005 : fc_build_coordinate_latitude_nocs + # 006 : fc_build_coordinate_longitude_nocs + # NOTES: + # build_coord_XXX_cs triggers, requires NO latlon/rotated mapping + # - but a non-ll mapping is 'ok'. + # TODO: not really clear why this is right ? + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ALBERS, + xco_name="longitude", + xco_units="degrees_east", + yco_name="latitude", + yco_units="degrees_north", + ) + self.check_result(result, cube_no_cs=True) + + def test_mapping__mismatch__latlon_coords_missing_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_latitude + # 003 : fc_provides_coordinate_longitude + # 004 : fc_build_coordinate_latitude_nocs + # 005 : fc_build_coordinate_longitude_nocs + # NOTES: + # same as nonll, except *NO* grid-mapping is detected, + # - which makes no practical difference + warning = "Missing.*grid mapping variable 'grid'" + result = self.run_testcase( + warning=warning, + gridmapvar_name="moved", + xco_name="longitude", + xco_units="degrees_east", + yco_name="latitude", + yco_units="degrees_north", + ) + self.check_result(result, cube_no_cs=True) + + def test_mapping__mismatch__rotated_coords_latlon_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_latitude_longitude + # 003 : fc_provides_coordinate_latitude + # 004 : fc_provides_coordinate_longitude + # NOTES: + # no build_coord triggers : requires NO latlon/rotated mapping + # hence no coords at all + result = self.run_testcase( + xco_name="grid_longitude", + xco_units="degrees", + yco_name="grid_latitude", + yco_units="degrees", + ) + self.check_result(result, cube_no_cs=True, cube_no_xycoords=True) + + def test_mapping__mismatch__rotated_coords_nonll_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_albers_equal_area + # 003 : fc_provides_coordinate_latitude + # 004 : fc_provides_coordinate_longitude + # 005 : fc_build_coordinate_latitude_nocs + # 006 : fc_build_coordinate_longitude_nocs + # NOTES: + # this is different from the previous + # build_coord.._nocs triggers : requires NO latlon/rotated mapping + # - which seems odd + inconsistent (with previous) ? + # TODO: should this change ?? + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ALBERS, + xco_name="grid_longitude", + xco_units="degrees", + yco_name="grid_latitude", + yco_units="degrees", + ) + self.check_result(result, cube_no_cs=True) + + def test_mapping__mismatch__rotated_coords_missing_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_latitude + # 003 : fc_provides_coordinate_longitude + # 004 : fc_build_coordinate_latitude_nocs + # 005 : fc_build_coordinate_longitude_nocs + # NOTES: + # as previous, but no grid-mapping (which makes no difference) + warning = "Missing.*grid mapping variable 'grid'" + result = self.run_testcase( + warning=warning, + gridmapvar_name="moved", + xco_name="grid_longitude", + xco_units="degrees", + yco_name="grid_latitude", + yco_units="degrees", + ) + self.check_result(result, cube_no_cs=True) + + def test_mapping__mismatch__nonll_coords_latlon_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_latitude_longitude + # 003 : fc_default_coordinate + # 004 : fc_default_coordinate + # NOTES: + # dim-coords built as "defaults" : dim-coords, but NO standard name + result = self.run_testcase( + xco_name="projection_x", + xco_units="m", + yco_name="projection_y", + yco_units="m", + ) + self.check_result( + result, cube_no_cs=True, xco_stdname=False, yco_stdname=False + ) + + def test_mapping__mismatch__nonll_coords_rotated_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_rotated_latitude_longitude + # 003 : fc_default_coordinate + # 004 : fc_default_coordinate + # NOTES: + # same as previous __mismatch__nonll_ + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ROTATED_LAT_LON, + xco_name="projection_x", + xco_units="m", + yco_name="projection_y", + yco_units="m", + ) + self.check_result( + result, cube_no_cs=True, xco_stdname=False, yco_stdname=False + ) + + def test_mapping__mismatch__nonll_coords_missing_system(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_default_coordinate + # 003 : fc_default_coordinate + # NOTES: + # effectively, just like previous 2 __mismatch__nonll_ + warning = "Missing.*grid mapping variable 'grid'" + result = self.run_testcase( + warning=warning, + gridmapvar_name="moved", + xco_name="projection_x", + xco_units="m", + yco_name="projection_y", + yco_units="m", + ) + self.check_result( + result, cube_no_cs=True, xco_stdname=False, yco_stdname=False + ) + + +class Test__grid_mapping__pyke_rules( + Mixin__grid_mapping__tests, tests.IrisTest +): + # Run grid-mapping tests with Pyke (rules) + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +from unittest import skip + + +@skip +class Test__grid_mapping__nonpyke_actions( + Mixin__grid_mapping__tests, tests.IrisTest +): + # Run grid-mapping tests with non-Pyke (actions) + use_pyke = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +class Mixin__aux_latlons(Mixin__grid_mapping): + # Testcases for translating auxiliary latitude+longitude variables + + def test_aux_lon(self): + # Change the name of xdim, and put xco on the coords list. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_latitude_longitude + # 003 : fc_provides_coordinate_latitude + # 004 : fc_build_auxiliary_coordinate_longitude + # 005 : fc_build_coordinate_latitude + result = self.run_testcase(xco_is_dim=False) + self.check_result(result, xco_is_aux=True, xco_no_cs=True) + + def test_aux_lat(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_latitude_longitude + # 003 : fc_provides_coordinate_longitude + # 004 : fc_build_auxiliary_coordinate_latitude + # 005 : fc_build_coordinate_longitude + result = self.run_testcase(yco_is_dim=False) + self.check_result(result, yco_is_aux=True, yco_no_cs=True) + + def test_aux_lat_and_lon(self): + # When *both* are aux, the grid-mapping is discarded. + # - as in this case there are then no dim-coords to reference it. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_latitude_longitude + # 003 : fc_build_auxiliary_coordinate_latitude + # 004 : fc_build_auxiliary_coordinate_longitude + result = self.run_testcase(xco_is_dim=False, yco_is_dim=False) + self.check_result( + result, xco_is_aux=True, yco_is_aux=True, cube_no_cs=True + ) + + def test_aux_lon_rotated(self): + # Same but with rotated-style lat + lon coords. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_rotated_latitude_longitude + # 003 : fc_provides_coordinate_latitude + # 004 : fc_build_auxiliary_coordinate_longitude_rotated + # 005 : fc_build_coordinate_latitude_rotated + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ROTATED_LAT_LON, + xco_is_dim=False, + ) + self.check_result(result, xco_is_aux=True, xco_no_cs=True) + + def test_aux_lat_rotated(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_grid_mapping_rotated_latitude_longitude + # 003 : fc_provides_coordinate_latitude + # 004 : fc_build_auxiliary_coordinate_longitude_rotated + # 005 : fc_build_coordinate_latitude_rotated + result = self.run_testcase( + mapping_type_name=hh.CF_GRID_MAPPING_ROTATED_LAT_LON, + xco_is_dim=False, + ) + self.check_result(result, xco_is_aux=True, xco_no_cs=True) + + +class Test__aux_latlons__pyke_rules(Mixin__aux_latlons, tests.IrisTest): + # Run aux-latlons tests with Pyke (rules) + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__hybrid_formulae.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__hybrid_formulae.py new file mode 100644 index 0000000000..c1a325925f --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__hybrid_formulae.py @@ -0,0 +1,225 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for the engine.activate() call within the +`iris.fileformats.netcdf._load_cube` function. + +Test rules activation relating to hybrid vertical coordinates. + +""" +import iris.tests as tests + +import iris.fileformats._nc_load_rules.helpers as hh +from iris.tests.unit.fileformats.netcdf.load_cube.load_cube__activate import ( + Mixin__nc_load_actions, +) + + +class Mixin__formulae_tests(Mixin__nc_load_actions): + def _make_testcase_cdl(self, formula_root_name=None, term_names=None): + """Construct a testcase CDL for data with hybrid vertical coords.""" + if formula_root_name is None: + formula_root_name = "atmosphere_hybrid_height_coordinate" + if term_names is None: + term_names = hh.CF_COORD_VERTICAL.get(formula_root_name) + if term_names is None: + # unsupported type : just make something up + term_names = ["term1"] + + terms_string = "" + phenom_coord_names = ["vert"] # always include the root variable + formula_term_strings = [] + for term_name in term_names: + term_varname = "v_" + term_name + phenom_coord_names.append(term_varname) + formula_term_strings.append(f"{term_name}: {term_varname}") + terms_string += f""" + double {term_varname}(h) ; + {term_varname}:long_name = "{term_name}_long_name" ; + {term_varname}:units = "m" ; +""" + + # remove the extra initial space from the formula terms string + phenom_coords_string = " ".join(phenom_coord_names) + formula_terms_string = " ".join(formula_term_strings) + # Create the main result string. + cdl_str = f""" +netcdf test {{ +dimensions: + h = 2 ; +variables: + double phenom(h) ; + phenom:standard_name = "air_temperature" ; + phenom:units = "K" ; + phenom:coordinates = "{phenom_coords_string}" ; + double vert(h) ; + vert:standard_name = "{formula_root_name}" ; + vert:long_name = "hybrid_vertical" ; + vert:units = "m" ; + vert:formula_terms = "{formula_terms_string}" ; +{terms_string} +}} +""" + return cdl_str + + def check_result(self, cube, factory_type="_auto", formula_terms="_auto"): + """Check the result of a cube load with a hybrid vertical coord.""" + if factory_type == "_auto": + # replace with our 'default', which is hybrid-height. + # N.B. 'None' is different: it means expect *no* factory. + factory_type = "atmosphere_hybrid_height_coordinate" + self.assertEqual(cube._formula_type_name, factory_type) + + if formula_terms == "_auto": + # Set default terms-expected, according to the expected factory + # type. + if factory_type is None: + # If no factory, expect no identified terms. + formula_terms = [] + else: + # Expect the correct ones defined for the factory type. + formula_terms = hh.CF_COORD_VERTICAL[factory_type] + + # Compare the formula_terms list with the 'expected' ones. + # N.B. first make the 'expected' list lower case, as the lists in + # hh.CF_COORD_VERTICAL include uppercase, but rules outputs don't. + formula_terms = [term.lower() for term in formula_terms] + + # N.B. the terms dictionary can be missing, if there were none + actual_terms = cube._formula_terms_byname or {} + self.assertEqual(sorted(formula_terms), sorted(actual_terms.keys())) + + # Check that there is an aux-coord of the expected name for each term + for var_name in actual_terms.values(): + coords = cube.coords(var_name=var_name, dim_coords=False) + self.assertEqual(len(coords), 1) + + # + # Actual testcase routines + # + + def test_basic_hybridheight(self): + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate + # 003 : fc_build_auxiliary_coordinate + # 004 : fc_build_auxiliary_coordinate + # 005 : fc_build_auxiliary_coordinate + # 006 : fc_build_auxiliary_coordinate + # 007 : fc_formula_type_atmosphere_hybrid_sigma_pressure_coordinate + # 008 : fc_formula_terms + # 009 : fc_formula_terms + # 010 : fc_formula_terms + # 011 : fc_formula_terms + result = self.run_testcase() + self.check_result(result) + + def test_missing_term(self): + # Check behaviour when a term is missing. + # For the test, omit "orography", which is common in practice. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate + # 003 : fc_build_auxiliary_coordinate + # 004 : fc_build_auxiliary_coordinate + # 005 : fc_formula_type_atmosphere_hybrid_height_coordinate + # 006 : fc_formula_terms + # 007 : fc_formula_terms + result = self.run_testcase( + term_names=["a", "b"] # missing the 'orog' term + ) + self.check_result(result, formula_terms=["a", "b"]) + + def test_no_terms(self): + # Check behaviour when *all* terms are missing. + # N.B. for any _actual_ type, this is probably invalid and would fail? + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate + result = self.run_testcase( + formula_root_name="atmosphere_hybrid_height_coordinate", + term_names=[], + ) + # This does *not* trigger + # 'fc_formula_type_atmosphere_hybrid_height_coordinate' + # This is because, within the 'assert_case_specific_facts' routine, + # formula_roots are only recognised by scanning the identified + # formula_terms. + self.check_result(result, factory_type=None) + + def test_unrecognised_verticaltype(self): + # Set the root variable name to something NOT a recognised hybrid type. + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate + # 003 : fc_build_auxiliary_coordinate + # 004 : fc_build_auxiliary_coordinate + # 005 : fc_formula_terms + # 006 : fc_formula_terms + result = self.run_testcase( + formula_root_name="unknown", term_names=["a", "b"] + ) + # Check that it picks up the terms, but *not* the factory root coord, + # which is simply discarded. + self.check_result(result, factory_type=None, formula_terms=["a", "b"]) + + +# Add in tests methods to exercise each (supported) vertical coordinate type +# individually. +# NOTE: hh.CF_COORD_VERTICAL lists all the valid types, but we don't yet +# support all of them. +_SUPPORTED_FORMULA_TYPES = ( + # NOTE: omit "atmosphere_hybrid_height_coordinate" : our basic testcase + "atmosphere_hybrid_sigma_pressure_coordinate", + "ocean_sigma_z_coordinate", + "ocean_sigma_coordinate", + "ocean_s_coordinate", + "ocean_s_coordinate_g1", + "ocean_s_coordinate_g2", +) +for hybrid_type in _SUPPORTED_FORMULA_TYPES: + + def construct_inner_func(hybrid_type): + term_names = hh.CF_COORD_VERTICAL[hybrid_type] + + def inner(self): + result = self.run_testcase( + formula_root_name=hybrid_type, term_names=term_names + ) + self.check_result( + result, factory_type=hybrid_type, formula_terms=term_names + ) + + return inner + + # Note: use an intermediate function to generate each test method, simply to + # generate a new local variable for 'hybrid_type' on each iteration. + # Otherwise all the test methods will refer to the *same* 'hybrid_type' + # variable, i.e. the loop variable, which does not work ! + method_name = f"test_{hybrid_type}_coord" + setattr( + Mixin__formulae_tests, method_name, construct_inner_func(hybrid_type) + ) + + +class Test__formulae__withpyke(Mixin__formulae_tests, tests.IrisTest): + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__miscellaneous.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__miscellaneous.py new file mode 100644 index 0000000000..d41b19e108 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__miscellaneous.py @@ -0,0 +1,216 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for the engine.activate() call within the +`iris.fileformats.netcdf._load_cube` function. + +Tests for rules activation relating to some isolated aspects : + * UKMO um-specific metadata + * label coordinates + * cell measures + * ancillary variables + +""" +import iris.tests as tests + +from iris.coords import AuxCoord, CellMeasure, AncillaryVariable +from iris.fileformats.pp import STASH + +from iris.tests.unit.fileformats.netcdf.load_cube.load_cube__activate import ( + Mixin__nc_load_actions, +) + + +class Mixin__ukmo_attributes(Mixin__nc_load_actions): + # Tests for handling of the special UM-specific data-var attributes. + def _make_testcase_cdl(self, **add_attrs): + phenom_attrs_string = "" + for key, value in add_attrs.items(): + phenom_attrs_string += f""" + phenom:{key} = "{value}" ; +""" + + cdl_string = f""" +netcdf test {{ + dimensions: + xdim = 2 ; + variables: + double phenom(xdim) ; + phenom:standard_name = "air_temperature" ; + phenom:units = "K" ; +{phenom_attrs_string} +}} +""" + return cdl_string + + def check_result(self, cube, stashcode=None, processflags=None): + cube_stashattr = cube.attributes.get("STASH") + cube_processflags = cube.attributes.get("ukmo__process_flags") + + if stashcode is not None: + self.assertIsInstance(cube_stashattr, STASH) + self.assertEqual(str(stashcode), str(cube_stashattr)) + else: + self.assertIsNone(cube_stashattr) + + if processflags is not None: + self.assertIsInstance(cube_processflags, tuple) + self.assertEqual(set(cube_processflags), set(processflags)) + else: + self.assertIsNone(cube_processflags) + + # + # Testcase routines + # + stashcode = "m01s02i034" # Just one valid STASH msi string for testing + + def test_stash(self): + cube = self.run_testcase(um_stash_source=self.stashcode) + self.check_result(cube, stashcode=self.stashcode) + + def test_stash_altname(self): + cube = self.run_testcase(ukmo__um_stash_source=self.stashcode) + self.check_result(cube, stashcode=self.stashcode) + + def test_stash_empty(self): + msg = "Expected STASH code MSI string" + with self.assertRaisesRegex(ValueError, msg): + self.run_testcase(ukmo__um_stash_source="") + + def test_stash_invalid(self): + msg = "Expected STASH code MSI string" + with self.assertRaisesRegex(ValueError, msg): + self.run_testcase(ukmo__um_stash_source="XXX") + + def test_processflags_single(self): + cube = self.run_testcase(ukmo__process_flags="this") + self.check_result(cube, processflags=["this"]) + + def test_processflags_multi_with_underscores(self): + flags_testinput = "this that_1 the_other_one x" + flags_expectresult = ["this", "that 1", "the other one", "x"] + cube = self.run_testcase(ukmo__process_flags=flags_testinput) + self.check_result(cube, processflags=flags_expectresult) + + def test_processflags_empty(self): + cube = self.run_testcase(ukmo__process_flags="") + expected_result = [""] # May seem odd, but that's what it does. + self.check_result(cube, processflags=expected_result) + + +class Test__ukmo_attributes__withpyke(Mixin__ukmo_attributes, tests.IrisTest): + use_pyke = True + + +class Mixin__labels_cellmeasures_ancils(Mixin__nc_load_actions): + # Tests for some simple rules that translate facts directly into cube data, + # with no alternative actions, complications or failure modes to test. + def _make_testcase_cdl( + self, + include_label=False, + include_cellmeasure=False, + include_ancil=False, + ): + + phenom_extra_attrs_string = "" + extra_vars_string = "" + + if include_label: + phenom_extra_attrs_string += """ + phenom:coordinates = "v_label" ; +""" + extra_vars_string += """ + char v_label(xdim, strdim) ; + v_label:long_name = "string data" ; +""" + + if include_cellmeasure: + # One simple case : a valid link + a variable definition. + phenom_extra_attrs_string += """ + phenom:cell_measures = "area: v_cellm" ; +""" + extra_vars_string += """ + double v_cellm(xdim) ; + v_cellm:long_name = "cell areas" ; +""" + + if include_ancil: + # One simple case : a valid link + a variable definition. + phenom_extra_attrs_string += """ + phenom:ancillary_variables = "v_ancil" ; +""" + extra_vars_string += """ + double v_ancil(xdim) ; + v_ancil:long_name = "ancillary values" ; +""" + cdl_string = f""" + netcdf test {{ + dimensions: + xdim = 2 ; + strdim = 5 ; + variables: + double phenom(xdim) ; + phenom:standard_name = "air_temperature" ; + phenom:units = "K" ; +{phenom_extra_attrs_string} +{extra_vars_string} + }} + """ + return cdl_string + + def check_result( + self, + cube, + expect_label=False, + expect_cellmeasure=False, + expect_ancil=False, + ): + label_coords = cube.coords(var_name="v_label") + if expect_label: + self.assertEqual(len(label_coords), 1) + (coord,) = label_coords + self.assertIsInstance(coord, AuxCoord) + self.assertEqual(coord.dtype.kind, "U") + else: + self.assertEqual(len(label_coords), 0) + + cell_measures = cube.cell_measures() + if expect_cellmeasure: + self.assertEqual(len(cell_measures), 1) + (cellm,) = cell_measures + self.assertIsInstance(cellm, CellMeasure) + else: + self.assertEqual(len(cell_measures), 0) + + ancils = cube.ancillary_variables() + if expect_ancil: + self.assertEqual(len(ancils), 1) + (ancil,) = ancils + self.assertIsInstance(ancil, AncillaryVariable) + else: + self.assertEqual(len(ancils), 0) + + def test_label(self): + cube = self.run_testcase(include_label=True) + self.check_result(cube, expect_label=True) + + def test_ancil(self): + cube = self.run_testcase(include_ancil=True) + self.check_result(cube, expect_ancil=True) + + def test_cellmeasure(self): + cube = self.run_testcase(include_cellmeasure=True) + self.check_result(cube, expect_cellmeasure=True) + + +class Test__labels_cellmeasures_ancils__withpyke( + Mixin__labels_cellmeasures_ancils, tests.IrisTest +): + use_pyke = True + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__time_coords.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__time_coords.py new file mode 100644 index 0000000000..fa010f446d --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/load_cube/load_cube__activate/test__time_coords.py @@ -0,0 +1,476 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Unit tests for the engine.activate() call within the +`iris.fileformats.netcdf._load_cube` function. + +Tests for rules activation relating to 'time' and 'time_period' coords. + +""" +import iris.tests as tests + +from iris.coords import AuxCoord, DimCoord + + +from iris.tests.unit.fileformats.netcdf.load_cube.load_cube__activate import ( + Mixin__nc_load_actions, +) + + +class Opts(dict): + # A dict-like thing which provides '.' access in place of indexing. + def __init__(self, **kwargs): + # Init like a dict + super().__init__(**kwargs) + # Alias contents "self['key']", as properties "self.key" + # See: https://stackoverflow.com/a/14620633/2615050 + self.__dict__ = self + + +# Per-coord options settings for testcase definitions. +_COORD_OPTIONS_TEMPLATE = { + "which": "", # set to "something" + "stdname": "_auto_which", # default = time / time_period + "varname": "_as_which", # default = time / period + "dimname": "_as_which", + "in_phenomvar_dims": True, + "in_phenomvar_coords": False, # set for an aux-coord + "values_all_zero": False, # set to block CFDimensionVariable identity + "units": "_auto_which", # specific to time/period +} + + +class Mixin__timecoords__common(Mixin__nc_load_actions): + def _make_testcase_cdl( + self, + phenom_dims="_auto", # =get from time+period opts + phenom_coords="_auto", # =get from time+period opts + time_opts=None, + period_opts=None, + timedim_name="time", + perioddim_name="period", + ): + opt_t = None + opt_p = None + if time_opts is not None: + # Replace 'True' with an options dict for 'time' options + opt_t = Opts(**_COORD_OPTIONS_TEMPLATE) + opt_t.update(which="time", **time_opts) + if period_opts is not None: + # Replace 'True' with an options dict for 'period' options + opt_p = Opts(**_COORD_OPTIONS_TEMPLATE) + opt_p.update(which="period", **period_opts) + + # Define the 'standard' dimensions which we will create + # NB we don't necessarily *use* either of these + dims_and_lens = {timedim_name: 2, perioddim_name: 3} + dims_string = "\n".join( + [ + f" {name} = {length} ;" + for name, length in dims_and_lens.items() + ] + ) + + phenom_auto_dims = [] + phenom_auto_coords = [] + coord_variables_string = "" + data_string = "" + for opt in (opt_t, opt_p): + # Handle computed defaults and common info for both coord options. + if opt: + if opt.which not in ("time", "period"): + raise ValueError(f"unrecognised opt.which={opt.which}") + + # Do computed defaults. + if opt.stdname == "_auto_which": + if opt.which == "time": + opt.stdname = "time" + else: + assert opt.which == "period" + opt.stdname = "forecast_period" + if opt.varname == "_as_which": + opt.varname = opt.which + if opt.dimname == "_as_which": + opt.dimname = opt.which + if opt.units == "_auto_which": + if opt.which == "time": + opt.units = "hours since 2000-01-01" + else: + assert opt.which == "period" + opt.units = "hours" + + # Build 'auto' lists of phenom dims and (aux) coordinates. + if opt.in_phenomvar_dims: + phenom_auto_dims.append(opt.dimname) + if opt.in_phenomvar_coords: + phenom_auto_coords.append(opt.varname) + + # Add a definition of the coord variable. + coord_variables_string += f""" + double {opt.varname}({opt.dimname}) ; + {opt.varname}:standard_name = "{opt.stdname}" ; + {opt.varname}:units = "{opt.units}" ; +""" + # NOTE: we don't bother with an 'axis' property. + # We can probe the behaviour we need without that, because we + # are *not* testing the cf.py categorisation code, or the + # helper "build_xxx" routines. + + # Define coord-var data values (so it can be a dimension). + varname = opt.varname + if opt.values_all_zero: + # Use 'values_all_zero' to prevent a dim-var from + # identifying as a CFDimensionCoordinate (as it is + # non-monotonic). + dim_vals = [0.0] * dims_and_lens[opt.dimname] + else: + # "otherwise", assign an ascending sequence. + dim_vals = range(dims_and_lens[opt.dimname]) + dimvals_string = ", ".join(f"{val:0.1f}" for val in dim_vals) + data_string += f"\n {varname} = {dimvals_string} ;" + + if phenom_dims == "_auto": + phenom_dims = phenom_auto_dims + if not phenom_dims: + phenom_dims_string = "" + else: + phenom_dims_string = ", ".join(phenom_dims) + + if phenom_coords == "_auto": + phenom_coords = phenom_auto_coords + if not phenom_coords: + phenom_coords_string = "" + else: + phenom_coords_string = " ".join(phenom_coords) + phenom_coords_string = ( + " " + f'phenom:coordinates = "{phenom_coords_string}" ; ' + ) + + # Create a testcase with time dims + coords. + cdl_string = f""" +netcdf test {{ + dimensions: +{dims_string} + variables: + double phenom({phenom_dims_string}) ; + phenom:standard_name = "air_temperature" ; + phenom:units = "K" ; +{phenom_coords_string} + +{coord_variables_string} + data: +{data_string} +}} +""" + return cdl_string + + def check_result(self, cube, time_is="dim", period_is="missing"): + """ + Check presence of expected dim/aux-coords in the result cube. + + Both of 'time_is' and 'period_is' can take values 'dim', 'aux' or + 'missing'. + + """ + options = ("dim", "aux", "missing") + msg = f'Invalid "{{name}}" = {{opt}} : Not one of {options!r}.' + if time_is not in options: + raise ValueError(msg.format(name="time_is", opt=time_is)) + if period_is not in options: + raise ValueError(msg.format(name="period_is", opt=period_is)) + + # Get the facts we want to check + time_name = "time" + period_name = "forecast_period" + time_dimcos = cube.coords(time_name, dim_coords=True) + time_auxcos = cube.coords(time_name, dim_coords=False) + period_dimcos = cube.coords(period_name, dim_coords=True) + period_auxcos = cube.coords(period_name, dim_coords=False) + + if time_is == "dim": + self.assertEqual(len(time_dimcos), 1) + self.assertEqual(len(time_auxcos), 0) + elif time_is == "aux": + self.assertEqual(len(time_dimcos), 0) + self.assertEqual(len(time_auxcos), 1) + else: + self.assertEqual(len(time_dimcos), 0) + self.assertEqual(len(time_auxcos), 0) + + if period_is == "dim": + self.assertEqual(len(period_dimcos), 1) + self.assertEqual(len(period_auxcos), 0) + elif period_is == "aux": + self.assertEqual(len(period_dimcos), 0) + self.assertEqual(len(period_auxcos), 1) + else: + self.assertEqual(len(period_dimcos), 0) + self.assertEqual(len(period_auxcos), 0) + + # Also check expected built Coord types. + if time_is == "dim": + self.assertIsInstance(time_dimcos[0], DimCoord) + elif time_is == "aux": + self.assertIsInstance(time_auxcos[0], AuxCoord) + + if period_is == "dim": + self.assertIsInstance(period_dimcos[0], DimCoord) + elif period_is == "aux": + self.assertIsInstance(period_auxcos[0], AuxCoord) + + +class Mixin__singlecoord__tests(Mixin__timecoords__common): + # Coordinate tests to be run for both 'time' and 'period' coordinate vars. + # Set (in inheritors) to select time/period testing. + which = None + + def run_testcase(self, coord_dim_name=None, **opts): + """ + Specialise 'run_testcase' for single-coord 'time' or 'period' testing. + """ + which = self.which + assert which in ("time", "period") + + # Separate the 'Opt' keywords from "others" : others are passed + # directly to the parent routine, whereas 'Opt' ones are passed to + # 'time_opts' / 'period_opts' keys accordingly. + general_opts = {} + for key, value in list(opts.items()): + if key not in _COORD_OPTIONS_TEMPLATE.keys(): + del opts[key] + general_opts[key] = value + + if coord_dim_name is not None: + # Translate this into one of timedim_name/perioddim_name + general_opts[f"{which}dim_name"] = coord_dim_name + + period_opts = None + time_opts = None + if which == "time": + time_opts = opts + else: + period_opts = opts + + result = super().run_testcase( + time_opts=time_opts, period_opts=period_opts, **general_opts + ) + + return result + + def check_result(self, cube, coord_is="dim"): + """ + Specialise 'check_result' for single-coord 'time' or 'period' testing. + """ + # Pass generic 'coord_is' option to parent as time/period options. + which = self.which + assert which in ("time", "period") + + if which == "time": + time_is = coord_is + period_is = "missing" + else: + period_is = coord_is + time_is = "missing" + + super().check_result(cube, time_is=time_is, period_is=period_is) + + # + # Generic single-coordinate testcases. + # ( these are repeated for both 'time' and 'time_period' ) + # + + def test_dimension(self): + # Coord is a normal dimension --> dimcoord + # + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_time + # 003 : fc_build_coordinate_time + result = self.run_testcase() + self.check_result(result, "dim") + + def test_dimension_in_phenom_coords(self): + # Dimension coord also present in phenom:coords. + # Strictly wrong but a common error in datafiles : must tolerate. + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_time + # 003 : fc_build_coordinate_time + result = self.run_testcase(in_phenomvar_coords=True) + self.check_result(result, "dim") + + def test_dim_nonmonotonic(self): + # Coord has all-zero values, which prevents it being a dimcoord. + # The rule has a special way of treating it as an aux coord + # -- even though it doesn't appear in the phenom coords. + # ( Done by the build_coord routine, so not really a rules issue). + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_time + # 003 : fc_build_coordinate_time + msg = "Failed to create.* dimension coordinate" + result = self.run_testcase(values_all_zero=True, warning=msg) + self.check_result(result, "aux") + + def test_dim_fails_typeident(self): + # The coord variable is identified as a CFDimensionCoordinate by cf.py, + # but having the wrong units causes it to fail the 'is_time' or + # 'is_period' test, so the 'provides_coord' rule fails to trigger. + # So it is built as a 'miscellaneous' dim-coord. + # N.B. this makes *no* practical difference, because a 'misc' dim + # coord is still a dim coord (albeit with bad units). + # ( N.B.#2 Not quite the same for lat/lon coords, where coord-specific + # 'build' rules always use a fixed standard-name ). + # Rules Triggered: + # #001 : fc_default + # #002 : fc_default_coordinate + result = self.run_testcase(units="1") + self.check_result(result, "dim") + + def test_aux(self): + # time/period is installed as an auxiliary coord. + # For this, rename both DIMENSIONS, so that the generated coords are + # not actually CF coordinates. + # For a valid case, we must *also* have a ref in phenom:coordinates + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate_time + result = self.run_testcase( + coord_dim_name="dim_renamed", + dimname="dim_renamed", + in_phenomvar_coords=True, + ) + self.check_result(result, "aux") + + def test_aux_not_in_phenom_coords(self): + # time/period is installed as an auxiliary coord, + # but we DIDN'T list it in phenom:coords -- otherwise as previous. + # Should have no result at all. + # Rules Triggered: + # 001 : fc_default + result = self.run_testcase( + coord_dim_name="dim_renamed", + dimname="dim_renamed", + in_phenomvar_coords=False, + ) # "should" be True for an aux-coord + self.check_result(result, "missing") + + def test_aux_fails_typeident(self): + # The coord variable is identified as a CFAuxiliaryCoordinate by cf.py, + # but having the wrong units causes it to fail the 'is_time' or + # 'is_period' test, so the 'provides_coord' rule fails to trigger. + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate + # Again, though it builds as a 'miscellaneous' rather than a recognised + # specific coord type, it makes no practical difference. + result = self.run_testcase( + coord_dim_name="dim_renamed", + dimname="dim_renamed", + in_phenomvar_coords=True, + units="1", + ) + self.check_result(result, "aux") + + def test_aux_no_coordsref(self): + # The coord variable is identified as a CFAuxiliaryCoordinate by cf.py, + # but having the wrong units causes it to fail the 'is_time' or + # 'is_period' test. + # Rules Triggered: + # 001 : fc_default + # 002 : fc_build_auxiliary_coordinate + # Again, though it builds as a 'miscellaneous' rather than a reocgnised + # specific coord type, it makes no practical difference. + result = self.run_testcase( + coord_dim_name="dim_renamed", + dimname="dim_renamed", + in_phenomvar_coords=True, + units="1", + ) + self.check_result(result, "aux") + + +class Test__time__withpyke(Mixin__singlecoord__tests, tests.IrisTest): + # Run 'time' coord tests + which = "time" + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +class Test__period__withpyke(Mixin__singlecoord__tests, tests.IrisTest): + # Run 'time_period' coord tests + which = "period" + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +class Mixin__dualcoord__tests(Mixin__timecoords__common): + # Coordinate tests for a combination of 'time' and 'time_period'. + # Not strictly necessary, as handling is independent, but a handy check + # on typical usage. + def test_time_and_period(self): + # Test case with both 'time' and 'period', with separate dims. + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_time + # 003 : fc_provides_coordinate_time_period + # 004 : fc_build_coordinate_time + # 005 : fc_build_coordinate_time_period + result = self.run_testcase(time_opts={}, period_opts={}) + self.check_result(result, time_is="dim", period_is="dim") + + def test_time_dim_period_aux(self): + # Test case with both 'time' and 'period' sharing a dim. + # Rules Triggered: + # 001 : fc_default + # Rules Triggered: + # 001 : fc_default + # 002 : fc_provides_coordinate_time + # 003 : fc_build_auxiliary_coordinate_time_period + # 004 : fc_build_coordinate_time + result = self.run_testcase( + time_opts={}, + period_opts=dict( + dimname="time", + in_phenomvar_dims=False, + in_phenomvar_coords=True, + ), + ) + self.check_result(result, time_is="dim", period_is="aux") + + +class Test__dualcoord_tests__withpyke(Mixin__dualcoord__tests, tests.IrisTest): + use_pyke = True + debug = False + + @classmethod + def setUpClass(cls): + super().setUpClass() + + @classmethod + def tearDownClass(cls): + super().tearDownClass() + + +if __name__ == "__main__": + tests.main() diff --git a/lib/iris/tests/unit/fileformats/netcdf/test__load_cube.py b/lib/iris/tests/unit/fileformats/netcdf/load_cube/test__load_cube.py similarity index 100% rename from lib/iris/tests/unit/fileformats/netcdf/test__load_cube.py rename to lib/iris/tests/unit/fileformats/netcdf/load_cube/test__load_cube.py