diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index 566a4cd3f0..d0579913c2 100644 --- a/doc/recipe/preprocessor.rst +++ b/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`. diff --git a/esmvalcore/cmor/check.py b/esmvalcore/cmor/check.py index 957f50282f..2ba287f238 100644 --- a/esmvalcore/cmor/check.py +++ b/esmvalcore/cmor/check.py @@ -511,7 +511,6 @@ def _check_alternative_dim_names(self, key): search for the correct coordinate version and then check against this. For ``cmor_strict=False`` project (like OBS) the check for requested values might be disabled. - """ table_type = self._cmor_var.table_type alternative_coord = None @@ -772,12 +771,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)) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 7e09570d6a..c4072d81af 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/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 diff --git a/tests/integration/preprocessor/_regrid/test_extract_levels.py b/tests/integration/preprocessor/_regrid/test_extract_levels.py index 0267913509..267c067566 100644 --- a/tests/integration/preprocessor/_regrid/test_extract_levels.py +++ b/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' diff --git a/tests/unit/cmor/test_cmor_check.py b/tests/unit/cmor/test_cmor_check.py index 4dabcab7f7..f528e1271d 100644 --- a/tests/unit/cmor/test_cmor_check.py +++ b/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. diff --git a/tests/unit/preprocessor/_regrid/test_extract_levels.py b/tests/unit/preprocessor/_regrid/test_extract_levels.py index 965948bb54..ad422b725e 100644 --- a/tests/unit/preprocessor/_regrid/test_extract_levels.py +++ b/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)