From 8c0a86e407be5913554427afb35bbef59e89eddb Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Mon, 14 Jun 2021 11:11:27 +0200 Subject: [PATCH 1/6] Automatically fix small deviations in vertical levels --- esmvalcore/cmor/check.py | 7 ++++++- esmvalcore/preprocessor/_regrid.py | 15 ++++++++++++--- .../preprocessor/_regrid/test_extract_levels.py | 8 ++++++++ tests/unit/cmor/test_cmor_check.py | 16 +++++++++++++++- 4 files changed, 41 insertions(+), 5 deletions(-) diff --git a/esmvalcore/cmor/check.py b/esmvalcore/cmor/check.py index 48a5113d35..317d179930 100644 --- a/esmvalcore/cmor/check.py +++ b/esmvalcore/cmor/check.py @@ -693,7 +693,12 @@ 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( + [float(val) for val in coord_info.requested]) + atol = 1e-7 * np.mean(coord.points) + if self.automatic_fixes and np.allclose( + coord.points, cmor_points, rtol=1e-7, atol=atol): + coord.points = cmor_points except ValueError: cmor_points = coord_info.requested coord_points = list(coord.points) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index b35c4b3a25..5604331364 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -652,7 +652,9 @@ def extract_levels(cube, levels, scheme, coordinate=None): levels : array 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', @@ -669,7 +671,9 @@ def extract_levels(cube, levels, scheme, coordinate=None): Returns ------- - cube + iris.cube.Cube + A cube with the requested vertical levels. + See Also -------- @@ -710,10 +714,15 @@ def extract_levels(cube, levels, scheme, coordinate=None): src_levels = cube.coord(axis='z', dim_coords=True) if (src_levels.shape == levels.shape - and np.allclose(src_levels.points, levels)): + and np.allclose(src_levels.points, + levels, + rtol=1e-7, + atol=1e-7 * np.mean(src_levels.points))): # 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. 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 fa1c314166..c1f53702c1 100644 --- a/tests/unit/cmor/test_cmor_check.py +++ b/tests/unit/cmor/test_cmor_check.py @@ -9,7 +9,7 @@ import numpy as np from cf_units import Unit -from esmvalcore.cmor.check import CMORCheck, CMORCheckError, CheckLevels +from esmvalcore.cmor.check import CheckLevels, CMORCheck, CMORCheckError class VariableInfoMock: @@ -589,6 +589,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. From 8b90f897f653facaf3da41334892985962e1e2ec Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Fri, 1 Oct 2021 12:36:58 +0200 Subject: [PATCH 2/6] Make rtol and atol configurable --- esmvalcore/preprocessor/_regrid.py | 32 ++++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index eb9b783c87..2eee6e1ed1 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -642,14 +642,19 @@ 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 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. If the requested levels are sufficiently close to the @@ -668,6 +673,16 @@ 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 ------- @@ -713,11 +728,12 @@ 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, - rtol=1e-7, - atol=1e-7 * np.mean(src_levels.points))): + 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 From e5e0ee1f4a53d25bf0a1413457ad8882aaf1d7a6 Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Fri, 1 Oct 2021 12:37:25 +0200 Subject: [PATCH 3/6] Refactor automatic fix --- esmvalcore/cmor/check.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/esmvalcore/cmor/check.py b/esmvalcore/cmor/check.py index ded30afe06..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,17 +771,22 @@ def _check_requested_values(self, coord, coord_info, var_name): """Check requested values.""" if coord_info.requested: try: - cmor_points = np.array( - [float(val) for val in coord_info.requested]) - atol = 1e-7 * np.mean(coord.points) - if self.automatic_fixes and np.allclose( - coord.points, cmor_points, rtol=1e-7, atol=atol): - coord.points = cmor_points + 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)) From aa3e6d1a8ff42a7d0de0d2a99a66707044135dc3 Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Fri, 1 Oct 2021 15:07:02 +0200 Subject: [PATCH 4/6] Reduce complexity of extract_levels --- esmvalcore/preprocessor/_regrid.py | 59 ++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 16 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 2eee6e1ed1..c4072d81af 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -642,6 +642,41 @@ def _vertical_interpolate(cube, src_levels, levels, interpolation, return _create_cube(cube, new_data, src_levels, levels.astype(float)) +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, @@ -694,20 +729,7 @@ def extract_levels(cube, -------- 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) @@ -752,8 +774,13 @@ def extract_levels(cube, 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 From 348c9946dd96321c2dd73b4ce8026e3d09a242f6 Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Fri, 1 Oct 2021 15:52:28 +0200 Subject: [PATCH 5/6] Add unit test --- .../_regrid/test_extract_levels.py | 23 +++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) 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) From 076243a42cc56c76202160173df77a268bd14c4a Mon Sep 17 00:00:00 2001 From: Bouwe Andela Date: Mon, 4 Oct 2021 16:55:40 +0200 Subject: [PATCH 6/6] Add documentation --- doc/recipe/preprocessor.rst | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) 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`.