Allowed special case for unit conversion of precipitation (`kg m-2 s-…

…1` <--> `mm day-1`) (#1574)

* Added special case for unit conversion of pr

* Only allowed special unit conversion when standard_name matches pre-defined cases

* Fixed docstring of test

* Suggestions from code review

Co-authored-by: sloosvel <[email protected]>

doc/recipe/preprocessor.rst

@@ -1966,11 +1966,33 @@ units.
 In these cases, having a unit conversion at the end of the processing
 will guarantee homogeneous input for the diagnostics.
 
+Conversion is only supported between compatible units!
+In other words, converting temperature units from ``degC`` to ``Kelvin`` works
+fine, while changing units from ``kg`` to ``m`` will not work.
+
+However, there are some well-defined exceptions from this rule in order to
+transform one quantity to another (physically related) quantity.
+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).
+Currently, the following special conversions are supported:
+
+* ``precipitation_flux`` (``kg m-2 s-1``) --
+  ``lwe_precipitation_rate`` (``mm day-1``)
+
+.. hint::
+   Names in the list correspond to ``standard_names`` of the input data.
+   Conversions are allowed from each quantity to any other quantity given in a
+   bullet point.
+   The corresponding target quantity is inferred from the desired target units.
+   In addition, any other units convertible to the ones given are also
+   supported (e.g., instead of ``mm day-1``, ``m s-1`` is also supported).
+
 .. note::
-   Conversion is only supported between compatible units! In other
-   words, converting temperature units from ``degC`` to ``Kelvin`` works
-   fine, changing precipitation units from a rate based unit to an
-   amount based unit is not supported at the moment.
+   For the transformation between the different precipitation variables, a
+   water density of ``1000 kg m-3`` is assumed.
 
 See also :func:`esmvalcore.preprocessor.convert_units`.
 

esmvalcore/preprocessor/_units.py

@@ -4,31 +4,113 @@
 """
 import logging
 
+from cf_units import Unit
 import iris
 import numpy as np
 
 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'),
+    ],
+]
+
+
+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_std_name == std_name:
+                        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.
 
+    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
+    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).
+
+    Currently, the following special conversions are supported:
+
+        * ``precipitation_flux`` (``kg m-2 s-1``) --
+          ``lwe_precipitation_rate`` (``mm day-1``)
+
+    Names in the list correspond to ``standard_names`` of the input data.
+    Conversions are allowed from each quantity to any other quantity given in a
+    bullet point. The corresponding target quantity is inferred from the
+    desired target units. In addition, any other units convertible to the ones
+    given are also supported (e.g., instead of ``mm day-1``, ``m s-1`` is also
+    supported).
+
+    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
+    cube: iris.cube.Cube
+        Input cube.
+    units: str
+        New units in udunits form.
 
     Returns
     -------
     iris.cube.Cube
         converted cube.
+
     """
-    cube.convert_units(units)
+    try:
+        cube.convert_units(units)
+    except ValueError:
+        if not _try_special_conversions(cube, units):
+            raise
+
     return cube
 
 

tests/unit/preprocessor/_units/test_convert_units.py

@@ -4,6 +4,7 @@
 
 import cf_units
 import iris
+import iris.fileformats
 import numpy as np
 
 import tests
@@ -43,6 +44,86 @@ def test_convert_compatible_units(self):
         self.assertEqual(result.units, expected_units)
         self.assert_array_equal(result.data, expected_data)
 
+    def test_convert_precipitation_flux(self):
+        """Test special conversion of precipitation_flux."""
+        self.arr.standard_name = 'precipitation_flux'
+        self.arr.units = 'kg m-2 s-1'
+        result = convert_units(self.arr, 'mm day-1')
+        self.assertEqual(result.standard_name, 'lwe_precipitation_rate')
+        self.assertEqual(result.units, 'mm day-1')
+        np.testing.assert_allclose(
+            result.data,
+            [[0.0, 86400.0], [172800.0, 259200.0]],
+        )
+
+    def test_convert_precipitation_flux_convertible(self):
+        """Test special conversion of precipitation_flux."""
+        self.arr.standard_name = 'precipitation_flux'
+        self.arr.units = 'g m-2 yr-1'
+        result = convert_units(self.arr, 'm yr-1')
+        self.assertEqual(result.standard_name, 'lwe_precipitation_rate')
+        self.assertEqual(result.units, 'm yr-1')
+        np.testing.assert_allclose(
+            result.data,
+            [[0.0, 1.0e-6], [2.0e-6, 3.0e-6]],
+        )
+
+    def test_convert_precipitation_flux_fail_invalid_name(self):
+        """Test special conversion of precipitation_flux."""
+        self.arr.units = 'kg m-2 s-1'
+        self.assertRaises(ValueError, convert_units, self.arr, 'mm day-1')
+
+    def test_convert_precipitation_flux_fail_invalid_source_units(self):
+        """Test special conversion of precipitation_flux."""
+        self.arr.standard_name = 'precipitation_flux'
+        self.assertRaises(ValueError, convert_units, self.arr, 'mm day-1')
+
+    def test_convert_precipitation_flux_fail_invalid_target_units(self):
+        """Test special conversion of precipitation_flux."""
+        self.arr.standard_name = 'precipitation_flux'
+        self.arr.units = 'kg m-2 s-1'
+        self.assertRaises(ValueError, convert_units, self.arr, 'K')
+
+    def test_convert_lwe_precipitation_rate(self):
+        """Test special conversion of lwe_precipitation_rate."""
+        self.arr.standard_name = 'lwe_precipitation_rate'
+        self.arr.units = 'mm s-1'
+        result = convert_units(self.arr, 'kg m-2 s-1')
+        self.assertEqual(result.standard_name, 'precipitation_flux')
+        self.assertEqual(result.units, 'kg m-2 s-1')
+        np.testing.assert_allclose(
+            result.data,
+            [[0.0, 1.0], [2.0, 3.0]],
+        )
+
+    def test_convert_lwe_precipitation_rate_convertible(self):
+        """Test special conversion of lwe_precipitation_rate."""
+        self.arr.standard_name = 'lwe_precipitation_rate'
+        self.arr.units = 'm yr-1'
+        result = convert_units(self.arr, 'g m-2 yr-1')
+        self.assertEqual(result.standard_name, 'precipitation_flux')
+        self.assertEqual(result.units, 'g m-2 yr-1')
+        np.testing.assert_allclose(
+            result.data,
+            [[0.0, 1.0e6], [2.0e6, 3.0e6]],
+        )
+
+    def test_convert_lwe_precipitation_rate_fail_invalid_name(self):
+        """Test special conversion of lwe_precipitation_rate."""
+        self.arr.units = 'mm s-1'
+        self.assertRaises(ValueError, convert_units, self.arr, 'kg m-2 s-1')
+
+    def test_convert_lwe_precipitation_rate_fail_invalid_source_units(self):
+        """Test special conversion of lwe_precipitation_rate."""
+        self.arr.standard_name = 'lwe_precipitation_rate'
+        self.assertRaises(ValueError, convert_units, self.arr, 'kg m-2 s-1')
+
+    def test_convert_lwe_precipitation_rate_fail_invalid_target_units(self):
+        """Test special conversion of lwe_precipitation_rate."""
+        self.arr.standard_name = 'lwe_precipitation_rate'
+        self.arr.units = 'mm s-1'
+        self.assertRaises(ValueError, convert_units, self.arr, 'K')
+
 
 class TestFluxToTotal(tests.Test):
     """Test class for _units."""