Extend preprocessor anomalies

doc/recipe/preprocessor.rst

@@ -1294,7 +1294,7 @@ The ``_time.py`` module contains the following preprocessor functions:
 * climate_statistics_: Compute statistics for the full period
 * resample_time_: Resample data
 * resample_hours_: Convert between N-hourly frequencies by resampling
-* anomalies_: Compute (standardized) anomalies
+* anomalies_: Compute (standardized or relative) anomalies
 * regrid_time_: Aligns the time coordinate of each dataset, against a standardized time axis.
 * timeseries_filter_: Allows application of a filter to the time-series data.
 * local_solar_time_: Convert cube with UTC time to local solar time.
@@ -1683,6 +1683,7 @@ Parameters:
     * reference: Time slice to use as the reference to compute the climatology
       on. Can be 'null' to use the full cube or a dictionary with the
       parameters from extract_time_. Default is null
+	* relative: if true, calculate relative (in percent) anomalies (default: false)
     * standardize: if true calculate standardized anomalies (default: false)
     * seasons: if period 'seasonal' or 'season' allows to set custom seasons.
       Default is '[DJF, MAM, JJA, SON]'

esmvalcore/preprocessor/_time.py

@@ -951,6 +951,7 @@ def anomalies(
     period: str,
     reference: dict | None = None,
     standardize: bool = False,
+    relative: bool = False,
     seasons: Iterable[str] = ("DJF", "MAM", "JJA", "SON"),
 ) -> Cube:
     """Compute anomalies using a mean with the specified granularity.
@@ -970,6 +971,8 @@ def anomalies(
         Period of time to use a reference, as needed for the
         :func:`~esmvalcore.preprocessor.extract_time` preprocessor function.
         If ``None``, all available data is used as a reference.
+    relative: optional
+        If ``True`` relative anomalies are calculated.
     standardize: optional
         If ``True`` standardized anomalies are calculated.
     seasons: optional
@@ -1002,33 +1005,60 @@ def anomalies(
             )
             cube = cube / cube_stddev
             cube.units = "1"
+        elif relative:
+            cube = cube / reference * 100.0
+            cube.metadata = metadata
+            cube.units = "%"
         return cube
 
     cube = _compute_anomalies(cube, reference, period, seasons)
 
-    # standardize the results if requested
+    # standardize results or compute relative anomalies if requested
+    if standardize or relative:
+        cube = _apply_scaling(cube, period, standardize, relative)
+
+    return cube
+
+
+def _apply_scaling(
+    cube: Cube,
+    period: str,
+    standardize: bool,
+    relative: bool,
+) -> Cube:
+    """Apply standardization or relative scaling."""
     if standardize:
-        cube_stddev = climate_statistics(
-            cube,
-            operator="std_dev",
-            period=period,
+        oper = "std_dev"
+    elif relative:
+        oper = "mean"
+
+    cube_div = climate_statistics(
+        cube,
+        operator=oper,
+        period=period,
+    )
+    tdim = cube.coord_dims("time")[0]
+    reps = cube.shape[tdim] / cube_div.shape[tdim]
+    if reps % 1 != 0:
+        msg = (
+            "Cannot safely apply preprocessor to this dataset, "
+            "since the full time period of this dataset is not "
+            f"a multiple of the period '{period}'"
         )
-        tdim = cube.coord_dims("time")[0]
-        reps = cube.shape[tdim] / cube_stddev.shape[tdim]
-        if reps % 1 != 0:
-            msg = (
-                "Cannot safely apply preprocessor to this dataset, "
-                "since the full time period of this dataset is not "
-                f"a multiple of the period '{period}'"
-            )
-            raise ValueError(
-                msg,
-            )
-        cube.data = cube.core_data() / da.concatenate(
-            [cube_stddev.core_data() for _ in range(int(reps))],
-            axis=tdim,
+        raise ValueError(
+            msg,
         )
+    cube.data = cube.core_data() / da.concatenate(
+        [cube_div.core_data() for _ in range(int(reps))],
+        axis=tdim,
+    )
+
+    if standardize:
         cube.units = "1"
+    elif relative:
+        cube.data = cube.core_data() * 100.0
+        cube.units = "%"
+
     return cube
 
 

tests/unit/preprocessor/_time/test_time.py

@@ -1855,6 +1855,27 @@ def make_map_data(number_years=2):
         )
 
 
+@pytest.mark.parametrize("period", ["full"])
+def test_relative_anomalies(period):
+    """Test relative ``anomalies``."""
+    cube = make_map_data(number_years=2)
+    result = anomalies(cube, period, relative=True)
+    if period == "full":
+        expected_anomalies = cube.data - np.mean(
+            cube.data,
+            axis=0,
+            keepdims=True,
+        )
+        expected_stdanomalies = expected_anomalies / np.mean(
+            cube.data,
+            axis=0,
+            keepdims=True,
+        ) * 100.
+        expected = np.ma.masked_invalid(expected_stdanomalies)
+        assert_array_equal(result.data, expected)
+        assert result.units == "%"
+
+
 @pytest.mark.parametrize("period", ["full"])
 def test_standardized_anomalies(period):
     """Test standardized ``anomalies``."""