Repository: ESMValGroup/ESMValCore

Issue:
New preprocessor: subtract global mean from field
A common operation in climate science is to subtract global means from a field. I will implement that.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
> Merging [#2189](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (3ee4d44) into [main](https://app.codecov.io/gh/ESMValGroup/ESMValCore/commit/4baa8aadc8afcbe12aa21f566666c0de11bc72fa?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (4baa8aa) will **increase** coverage by `0.02%`.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff             @@
##             main    #2189      +/-   ##
==========================================
+ Coverage   93.14%   93.16%   +0.02%     
==========================================
  Files         237      237              
  Lines       12828    12858      +30     
==========================================
+ Hits        11948    11979      +31     
+ Misses        880      879       -1     
```


| [Files Changed](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) | Coverage Δ | |
|---|---|---|
| [esmvalcore/preprocessor/\_\_init\_\_.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9wcmVwcm9jZXNzb3IvX19pbml0X18ucHk=) | `95.57% <ø> (+0.34%)` | :arrow_up: |
| [esmvalcore/\_recipe/check.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9fcmVjaXBlL2NoZWNrLnB5) | `92.45% <100.00%> (+0.42%)` | :arrow_up: |
| [esmvalcore/\_recipe/recipe.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9fcmVjaXBlL3JlY2lwZS5weQ==) | `99.01% <100.00%> (+<0.01%)` | :arrow_up: |
| [esmvalcore/preprocessor/\_area.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2189?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9wcmVwcm9jZXNzb3IvX2FyZWEucHk=) | `94.78% <100.00%> (+0.23%)` | :arrow_up: |


valeriupredoi: @bouweandela would you be a darling a merge this, please? (I've been merging your stuff all morning :rofl: )
bouweandela: @schlunma Would it make sense to apply the same changes here as in #2191? I.e. remove the nesting from the function signature and it to:
```python
def normalize_with_area_statistics(
    cube: Cube,
    operation: str = 'subtract',
    operator: str,
    **kwargs: Any,
) -> Cube:
```
? (Or, if you think having bot operation and operator is confusing, rename one of the two?)
schlunma: I honestly think it is better the way it is. The keywords here are only relevant for the `area_statistics` part of this preprocessor, whereas in #2191 they are relevant for the actual preprocessors.

I have seen this kind of API in other packages as well, e.g., the `line_kws` in [`seaborn.regplot`](https://seaborn.pydata.org/generated/seaborn.regplot.html) which are only relevant for a part of the function.


schlunma: I would be nice to merge this PR after #2191; I want to add a proper check for the `area_statistics` keywords here introduced with #2191.
bouweandela: Thinking ahead a bit, this kind of feature could be implemented by adding just one extra keyword argument to all of the current `*_statistics` preprocessor functions. Maybe people would also be interested in knowing the spread of an ensemble for the `ensemble_statistics` preprocessor function, or normalizing by `volume_statistics`, etc. If we add new preprocessor functions for all of those we end up with many preprocessor functions. Are you really keen on getting this included into v2.10 (today is the deadline for the feature freeze) or would you be willing to think about it a bit more?
schlunma: I like the idea of generalizing this, but I don't think we should simply add a keyword to the existing `statistics` preprocessors. I think it would be weird if `area_statistics` somehow returned a cube that still has `latitude/longitude` dimensions.

It would probably be clearer to add a `normalize_with_statistics` preprocessor that takes the statistic as argument. Probably easy to implement for most preprocessors except for the multi-model ones since their call signature is different.  I can try to implement this for all preprocessors except for the multi-model ones today. Would that work for you?


bouweandela: That would probably be nicer than this approach, right?
schlunma: One issue with this generalized approach is that `normalize_with_statistics` requires different supplementary variables dependent on the statistic used (this is given as an argument to the function). Currently, the type of supplementary variable that a preperocessor uses is given with a [decorator](https://github.com/ESMValGroup/ESMValCore/blob/87e3f9fe95d2496222fb4029ce1e13feed1752b1/esmvalcore/preprocessor/_normalize.py#L18-L21), which knows nothing about the preprocessor arguments. In addition, if this `normalize_with_statistics` is used multiple times in the recipe, it might require different supplementary variables. This is also a problem for the current implementation, since the necessary supplementaries are [stored in a dictionary with the name of the preprocessor as key](https://github.com/ESMValGroup/ESMValCore/blob/87e3f9fe95d2496222fb4029ce1e13feed1752b1/esmvalcore/preprocessor/_supplementary_vars.py#L36).

@bouweandela any ideas how to solve this smartly without changing our entire code base?
bouweandela: In that case, adding the normalize option as an extra keyword argument to existing preprocessor functions may be a more attractive solution because it avoids all these technical issues. I agree that it may feel a bit counter-intuitive that the returned cube can have different dimensions based on the arguments though. Maybe we could ask the community for their opinion at the next monthly community meeting?
bouweandela: I just came across discussion https://github.com/ESMValGroup/ESMValTool/discussions/3353 where this seems relevant.
schlunma: > In that case, adding the normalize option as an extra keyword argument to existing preprocessor functions may be a more attractive solution because it avoids all these technical issues.

Agreed, that would be much simpler. I'll try to implement something like this for v2.11 because I need this new feature.

> I just came across discussion [ESMValGroup/ESMValTool#3353](https://github.com/ESMValGroup/ESMValTool/discussions/3353) where this seems relevant.

I don't think that's relevant if I understand the problem correctly (weighting != normalization).


valeriupredoi: oh sorry, Manu - I got alerted I am up for providing a review - just saw it's still a draft - carry on, then
schlunma: This is ready for review now
valeriupredoi: not willing to approve this yet another time, I'll pass the ball to @bouweandela - if you are happy with all the stuff here, please merge, bud :beer: 
valeriupredoi: frenly ping @bouweandela :baby_chick: 

Patch:
diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst
index 5c994ff00..638a9c6c1 100644
--- a/doc/recipe/preprocessor.rst
+++ b/doc/recipe/preprocessor.rst
@@ -1962,11 +1962,15 @@ along the longitude coordinate.
 Parameters:
     * `operator`: Operation to apply.
       See :ref:`stat_preprocs` for more details on supported statistics.
+    * `normalize`: If given, do not return the statistics cube itself, but
+      rather, the input cube, normalized with the statistics cube. Can either
+      be `subtract` (statistics cube is subtracted from the input cube) or
+      `divide` (input cube is divided by the statistics cube).
     * Other parameters are directly passed to the `operator` as keyword
       arguments.
       See :ref:`stat_preprocs` for more details.
 
-See also :func:`esmvalcore.preprocessor.zonal_means`.
+See also :func:`esmvalcore.preprocessor.zonal_statistics`.
 
 
 ``meridional_statistics``
@@ -1979,11 +1983,15 @@ argument:
 Parameters:
     * `operator`: Operation to apply.
       See :ref:`stat_preprocs` for more details on supported statistics.
+    * `normalize`: If given, do not return the statistics cube itself, but
+      rather, the input cube, normalized with the statistics cube. Can either
+      be `subtract` (statistics cube is subtracted from the input cube) or
+      `divide` (input cube is divided by the statistics cube).
     * Other parameters are directly passed to the `operator` as keyword
       arguments.
       See :ref:`stat_preprocs` for more details.
 
-See also :func:`esmvalcore.preprocessor.meridional_means`.
+See also :func:`esmvalcore.preprocessor.meridional_statistics`.
 
 
 .. _area_statistics:
@@ -1991,10 +1999,6 @@ See also :func:`esmvalcore.preprocessor.meridional_means`.
 ``area_statistics``
 -------------------
 
-This function calculates statistics over a region.
-It takes one argument, ``operator``, which is the name of the operation to
-apply.
-
 This function can be used to apply several different operations in the
 horizontal plane: for example, mean, sum, standard deviation, median, variance,
 minimum, maximum and root mean square.
@@ -2013,6 +2017,33 @@ The required supplementary variable, either ``areacella`` for atmospheric
 variables or ``areacello`` for ocean variables, can be attached to the main
 dataset as described in :ref:`supplementary_variables`.
 
+Parameters:
+    * `operator`: Operation to apply.
+      See :ref:`stat_preprocs` for more details on supported statistics.
+    * `normalize`: If given, do not return the statistics cube itself, but
+      rather, the input cube, normalized with the statistics cube. Can either
+      be `subtract` (statistics cube is subtracted from the input cube) or
+      `divide` (input cube is divided by the statistics cube).
+    * Other parameters are directly passed to the `operator` as keyword
+      arguments.
+      See :ref:`stat_preprocs` for more details.
+
+Examples:
+* Calculate global mean:
+
+  .. code-block:: yaml
+
+    area_statistics:
+      operator: mean
+
+* Subtract global mean from dataset:
+
+  .. code-block:: yaml
+
+    area_statistics:
+      operator: mean
+      normalize: subtract
+
 See also :func:`esmvalcore.preprocessor.area_statistics`.
 
 
@@ -2100,9 +2131,6 @@ See also :func:`esmvalcore.preprocessor.extract_volume`.
 This function calculates the volume-weighted average across three dimensions,
 but maintains the time dimension.
 
-This function takes the argument: `operator`, which defines the operation to
-apply over the volume.
-At the moment, only `mean` is supported.
 By default, the `mean` operation is weighted by the grid cell volumes.
 
 For weighted statistics, this function requires a cell volume `cell measure`_,
@@ -2113,6 +2141,18 @@ dataset as described in :ref:`supplementary_variables`.
 
 No depth coordinate is required as this is determined by Iris.
 
+Parameters:
+    * `operator`: Operation to apply.
+      At the moment, only `mean` is supported.
+      See :ref:`stat_preprocs` for more details on supported statistics.
+    * `normalize`: If given, do not return the statistics cube itself, but
+      rather, the input cube, normalized with the statistics cube. Can either
+      be `subtract` (statistics cube is subtracted from the input cube) or
+      `divide` (input cube is divided by the statistics cube).
+    * Other parameters are directly passed to the `operator` as keyword
+      arguments.
+      See :ref:`stat_preprocs` for more details.
+
 See also :func:`esmvalcore.preprocessor.volume_statistics`.
 
 .. _axis_statistics:
@@ -2128,6 +2168,10 @@ Takes arguments:
       Possible values for the axis are `x`, `y`, `z`, `t`.
     * `operator`: Operation to apply.
       See :ref:`stat_preprocs` for more details on supported statistics.
+    * `normalize`: If given, do not return the statistics cube itself, but
+      rather, the input cube, normalized with the statistics cube. Can either
+      be `subtract` (statistics cube is subtracted from the input cube) or
+      `divide` (input cube is divided by the statistics cube).
     * Other parameters are directly passed to the `operator` as keyword
       arguments.
       See :ref:`stat_preprocs` for more details.
diff --git a/esmvalcore/preprocessor/_area.py b/esmvalcore/preprocessor/_area.py
index 145b87ec3..7d7bbe0c9 100644
--- a/esmvalcore/preprocessor/_area.py
+++ b/esmvalcore/preprocessor/_area.py
@@ -8,7 +8,7 @@ from __future__ import annotations
 import logging
 import warnings
 from pathlib import Path
-from typing import TYPE_CHECKING, Iterable, Optional
+from typing import TYPE_CHECKING, Iterable, Literal, Optional
 
 import fiona
 import iris
@@ -20,8 +20,13 @@ from iris.coords import AuxCoord, CellMeasure
 from iris.cube import Cube, CubeList
 from iris.exceptions import CoordinateMultiDimError, CoordinateNotFoundError
 
-from ._shared import get_iris_aggregator, guess_bounds, update_weights_kwargs
-from ._supplementary_vars import (
+from esmvalcore.preprocessor._shared import (
+    get_iris_aggregator,
+    get_normalized_cube,
+    guess_bounds,
+    update_weights_kwargs,
+)
+from esmvalcore.preprocessor._supplementary_vars import (
     add_ancillary_variable,
     add_cell_measure,
     register_supplementaries,
@@ -187,6 +192,7 @@ def _extract_irregular_region(cube, start_longitude, end_longitude,
 def zonal_statistics(
     cube: Cube,
     operator: str,
+    normalize: Optional[Literal['subtract', 'divide']] = None,
     **operator_kwargs
 ) -> Cube:
     """Compute zonal statistics.
@@ -199,6 +205,11 @@ def zonal_statistics(
         The operation. Used to determine the :class:`iris.analysis.Aggregator`
         object used to calculate the statistics. Allowed options are given in
         :ref:`this table <supported_stat_operator>`.
+    normalize:
+        If given, do not return the statistics cube itself, but rather, the
+        input cube, normalized with the statistics cube. Can either be
+        `subtract` (statistics cube is subtracted from the input cube) or
+        `divide` (input cube is divided by the statistics cube).
     **operator_kwargs:
         Optional keyword arguments for the :class:`iris.analysis.Aggregator`
         object defined by `operator`.
@@ -206,7 +217,8 @@ def zonal_statistics(
     Returns
     -------
     iris.cube.Cube
-        Zonal statistics cube.
+        Zonal statistics cube or input cube normalized by statistics cube (see
+        `normalize`).
 
     Raises
     ------
@@ -220,14 +232,17 @@ def zonal_statistics(
             "Zonal statistics on irregular grids not yet implemented"
         )
     (agg, agg_kwargs) = get_iris_aggregator(operator, **operator_kwargs)
-    cube = cube.collapsed('longitude', agg, **agg_kwargs)
-    cube.data = cube.core_data().astype(np.float32, casting='same_kind')
-    return cube
+    result = cube.collapsed('longitude', agg, **agg_kwargs)
+    if normalize is not None:
+        result = get_normalized_cube(cube, result, normalize)
+    result.data = result.core_data().astype(np.float32, casting='same_kind')
+    return result
 
 
 def meridional_statistics(
     cube: Cube,
     operator: str,
+    normalize: Optional[Literal['subtract', 'divide']] = None,
     **operator_kwargs,
 ) -> Cube:
     """Compute meridional statistics.
@@ -240,6 +255,11 @@ def meridional_statistics(
         The operation. Used to determine the :class:`iris.analysis.Aggregator`
         object used to calculate the statistics. Allowed options are given in
         :ref:`this table <supported_stat_operator>`.
+    normalize:
+        If given, do not return the statistics cube itself, but rather, the
+        input cube, normalized with the statistics cube. Can either be
+        `subtract` (statistics cube is subtracted from the input cube) or
+        `divide` (input cube is divided by the statistics cube).
     **operator_kwargs:
         Optional keyword arguments for the :class:`iris.analysis.Aggregator`
         object defined by `operator`.
@@ -261,9 +281,11 @@ def meridional_statistics(
             "Meridional statistics on irregular grids not yet implemented"
         )
     (agg, agg_kwargs) = get_iris_aggregator(operator, **operator_kwargs)
-    cube = cube.collapsed('latitude', agg, **agg_kwargs)
-    cube.data = cube.core_data().astype(np.float32, casting='same_kind')
-    return cube
+    result = cube.collapsed('latitude', agg, **agg_kwargs)
+    if normalize is not None:
+        result = get_normalized_cube(cube, result, normalize)
+    result.data = result.core_data().astype(np.float32, casting='same_kind')
+    return result
 
 
 def compute_area_weights(cube):
@@ -348,6 +370,7 @@ def _try_adding_calculated_cell_area(cube: Cube) -> None:
 def area_statistics(
     cube: Cube,
     operator: str,
+    normalize: Optional[Literal['subtract', 'divide']] = None,
     **operator_kwargs,
 ) -> Cube:
     """Apply a statistical operator in the horizontal plane.
@@ -370,6 +393,11 @@ def area_statistics(
         The operation. Used to determine the :class:`iris.analysis.Aggregator`
         object used to calculate the statistics. Allowed options are given in
         :ref:`this table <supported_stat_operator>`.
+    normalize:
+        If given, do not return the statistics cube itself, but rather, the
+        input cube, normalized with the statistics cube. Can either be
+        `subtract` (statistics cube is subtracted from the input cube) or
+        `divide` (input cube is divided by the statistics cube).
     **operator_kwargs:
         Optional keyword arguments for the :class:`iris.analysis.Aggregator`
         object defined by `operator`.
@@ -396,6 +424,8 @@ def area_statistics(
     )
 
     result = cube.collapsed(['latitude', 'longitude'], agg, **agg_kwargs)
+    if normalize is not None:
+        result = get_normalized_cube(cube, result, normalize)
 
     # Make sure to preserve dtype
     new_dtype = result.dtype
diff --git a/esmvalcore/preprocessor/_shared.py b/esmvalcore/preprocessor/_shared.py
index 09e98d693..56bfc4256 100644
--- a/esmvalcore/preprocessor/_shared.py
+++ b/esmvalcore/preprocessor/_shared.py
@@ -7,8 +7,9 @@ import logging
 import re
 import warnings
 from collections.abc import Callable
-from typing import Any, Optional
+from typing import Any, Literal, Optional
 
+import dask.array as da
 import iris.analysis
 import numpy as np
 from iris.coords import DimCoord
@@ -181,3 +182,59 @@ def update_weights_kwargs(
     else:
         kwargs.pop('weights', None)
     return kwargs
+
+
+def get_normalized_cube(
+    cube: Cube,
+    statistics_cube: Cube,
+    normalize: Literal['subtract', 'divide'],
+) -> Cube:
+    """Get cube normalized with statistics cube.
+
+    Parameters
+    ----------
+    cube:
+        Input cube that will be normalized.
+    statistics_cube:
+        Cube that is used to normalize the input cube. Needs to be
+        broadcastable to the input cube's shape according to iris' rich
+        broadcasting rules enabled by the use of named dimensions (see also
+        https://scitools-iris.readthedocs.io/en/latest/userguide/cube_maths.
+        html#calculating-a-cube-anomaly). This is usually ensure by using
+        :meth:`iris.cube.Cube.collapsed` to calculate the statistics cube.
+    normalize:
+        Normalization operation. Can either be `subtract` (statistics cube is
+        subtracted from the input cube) or `divide` (input cube is divided by
+        the statistics cube).
+
+    Returns
+    -------
+    Cube
+        Input cube normalized with statistics cube.
+
+    """
+    if normalize == 'subtract':
+        normalized_cube = cube - statistics_cube
+
+    elif normalize == 'divide':
+        normalized_cube = cube / statistics_cube
+
+        # Iris sometimes masks zero-divisions, sometimes not
+        # (https://github.com/SciTools/iris/issues/5523). Make sure to
+        # consistently mask them here.
+        normalized_cube.data = da.ma.masked_invalid(
+            normalized_cube.core_data()
+        )
+
+    else:
+        raise ValueError(
+            f"Expected 'subtract' or 'divide' for `normalize`, got "
+            f"'{normalize}'"
+        )
+
+    # Keep old metadata except for units
+    new_units = normalized_cube.units
+    normalized_cube.metadata = cube.metadata
+    normalized_cube.units = new_units
+
+    return normalized_cube
diff --git a/esmvalcore/preprocessor/_volume.py b/esmvalcore/preprocessor/_volume.py
index e0c376585..8eb03e11d 100644
--- a/esmvalcore/preprocessor/_volume.py
+++ b/esmvalcore/preprocessor/_volume.py
@@ -7,7 +7,7 @@ from __future__ import annotations
 
 import logging
 import warnings
-from typing import Iterable, Sequence
+from typing import Iterable, Literal, Optional, Sequence
 
 import dask.array as da
 import iris
@@ -17,7 +17,11 @@ from iris.cube import Cube
 from iris.exceptions import CoordinateMultiDimError
 
 from ._area import compute_area_weights
-from ._shared import get_iris_aggregator, update_weights_kwargs
+from ._shared import (
+    get_iris_aggregator,
+    get_normalized_cube,
+    update_weights_kwargs,
+)
 from ._supplementary_vars import register_supplementaries
 
 logger = logging.getLogger(__name__)
@@ -182,6 +186,7 @@ def _try_adding_calculated_ocean_volume(cube: Cube) -> None:
 def volume_statistics(
     cube: Cube,
     operator: str,
+    normalize: Optional[Literal['subtract', 'divide']] = None,
     **operator_kwargs,
 ) -> Cube:
     """Apply a statistical operation over a volume.
@@ -201,6 +206,11 @@ def volume_statistics(
         The operation. Used to determine the :class:`iris.analysis.Aggregator`
         object used to calculate the statistics. Currently, only `mean` is
         allowed.
+    normalize:
+        If given, do not return the statistics cube itself, but rather, the
+        input cube, normalized with the statistics cube. Can either be
+        `subtract` (statistics cube is subtracted from the input cube) or
+        `divide` (input cube is divided by the statistics cube).
     **operator_kwargs:
         Optional keyword arguments for the :class:`iris.analysis.Aggregator`
         object defined by `operator`.
@@ -232,6 +242,8 @@ def volume_statistics(
         agg,
         **agg_kwargs,
     )
+    if normalize is not None:
+        result = get_normalized_cube(cube, result, normalize)
 
     # Make sure input cube has not been modified
     if not has_cell_measure and cube.cell_measures('ocean_volume'):
@@ -244,6 +256,7 @@ def axis_statistics(
     cube: Cube,
     axis: str,
     operator: str,
+    normalize: Optional[Literal['subtract', 'divide']] = None,
     **operator_kwargs,
 ) -> Cube:
     """Perform statistics along a given axis.
@@ -267,6 +280,11 @@ def axis_statistics(
         The operation. Used to determine the :class:`iris.analysis.Aggregator`
         object used to calculate the statistics. Allowed options are given in
         :ref:`this table <supported_stat_operator>`.
+    normalize:
+        If given, do not return the statistics cube itself, but rather, the
+        input cube, normalized with the statistics cube. Can either be
+        `subtract` (statistics cube is subtracted from the input cube) or
+        `divide` (input cube is divided by the statistics cube).
     **operator_kwargs:
         Optional keyword arguments for the :class:`iris.analysis.Aggregator`
         object defined by `operator`.
@@ -320,6 +338,9 @@ def axis_statistics(
         )
         result = cube.collapsed(coord, agg, **agg_kwargs)
 
+    if normalize is not None:
+        result = get_normalized_cube(cube, result, normalize)
+
     # Make sure input and output cubes do not have auxiliary coordinate
     if cube.coords('_axis_statistics_weights_'):
         cube.remove_coord('_axis_statistics_weights_')


Test patch:
diff --git a/tests/unit/preprocessor/_area/test_area.py b/tests/unit/preprocessor/_area/test_area.py
index 3939437f3..8835014ff 100644
--- a/tests/unit/preprocessor/_area/test_area.py
+++ b/tests/unit/preprocessor/_area/test_area.py
@@ -188,6 +188,41 @@ class Test(tests.Test):
         self.assert_array_equal(result.data, expected)
         self.assertEqual(result.units, 'kg m-2 s-1')
 
+    def test_area_statistics_subtract_mean(self):
+        """Test for area average of a 2D field."""
+        input_data = self.grid.copy()
+        self.assertFalse(input_data.cell_measures('cell_area'))
+
+        result = area_statistics(input_data, 'mean', normalize='subtract')
+
+        self.assertEqual(input_data, self.grid)
+        self.assertEqual(result.shape, input_data.shape)
+        expected = np.ma.zeros((2, 5, 5), dtype=np.float32)
+        self.assert_array_equal(result.data, expected)
+        self.assertFalse(result.cell_measures('cell_area'))
+        self.assertEqual(result.metadata, self.grid.metadata)
+        for coord in self.grid.coords():
+            self.assertEqual(result.coord(coord.name()), coord)
+
+    def test_area_statistics_cell_measure_subtract_mean(self):
+        """Test for area average of a 2D field.
+
+        The area measure is pre-loaded in the cube.
+        """
+        self._add_cell_measure_to_grid()
+        input_data = self.grid.copy()
+
+        result = area_statistics(input_data, 'mean', normalize='subtract')
+
+        self.assertEqual(input_data, self.grid)
+        self.assertEqual(result.shape, input_data.shape)
+        expected = np.ma.zeros((2, 5, 5), dtype=np.float32)
+        self.assert_array_equal(result.data, expected)
+        self.assertFalse(result.cell_measures('cell_area'))
+        self.assertEqual(result.metadata, self.grid.metadata)
+        for coord in self.grid.coords():
+            self.assertEqual(result.coord(coord.name()), coord)
+
     def test_extract_region(self):
         """Test for extracting a region from a 2D field."""
         result = extract_region(self.grid, 1.5, 2.5, 1.5, 2.5)
@@ -1270,6 +1305,39 @@ def test_zonal_statistics(make_testcube):
     assert res.dtype == np.float32
 
 
+def test_zonal_statistics_divide_by_min(make_testcube):
+    """Test ``zonal_statistics``."""
+    make_testcube.data = np.ones(make_testcube.shape, dtype=np.float32)
+    make_testcube.data[0, 0] = 0.0
+    make_testcube.data[1, 0] = -1.0
+    make_testcube.data[2, 0] = -0.5
+    make_testcube.units = 'K'
+    input_data = make_testcube.copy()
+
+    res = zonal_statistics(input_data, 'min', normalize='divide')
+
+    assert input_data == make_testcube
+    assert res.shape == input_data.shape
+    expected = np.ma.masked_invalid(
+        [
+            [np.nan, np.nan, np.nan, np.nan, np.nan],
+            [1.0, -1.0, -1.0, -1.0, -1.0],
+            [1.0, -2.0, -2.0, -2.0, -2.0],
+            [1.0, 1.0, 1.0, 1.0, 1.0],
+            [1.0, 1.0, 1.0, 1.0, 1.0],
+        ],
+    )
+    tests.assert_array_equal(res.data, expected)
+    assert res.dtype == np.float32
+
+    assert res.standard_name == input_data.standard_name
+    assert res.var_name == input_data.var_name
+    assert res.long_name == input_data.long_name
+    assert res.cell_methods == input_data.cell_methods
+    assert res.attributes == input_data.attributes
+    assert res.units == '1'
+
+
 def test_zonal_statistics_2d_lon_fail(irreg_extract_shape_cube):
     """Test ``zonal_statistics``."""
     with pytest.raises(ValueError):
@@ -1286,11 +1354,50 @@ def test_meridional_statistics(make_testcube):
     assert res.dtype == np.float32
 
 
+def test_meridional_statistics_divide_by_max(make_testcube):
+    """Test ``meridional_statistics``."""
+    make_testcube.data = np.ones(make_testcube.shape, dtype=np.float32)
+    make_testcube.data[0, 0] = 0.25
+    make_testcube.data[0, 1] = 2.0
+    make_testcube.units = 'K'
+    input_data = make_testcube.copy()
+
+    res = meridional_statistics(input_data, 'max', normalize='divide')
+
+    assert input_data == make_testcube
+    assert res.shape == input_data.shape
+    expected = np.ma.masked_invalid(
+        [
+            [0.25, 1.0, 1.0, 1.0, 1.0],
+            [1.0, 0.5, 1.0, 1.0, 1.0],
+            [1.0, 0.5, 1.0, 1.0, 1.0],
+            [1.0, 0.5, 1.0, 1.0, 1.0],
+            [1.0, 0.5, 1.0, 1.0, 1.0],
+        ],
+    )
+    tests.assert_array_equal(res.data, expected)
+    assert res.dtype == np.float32
+
+    assert res.standard_name == input_data.standard_name
+    assert res.var_name == input_data.var_name
+    assert res.long_name == input_data.long_name
+    assert res.cell_methods == input_data.cell_methods
+    assert res.attributes == input_data.attributes
+    assert res.units == '1'
+
+
 def test_meridional_statistics_2d_lon_fail(irreg_extract_shape_cube):
     """Test ``meridional_statistics``."""
     with pytest.raises(ValueError):
         meridional_statistics(irreg_extract_shape_cube, 'sum')
 
 
+def test_meridional_statistics_invalid_norm_fail(make_testcube):
+    """Test ``meridional_statistics``."""
+    msg = "Expected 'subtract' or 'divide' for `normalize`"
+    with pytest.raises(ValueError, match=msg):
+        meridional_statistics(make_testcube, 'sum', normalize='x')
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/unit/preprocessor/_volume/test_volume.py b/tests/unit/preprocessor/_volume/test_volume.py
index 2e35f98d3..f0b730561 100644
--- a/tests/unit/preprocessor/_volume/test_volume.py
+++ b/tests/unit/preprocessor/_volume/test_volume.py
@@ -173,6 +173,34 @@ class Test(tests.Test):
         self.assert_array_equal(result.data, expected)
         self.assertEqual(result.units, 'kg m-2')
 
+    def test_axis_statistics_subtract_mean(self):
+        """Test axis statistics with operator mean."""
+        data = np.ma.arange(1, 25).reshape(2, 3, 2, 2)
+        self.grid_4d.data = data
+        result = axis_statistics(
+            self.grid_4d, 'z', 'mean', normalize='subtract'
+        )
+        bounds = self.grid_4d.coord(axis='z').bounds
+        weights = (bounds[:, 1] - bounds[:, 0])
+        expected = (
+            data - np.average(data, axis=1, weights=weights, keepdims=True)
+        )
+        self.assert_array_equal(result.data, expected)
+        self.assertEqual(result.units, 'kg m-3')
+        self.assertFalse(self.grid_4d.coords('_axis_statistics_weights_'))
+        self.assertFalse(result.coords('_axis_statistics_weights_'))
+
+    def test_axis_statistics_divide_min(self):
+        """Test axis statistics with operator sum."""
+        data = np.ma.arange(1, 25).reshape(2, 3, 2, 2)
+        self.grid_4d.data = data
+        result = axis_statistics(self.grid_4d, 'z', 'min', normalize='divide')
+        expected = data / np.min(data, axis=1, keepdims=True)
+        self.assert_array_equal(result.data, expected)
+        self.assertEqual(result.units, '1')
+        self.assertFalse(self.grid_4d.coords('_axis_statistics_weights_'))
+        self.assertFalse(result.coords('_axis_statistics_weights_'))
+
     def test_wrong_axis_statistics_fail(self):
         """Test raises error when axis is not found in cube."""
         with self.assertRaises(ValueError) as err:
@@ -360,6 +388,36 @@ class Test(tests.Test):
         self.assert_array_equal(result.data, expected)
         self.assertEqual(result.units, 'kg m-3')
 
+    def test_volume_statistics_subtract_mean(self):
+        """Test to take the volume weighted average of a (2,3,2,2) cube."""
+        self.assertFalse(self.grid_4d.cell_measures('ocean_volume'))
+
+        result = volume_statistics(self.grid_4d, 'mean', normalize='subtract')
+
+        expected = np.ma.zeros((2, 3, 2, 2))
+        self.assert_array_equal(result.data, expected)
+        self.assertEqual(result.units, 'kg m-3')
+        self.assertFalse(self.grid_4d.cell_measures('ocean_volume'))
+        self.assertFalse(result.cell_measures('ocean_volume'))
+
+    def test_volume_statistics_cell_measure_divide_mean(self):
+        """Test to take the volume weighted average of a (2,3,2,2) cube.
+
+        The volume measure is pre-loaded in the cube.
+        """
+        grid_volume = calculate_volume(self.grid_4d)
+        measure = iris.coords.CellMeasure(grid_volume,
+                                          standard_name='ocean_volume',
+                                          units='m3',
+                                          measure='volume')
+        self.grid_4d.add_cell_measure(measure, range(0, measure.ndim))
+
+        result = volume_statistics(self.grid_4d, 'mean', normalize='divide')
+
+        expected = np.ma.ones((2, 3, 2, 2))
+        self.assert_array_equal(result.data, expected)
+        self.assertEqual(result.units, '1')
+
     def test_volume_statistics_weights(self):
         """Test to take the volume weighted average of a (2,3,2,2) cube.