Automatically fix small deviations in vertical levels (#1177)

* Automatically fix small deviations in vertical levels

* Make rtol and atol configurable

* Refactor automatic fix

* Reduce complexity of extract_levels

* Add unit test

* Add documentation

doc/recipe/preprocessor.rst

@@ -150,6 +150,7 @@ following:
     - Incorrect units, if they can be converted to the correct ones.
     - Direction of coordinates.
     - Automatic clipping of longitude to 0 - 360 interval.
+    - Minute differences between the required and actual vertical coordinate values
 
 
 Dataset specific fixes
@@ -349,6 +350,23 @@ to height levels and vice versa using the US standard atmosphere. E.g.
 ``coordinate = air_pressure`` will convert existing height levels
 (altitude) to pressure levels (air_pressure).
 
+If the requested levels are very close to the values in the input data,
+the function will just select the available levels instead of interpolating.
+The meaning of 'very close' can be changed by providing the parameters:
+
+* ``rtol``
+    Relative tolerance for comparing the levels in the input data to the requested
+    levels. If the levels are sufficiently close, the requested levels
+    will be assigned to the vertical coordinate and no interpolation will take place.
+    The default value is 10^-7.
+* ``atol``
+    Absolute tolerance for comparing the levels in the input data to the requested
+    levels. If the levels are sufficiently close, the requested levels
+    will be assigned to the vertical coordinate and no interpolation will take place.
+    By default, `atol` will be set to 10^-7 times the mean value of
+    of the available levels.
+
+
 * See also :func:`esmvalcore.preprocessor.extract_levels`.
 * See also :func:`esmvalcore.preprocessor.get_cmor_levels`.
 

esmvalcore/cmor/check.py

@@ -784,12 +784,22 @@ def _check_requested_values(self, coord, coord_info, var_name):
         """Check requested values."""
         if coord_info.requested:
             try:
-                cmor_points = [float(val) for val in coord_info.requested]
+                cmor_points = np.array(coord_info.requested, dtype=float)
             except ValueError:
                 cmor_points = coord_info.requested
-            coord_points = list(coord.points)
+            else:
+                atol = 1e-7 * np.mean(cmor_points)
+                if (self.automatic_fixes
+                        and coord.points.shape == cmor_points.shape
+                        and np.allclose(
+                            coord.points,
+                            cmor_points,
+                            rtol=1e-7,
+                            atol=atol,
+                        )):
+                    coord.points = cmor_points
             for point in cmor_points:
-                if point not in coord_points:
+                if point not in coord.points:
                     self.report_warning(self._contain_msg, var_name,
                                         str(point), str(coord.units))
 

esmvalcore/preprocessor/_regrid.py

@@ -642,17 +642,59 @@ def _vertical_interpolate(cube, src_levels, levels, interpolation,
     return _create_cube(cube, new_data, src_levels, levels.astype(float))
 
 
-def extract_levels(cube, levels, scheme, coordinate=None):
+def parse_vertical_scheme(scheme):
+    """Parse the scheme provided for level extraction.
+
+    Parameters
+    ----------
+    scheme : str
+        The vertical interpolation scheme to use. Choose from
+        'linear',
+        'nearest',
+        'nearest_horizontal_extrapolate_vertical',
+        'linear_horizontal_extrapolate_vertical'.
+
+    Returns
+    -------
+    (str, str)
+        A tuple containing the interpolation and extrapolation scheme.
+    """
+    if scheme not in VERTICAL_SCHEMES:
+        emsg = 'Unknown vertical interpolation scheme, got {!r}. '
+        emsg += 'Possible schemes: {!r}'
+        raise ValueError(emsg.format(scheme, VERTICAL_SCHEMES))
+
+    # This allows us to put level 0. to load the ocean surface.
+    extrap_scheme = 'nan'
+    if scheme == 'nearest_horizontal_extrapolate_vertical':
+        scheme = 'nearest'
+        extrap_scheme = 'nearest'
+
+    if scheme == 'linear_horizontal_extrapolate_vertical':
+        scheme = 'linear'
+        extrap_scheme = 'nearest'
+
+    return scheme, extrap_scheme
+
+
+def extract_levels(cube,
+                   levels,
+                   scheme,
+                   coordinate=None,
+                   rtol=1e-7,
+                   atol=None):
     """Perform vertical interpolation.
 
     Parameters
     ----------
-    cube : cube
+    cube : iris.cube.Cube
         The source cube to be vertically interpolated.
-    levels : array
+    levels : ArrayLike
         One or more target levels for the vertical interpolation. Assumed
         to be in the same S.I. units of the source cube vertical dimension
-        coordinate.
+        coordinate. If the requested levels are sufficiently close to the
+        levels of the cube, cube slicing will take place instead of
+        interpolation.
     scheme : str
         The vertical interpolation scheme to use. Choose from
         'linear',
@@ -666,29 +708,28 @@ def extract_levels(cube, levels, scheme, coordinate=None):
         existing pressure levels (air_pressure) to height levels (altitude);
         'coordinate = air_pressure' will convert existing height levels
         (altitude) to pressure levels (air_pressure).
+    rtol : float
+        Relative tolerance for comparing the levels in `cube` to the requested
+        levels. If the levels are sufficiently close, the requested levels
+        will be assigned to the cube and no interpolation will take place.
+    atol : float
+        Absolute tolerance for comparing the levels in `cube` to the requested
+        levels. If the levels are sufficiently close, the requested levels
+        will be assigned to the cube and no interpolation will take place.
+        By default, `atol` will be set to 10^-7 times the mean value of
+        the levels on the cube.
 
     Returns
     -------
-    cube
+    iris.cube.Cube
+        A cube with the requested vertical levels.
+
 
     See Also
     --------
     regrid : Perform horizontal regridding.
     """
-    if scheme not in VERTICAL_SCHEMES:
-        emsg = 'Unknown vertical interpolation scheme, got {!r}. '
-        emsg += 'Possible schemes: {!r}'
-        raise ValueError(emsg.format(scheme, VERTICAL_SCHEMES))
-
-    # This allows us to put level 0. to load the ocean surface.
-    extrap_scheme = 'nan'
-    if scheme == 'nearest_horizontal_extrapolate_vertical':
-        scheme = 'nearest'
-        extrap_scheme = 'nearest'
-
-    if scheme == 'linear_horizontal_extrapolate_vertical':
-        scheme = 'linear'
-        extrap_scheme = 'nearest'
+    interpolation, extrapolation = parse_vertical_scheme(scheme)
 
     # Ensure we have a non-scalar array of levels.
     levels = np.array(levels, ndmin=1)
@@ -709,11 +750,17 @@ def extract_levels(cube, levels, scheme, coordinate=None):
     else:
         src_levels = cube.coord(axis='z', dim_coords=True)
 
-    if (src_levels.shape == levels.shape
-            and np.allclose(src_levels.points, levels)):
+    if (src_levels.shape == levels.shape and np.allclose(
+            src_levels.points,
+            levels,
+            rtol=rtol,
+            atol=1e-7 * np.mean(src_levels.points) if atol is None else atol,
+    )):
         # Only perform vertical extraction/interpolation if the source
         # and target levels are not "similar" enough.
         result = cube
+        # Set the levels to the requested values
+        src_levels.points = levels
     elif len(src_levels.shape) == 1 and \
             set(levels).issubset(set(src_levels.points)):
         # If all target levels exist in the source cube, simply extract them.
@@ -727,8 +774,13 @@ def extract_levels(cube, levels, scheme, coordinate=None):
             raise ValueError(emsg.format(list(levels), name))
     else:
         # As a last resort, perform vertical interpolation.
-        result = _vertical_interpolate(cube, src_levels, levels, scheme,
-                                       extrap_scheme)
+        result = _vertical_interpolate(
+            cube,
+            src_levels,
+            levels,
+            interpolation,
+            extrapolation,
+        )
 
     return result
 

tests/integration/preprocessor/_regrid/test_extract_levels.py

@@ -45,6 +45,14 @@ def test_nop__levels_match(self):
         self.assertEqual(result, self.cube)
         self.assertEqual(id(result), id(self.cube))
 
+    def test_levels_almost_match(self):
+        vcoord = self.cube.coord(axis='z', dim_coords=True)
+        levels = np.array(vcoord.points, dtype=float)
+        vcoord.points = vcoord.points + 1.e-7
+        result = extract_levels(self.cube, levels, 'linear')
+        self.assert_array_equal(vcoord.points, levels)
+        self.assertTrue(result is self.cube)
+
     def test_interpolation__linear(self):
         levels = [0.5, 1.5]
         scheme = 'linear'

tests/unit/cmor/test_cmor_check.py

@@ -641,6 +641,20 @@ def test_non_requested(self):
         checker.check_metadata()
         self.assertTrue(checker.has_warnings())
 
+    def test_almost_requested(self):
+        """Automatically fix tiny deviations from the requested values."""
+        coord = self.cube.coord('air_pressure')
+        values = np.array(coord.points)
+        values[0] += 1.e-7
+        values[-1] += 1.e-7
+        self._update_coordinate_values(self.cube, coord, values)
+        checker = CMORCheck(self.cube, self.var_info, automatic_fixes=True)
+        checker.check_metadata()
+        reference = np.array(
+            self.var_info.coordinates['air_pressure'].requested, dtype=float)
+        np.testing.assert_array_equal(
+            self.cube.coord('air_pressure').points, reference)
+
     def test_requested_str_values(self):
         """
         Warning if requested values are not present.

tests/unit/preprocessor/_regrid/test_extract_levels.py

@@ -1,21 +1,25 @@
 """Unit tests for :func:`esmvalcore.preprocessor.regrid.extract_levels`."""
 import os
-
+import tempfile
 import unittest
 from unittest import mock
 
 import iris
 import numpy as np
 from numpy import ma
-import tempfile
 
 import tests
-from esmvalcore.preprocessor._regrid import (_MDI, VERTICAL_SCHEMES,
-                                             extract_levels)
+from esmvalcore.preprocessor._regrid import (
+    _MDI,
+    VERTICAL_SCHEMES,
+    extract_levels,
+    parse_vertical_scheme,
+)
 from tests.unit.preprocessor._regrid import _make_cube, _make_vcoord
 
 
 class Test(tests.Test):
+
     def setUp(self):
         self.shape = (3, 2, 1)
         self.z = self.shape[0]
@@ -44,6 +48,17 @@ def test_invalid_scheme__unknown(self):
     def test_vertical_schemes(self):
         self.assertEqual(set(VERTICAL_SCHEMES), set(self.schemes))
 
+    def test_parse_vertical_schemes(self):
+        reference = {
+            'linear': ('linear', 'nan'),
+            'nearest': ('nearest', 'nan'),
+            'linear_horizontal_extrapolate_vertical': ('linear', 'nearest'),
+            'nearest_horizontal_extrapolate_vertical': ('nearest', 'nearest'),
+        }
+        for scheme in self.schemes:
+            interpolation, extrapolation = parse_vertical_scheme(scheme)
+            assert interpolation, extrapolation == reference[scheme]
+
     def test_nop__levels_match(self):
         vcoord = _make_vcoord(self.z, dtype=self.dtype)
         self.assertEqual(self.cube.coord(axis='z', dim_coords=True), vcoord)