diff --git a/esmvalcore/cmor/_fixes/fix.py b/esmvalcore/cmor/_fixes/fix.py index 973ac57d0b..4d3e297e3a 100644 --- a/esmvalcore/cmor/_fixes/fix.py +++ b/esmvalcore/cmor/_fixes/fix.py @@ -27,7 +27,7 @@ ) from esmvalcore.cmor.fixes import get_time_bounds from esmvalcore.cmor.table import get_var_info -from esmvalcore.iris_helpers import has_unstructured_grid +from esmvalcore.iris_helpers import has_unstructured_grid, safe_convert_units if TYPE_CHECKING: from esmvalcore.cmor.table import CoordinateInfo, VariableInfo @@ -455,7 +455,7 @@ def _fix_units(self, cube: Cube) -> Cube: if str(cube.units) != units: old_units = cube.units try: - cube.convert_units(units) + safe_convert_units(cube, units) except (ValueError, UnitConversionError): self._warning_msg( cube, diff --git a/esmvalcore/cmor/_fixes/native6/era5.py b/esmvalcore/cmor/_fixes/native6/era5.py index aa4beae984..79fa2856fd 100644 --- a/esmvalcore/cmor/_fixes/native6/era5.py +++ b/esmvalcore/cmor/_fixes/native6/era5.py @@ -6,7 +6,7 @@ import iris import numpy as np -from esmvalcore.iris_helpers import date2num +from esmvalcore.iris_helpers import date2num, safe_convert_units from ...table import CMOR_TABLES from ..fix import Fix @@ -414,10 +414,6 @@ def _fix_monthly_time_coord(cube): coord.points = 0.5 * (start + end) coord.bounds = np.column_stack([start, end]) - def _fix_units(self, cube): - """Fix units.""" - cube.convert_units(self.vardef.units) - def fix_metadata(self, cubes): """Fix metadata.""" fixed_cubes = iris.cube.CubeList() @@ -428,7 +424,7 @@ def fix_metadata(self, cubes): cube.long_name = self.vardef.long_name cube = self._fix_coordinates(cube) - self._fix_units(cube) + cube = safe_convert_units(cube, self.vardef.units) cube.data = cube.core_data().astype("float32") year = datetime.datetime.now().year diff --git a/esmvalcore/cmor/_fixes/native_datasets.py b/esmvalcore/cmor/_fixes/native_datasets.py index 6b48d92f54..da6f470816 100644 --- a/esmvalcore/cmor/_fixes/native_datasets.py +++ b/esmvalcore/cmor/_fixes/native_datasets.py @@ -5,6 +5,8 @@ from iris import NameConstraint +from esmvalcore.iris_helpers import safe_convert_units + from ..fix import Fix from .shared import ( add_scalar_height_coord, @@ -77,7 +79,7 @@ def fix_var_metadata(self, cube): f"Failed to fix invalid units '{invalid_units}' for " f"variable '{self.vardef.short_name}'" ) from exc - cube.convert_units(self.vardef.units) + safe_convert_units(cube, self.vardef.units) # Fix attributes if self.vardef.positive != "": diff --git a/esmvalcore/iris_helpers.py b/esmvalcore/iris_helpers.py index 4162233eec..8d1c676682 100644 --- a/esmvalcore/iris_helpers.py +++ b/esmvalcore/iris_helpers.py @@ -9,6 +9,7 @@ import iris.cube import iris.util import numpy as np +from cf_units import Unit from iris.coords import Coord from iris.cube import Cube from iris.exceptions import CoordinateMultiDimError, CoordinateNotFoundError @@ -92,8 +93,8 @@ def date2num(date, unit, dtype=np.float64): This is a custom version of :meth:`cf_units.Unit.date2num` that guarantees the correct dtype for the return value. - Arguments - --------- + Parameters + ---------- date : :class:`datetime.datetime` or :class:`cftime.datetime` unit : :class:`cf_units.Unit` dtype : a numpy dtype @@ -358,3 +359,121 @@ def has_unstructured_grid(cube: Cube) -> bool: if cube.coord_dims(lat) != cube.coord_dims(lon): return False return True + + +# List containing special cases for unit conversion. Each list item is another +# list. Each of these sublists defines one special conversion. Each element in +# the sublists is a tuple (standard_name, units). Note: All units for a single +# special case need to be "physically identical", e.g., 1 kg m-2 s-1 "equals" 1 +# mm s-1 for precipitation +_SPECIAL_UNIT_CONVERSIONS = [ + [ + ("precipitation_flux", "kg m-2 s-1"), + ("lwe_precipitation_rate", "mm s-1"), + ], + [ + ("equivalent_thickness_at_stp_of_atmosphere_ozone_content", "m"), + ("equivalent_thickness_at_stp_of_atmosphere_ozone_content", "1e5 DU"), + ], +] + + +def _try_special_unit_conversions(cube: Cube, units: str | Unit) -> bool: + """Try special unit conversion (in-place). + + Parameters + ---------- + cube: + Input cube (modified in place). + units: + New units + + Returns + ------- + bool + ``True`` if special unit conversion was successful, ``False`` if not. + + """ + for special_case in _SPECIAL_UNIT_CONVERSIONS: + for std_name, special_units in special_case: + # Special unit conversion only works if all of the following + # criteria are met: + # - the cube's standard_name is one of the supported + # standard_names + # - the cube's units are convertible to the ones defined for + # that given standard_name + # - the desired target units are convertible to the units of + # one of the other standard_names in that special case + + # Step 1: find suitable source name and units + if cube.standard_name == std_name and cube.units.is_convertible( + special_units + ): + for target_std_name, target_units in special_case: + if target_units == special_units: + continue + + # Step 2: find suitable target name and units + if Unit(units).is_convertible(target_units): + cube.standard_name = target_std_name + + # In order to avoid two calls to cube.convert_units, + # determine the conversion factor between the cube's + # units and the source units first and simply add this + # factor to the target units (remember that the source + # units and the target units should be "physically + # identical"). + factor = cube.units.convert(1.0, special_units) + cube.units = f"{factor} {target_units}" + cube.convert_units(units) + return True + + # If no special case has been detected, return False + return False + + +def safe_convert_units(cube: Cube, units: str | Unit) -> Cube: + """Safe unit conversion (change of `standard_name` not allowed; in-place). + + This is a safe version of :func:`esmvalcore.preprocessor.convert_units` + that will raise an error if the input cube's + :attr:`~iris.cube.Cube.standard_name` has been changed. + + Parameters + ---------- + cube: + Input cube (modified in place). + units: + New units. + + Returns + ------- + iris.cube.Cube + Converted cube. Just returned for convenience; input cube is modified + in place. + + Raises + ------ + iris.exceptions.UnitConversionError + Old units are unknown. + ValueError + Old units are not convertible to new units or unit conversion required + change of `standard_name`. + + """ + old_units = cube.units + old_standard_name = cube.standard_name + + try: + cube.convert_units(units) + except ValueError: + if not _try_special_unit_conversions(cube, units): + raise + + if cube.standard_name != old_standard_name: + raise ValueError( + f"Cannot safely convert units from '{old_units}' to '{units}'; " + f"standard_name changed from '{old_standard_name}' to " + f"'{cube.standard_name}'" + ) + return cube diff --git a/esmvalcore/preprocessor/_units.py b/esmvalcore/preprocessor/_units.py index 23ebe23cc2..0ef2c133b6 100644 --- a/esmvalcore/preprocessor/_units.py +++ b/esmvalcore/preprocessor/_units.py @@ -11,77 +11,23 @@ import iris import numpy as np from cf_units import Unit +from iris.coords import AuxCoord, DimCoord +from iris.cube import Cube + +from esmvalcore.iris_helpers import _try_special_unit_conversions logger = logging.getLogger(__name__) -# List containing special cases for convert_units. Each list item is another -# list. Each of these sublists defines one special conversion. Each element in -# the sublists is a tuple (standard_name, units). Note: All units for a single -# special case need to be "physically identical", e.g., 1 kg m-2 s-1 "equals" 1 -# mm s-1 for precipitation -SPECIAL_CASES = [ - [ - ("precipitation_flux", "kg m-2 s-1"), - ("lwe_precipitation_rate", "mm s-1"), - ], - [ - ("equivalent_thickness_at_stp_of_atmosphere_ozone_content", "m"), - ("equivalent_thickness_at_stp_of_atmosphere_ozone_content", "1e5 DU"), - ], -] - - -def _try_special_conversions(cube, units): - """Try special conversion.""" - for special_case in SPECIAL_CASES: - for std_name, special_units in special_case: - # Special unit conversion only works if all of the following - # criteria are met: - # - the cube's standard_name is one of the supported - # standard_names - # - the cube's units are convertible to the ones defined for - # that given standard_name - # - the desired target units are convertible to the units of - # one of the other standard_names in that special case - - # Step 1: find suitable source name and units - if cube.standard_name == std_name and cube.units.is_convertible( - special_units - ): - for target_std_name, target_units in special_case: - if target_units == special_units: - continue - - # Step 2: find suitable target name and units - if Unit(units).is_convertible(target_units): - cube.standard_name = target_std_name - - # In order to avoid two calls to cube.convert_units, - # determine the conversion factor between the cube's - # units and the source units first and simply add this - # factor to the target units (remember that the source - # units and the target units should be "physically - # identical"). - factor = cube.units.convert(1.0, special_units) - cube.units = f"{factor} {target_units}" - cube.convert_units(units) - return True - - # If no special case has been detected, return False - return False - - -def convert_units(cube, units): - """Convert the units of a cube to new ones. - - This converts units of a cube. +def convert_units(cube: Cube, units: str | Unit) -> Cube: + """Convert the units of a cube to new ones (in-place). Note ---- Allows special unit conversions which transforms one quantity to another - (physically related) quantity. These quantities are identified via their - ``standard_name`` and their ``units`` (units convertible to the ones + (physically related) quantity, which may also change the input cube's + :attr:`~iris.cube.Cube.standard_name`. These quantities are identified via + their ``standard_name`` and their ``units`` (units convertible to the ones defined are also supported). For example, this enables conversions between precipitation fluxes measured in ``kg m-2 s-1`` and precipitation rates measured in ``mm day-1`` (and vice versa). @@ -103,32 +49,40 @@ def convert_units(cube, units): Note that for precipitation variables, a water density of ``1000 kg m-3`` is assumed. - Arguments - --------- - cube: iris.cube.Cube - Input cube. - units: str - New units in udunits form. + Parameters + ---------- + cube: + Input cube (modified in place). + units: + New units. Returns ------- iris.cube.Cube - converted cube. + Converted cube. Just returned for convenience; input cube is modified + in place. + + Raises + ------ + iris.exceptions.UnitConversionError + Old units are unknown. + ValueError + Old units are not convertible to new units. """ try: cube.convert_units(units) except ValueError: - if not _try_special_conversions(cube, units): + if not _try_special_unit_conversions(cube, units): raise return cube def accumulate_coordinate( - cube: iris.cube.Cube, - coordinate: str | iris.coords.DimCoord | iris.coords.AuxCoord, -) -> iris.cube.Cube: + cube: Cube, + coordinate: str | DimCoord | AuxCoord, +) -> Cube: """Weight data using the bounds from a given coordinate. The resulting cube will then have units given by @@ -137,7 +91,7 @@ def accumulate_coordinate( Parameters ---------- cube: - Data cube for the flux + Data cube for the flux. coordinate: Name of the coordinate that will be used as weights. @@ -145,7 +99,7 @@ def accumulate_coordinate( Returns ------- iris.cube.Cube - Cube with the aggregated data + Cube with the aggregated data. Raises ------ @@ -169,7 +123,7 @@ def accumulate_coordinate( ) array_module = da if coord.has_lazy_bounds() else np - factor = iris.coords.AuxCoord( + factor = AuxCoord( array_module.diff(coord.core_bounds())[..., -1], var_name=coord.var_name, long_name=coord.long_name, diff --git a/tests/unit/test_iris_helpers.py b/tests/unit/test_iris_helpers.py index 1b5066ae51..fe56d2d068 100644 --- a/tests/unit/test_iris_helpers.py +++ b/tests/unit/test_iris_helpers.py @@ -18,7 +18,7 @@ DimCoord, ) from iris.cube import Cube, CubeList -from iris.exceptions import CoordinateMultiDimError +from iris.exceptions import CoordinateMultiDimError, UnitConversionError from esmvalcore.iris_helpers import ( add_leading_dim_to_cube, @@ -28,6 +28,7 @@ has_unstructured_grid, merge_cube_attributes, rechunk_cube, + safe_convert_units, ) @@ -571,3 +572,59 @@ def test_has_unstructured_grid_true(lat_coord_1d, lon_coord_1d): aux_coords_and_dims=[(lat_coord_1d, 0), (lon_coord_1d, 0)], ) assert has_unstructured_grid(cube) is True + + +@pytest.mark.parametrize( + "old_units,new_units,old_standard_name,new_standard_name,err_msg", + [ + ("m", "km", "altitude", "altitude", None), + ("Pa", "hPa", "air_pressure", "air_pressure", None), + ( + "m", + "DU", + "equivalent_thickness_at_stp_of_atmosphere_ozone_content", + "equivalent_thickness_at_stp_of_atmosphere_ozone_content", + None, + ), + ( + "m", + "s", + "altitude", + ValueError, + r"Unable to convert from 'Unit\('m'\)' to 'Unit\('s'\)'", + ), + ( + "unknown", + "s", + "air_temperature", + UnitConversionError, + r"Cannot convert from unknown units", + ), + ( + "kg m-2 s-1", + "mm day-1", + "precipitation_flux", + ValueError, + r"Cannot safely convert units from 'kg m-2 s-1' to 'mm day-1'; " + r"standard_name changed from 'precipitation_flux' to " + r"'lwe_precipitation_rate'", + ), + ], +) +def test_safe_convert_units( + old_units, new_units, old_standard_name, new_standard_name, err_msg +): + """Test ``esmvalcore.preprocessor._units.safe_convert_units``.""" + cube = Cube(0, standard_name=old_standard_name, units=old_units) + + # Exceptions + if isinstance(new_standard_name, type): + with pytest.raises(new_standard_name, match=err_msg): + safe_convert_units(cube, new_units) + return + + # Regular test cases + new_cube = safe_convert_units(cube, new_units) + assert new_cube is cube + assert new_cube.standard_name == new_standard_name + assert new_cube.units == new_units