Merge pull request #218 from ESMValGroup/detrend

Add detrend preprocessor function

doc/esmvalcore/preprocessor.rst

@@ -17,6 +17,7 @@ following the default order in which they are applied:
 * :ref:`Time operations`
 * :ref:`Area operations`
 * :ref:`Volume operations`
+* :ref:`Detrend`
 * :ref:`Unit conversion`
 
 Overview
@@ -859,6 +860,28 @@ Note that this function uses the expensive ``interpolate`` method from
 
 See also :func:`esmvalcore.preprocessor.extract_trajectory`.
 
+.. _detrend:
+
+Detrend
+=======
+
+ESMValTool also supports detrending along any dimension using
+the preprocessor function 'detrend'.
+This function has two parameters:
+
+* ``dimension``: dimension to apply detrend on. Default: "time"
+* ``method``: It can be ``linear`` or ``constant``. Default: ``linear``
+
+If method is ``linear``, detrend will calculate the linear trend along the
+selected axis and substract it to the data. For example, this can be used to
+remove the linear trend caused by climate change on some variables is selected
+dimension is time.
+
+If method is ``constant``, detrend will compute the mean along that dimension
+and substract it from the data
+
+See also :func:`esmvalcore.preprocessor.detrend`.
+
 .. _unit conversion:
 
 Unit conversion

esmvalcore/preprocessor/__init__.py

@@ -10,6 +10,7 @@
 from ._area import (area_statistics, extract_named_regions, extract_region,
                     zonal_means)
 from ._derive import derive
+from ._detrend import detrend
 from ._download import download
 from ._io import (_get_debug_filename, cleanup, concatenate, load, save,
                   write_metadata)
@@ -68,6 +69,7 @@
     'extract_transect',
     # 'average_zone': average_zone,
     # 'cross_section': cross_section,
+    'detrend',
     'multi_model_statistics',
     # Grid-point operations
     'extract_named_regions',

esmvalcore/preprocessor/_detrend.py

@@ -0,0 +1,38 @@
+"""Preprocessor functions that remove trends from the data."""
+import logging
+
+import dask.array as da
+import scipy.signal
+
+logger = logging.getLogger(__name__)
+
+
+def detrend(cube, dimension='time', method='linear'):
+    """
+    Detrend data along a given dimension.
+
+    Parameters
+    ----------
+    cube: iris.cube.Cube
+        input cube.
+    dimension: str
+        Dimension to detrend
+    method: str
+        Method to detrend. Available: linear, constant. See documentation of
+        'scipy.signal.detrend' for details
+
+    Returns
+    -------
+    iris.cube.Cube
+        Detrended cube
+    """
+    coord = cube.coord(dimension)
+    axis = cube.coord_dims(coord)[0]
+    detrended = da.apply_along_axis(
+        scipy.signal.detrend,
+        axis=axis,
+        arr=cube.lazy_data(),
+        type=method,
+        shape=(cube.shape[axis],)
+    )
+    return cube.copy(detrended)

tests/unit/preprocessor/_detrend/test_detrend.py

@@ -0,0 +1,57 @@
+"""Unit tests for the :func:`esmvalcore.preprocessor._detrend` module."""
+
+import unittest
+
+import iris
+import iris.coords
+from iris.cube import Cube
+import numpy as np
+import pytest
+from cf_units import Unit
+
+from numpy.testing import assert_array_almost_equal
+
+from esmvalcore.preprocessor._detrend import detrend
+
+
+def _create_sample_cube():
+    cube = Cube(
+        np.array((np.arange(1, 25), np.arange(25, 49))),
+        var_name='co2',
+        units='J'
+    )
+    cube.add_dim_coord(
+        iris.coords.DimCoord(
+            np.arange(15., 720., 30.),
+            standard_name='time',
+            units=Unit('days since 1950-01-01 00:00:00', calendar='gregorian'),
+        ),
+        1,
+    )
+    cube.add_dim_coord(
+        iris.coords.DimCoord(
+            np.arange(1, 3),
+            standard_name='latitude',
+        ),
+        0,
+    )
+    return cube
+
+
+@pytest.mark.parametrize('method', ['linear', 'constant'])
+def test_decadal_average(method):
+    """Test for decadal average."""
+    cube = _create_sample_cube()
+
+    result = detrend(cube, 'time', method)
+    if method == 'linear':
+        expected = np.zeros([2, 24])
+    else:
+        expected = np.array(
+            (np.arange(1, 25) - 12.5, np.arange(25, 49) - 36.5)
+        )
+    assert_array_almost_equal(result.data, expected)
+
+
+if __name__ == '__main__':
+    unittest.main()