Added new Extract surface preprocessor

doc/recipe/preprocessor.rst

@@ -1326,6 +1326,7 @@ The ``_volume.py`` module contains the following preprocessor functions:
 * ``extract_transect``: Extract data along a line of constant latitude or
   longitude.
 * ``extract_trajectory``: Extract data along a specified trajectory.
+* ``extract_surface``: Extract the surface layer from a cube.
 
 
 ``extract_volume``
@@ -1421,6 +1422,24 @@ Note that this function uses the expensive ``interpolate`` method from
 See also :func:`esmvalcore.preprocessor.extract_trajectory`.
 
 
+``extract_surface``
+-------------------
+
+This function extracts the surface layer from a dataset. 
+
+The surface layer is defined as the minimum
+of the absolute value of the Z-dimension. 
+This is typically the case for ocean models, but might not be the
+case for atmospheric models. 
+
+The same functionality exists in the extract_levels preprocessor,
+and this function should be used for more complex datasets
+if extract_surface fails. However, extract_levels also 
+has a higher computational cost, and may be slower. 
+
+See also :func:`esmvalcore.preprocessor.extract_levels`.
+
+
 .. _cycles:
 
 Cycles

esmvalcore/preprocessor/__init__.py

@@ -68,6 +68,7 @@
     extract_trajectory,
     extract_transect,
     extract_volume,
+    extract_surface,
     volume_statistics,
 )
 from ._weighting import weighting_landsea_fraction
@@ -99,6 +100,7 @@
     'resample_hours',
     'resample_time',
     # Level extraction
+    'extract_surface',
     'extract_levels',
     # Weighting
     'weighting_landsea_fraction',

esmvalcore/preprocessor/_volume.py

@@ -295,6 +295,68 @@ def volume_statistics(
     return _create_cube_time(src_cube, result, times)
 
 
+def extract_surface(cube):
+    """
+    Extact the surface layer from a 3D cube.
+
+    Extracts the entire layer along the z axis where the z coordinate is
+    the minimum of the absolute value of the z-axis dimensions points.
+
+    Requires a cube with a z axis.
+
+    This assumes that the surface is the closest layer to zero.
+    This assumption is usually true in the ocean, but may not be the
+    case in all models for all z-axes.
+
+    In the case of temporally or spatially varying depth grid,
+    (ie the depth array is 3D or 4D), the surface layer is determined
+    using the mean of the depth points along the time, latitude and
+    longitude axes. This preprocessor may also behave strangely
+    if the z axis data are not monotonic, or unusual in some other way.
+
+    The preprocessor extract_layer can also do this, but it may be
+    much slower and more memory intensie as it regrids.
+    It is more robust to variability in the grid definitions.
+
+    Arguments
+    ---------
+    cube: iris.cube.Cube
+        input cube.
+
+    Returns
+    -------
+    iris.cube.Cube
+        collapsed cube.
+
+    """
+    # Get the z axis.
+    zcoord = cube.coord(axis='Z')
+    positive = zcoord.attributes['positive']
+
+    # Get the Z points:
+    if zcoord.points.ndim == 1:
+        points = zcoord.points
+    elif zcoord.points.ndim == 3:
+        points = zcoord.points.mean(axis=[1, 2])
+    elif zcoord.points.ndim == 4:
+        points = zcoord.points.mean(axis=[0, 2, 3])
+
+    # Calculate the surface layer index:
+    surf = np.abs(points).argmin()
+
+    # Get the z axis dimension in the cude:
+    zcoord_dim = cube.coord_dims(zcoord)
+    if zcoord_dim in [0, (0,)]:
+        return cube[surf]
+
+    if zcoord_dim in [1, (1,)]:
+        return cube[:, surf]
+
+    logger.error('Condition not recognised: postive: %s , zcoord_dim: %s, '
+                 'surface level: %s', positive, zcoord_dim, surf)
+    return cube
+
+
 def depth_integration(cube):
     """
     Determine the total sum over the vertical component.

tests/unit/preprocessor/_volume/test_volume.py

@@ -10,7 +10,9 @@
 from esmvalcore.preprocessor._volume import (volume_statistics,
                                              depth_integration,
                                              extract_trajectory,
-                                             extract_transect, extract_volume)
+                                             extract_transect,
+                                             extract_volume,
+                                             extract_surface)
 
 
 class Test(tests.Test):
@@ -19,6 +21,7 @@ class Test(tests.Test):
     def setUp(self):
         """Prepare tests"""
         coord_sys = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
+        data0 = np.ones((2, 2))
         data1 = np.ones((3, 2, 2))
         data2 = np.ma.ones((2, 3, 2, 2))
         data3 = np.ma.ones((4, 3, 2, 2))
@@ -61,6 +64,9 @@ def setUp(self):
                                      units='degrees_north',
                                      coord_system=coord_sys)
 
+        coords_spec2 = [(lats2, 0), (lons2, 1)]
+        self.grid_2d = iris.cube.Cube(data0, dim_coords_and_dims=coords_spec2)
+
         coords_spec3 = [(zcoord, 0), (lats2, 1), (lons2, 2)]
         self.grid_3d = iris.cube.Cube(data1, dim_coords_and_dims=coords_spec3)
 
@@ -152,6 +158,13 @@ def test_extract_trajectory(self):
         expected = np.ones((3, 2))
         self.assert_array_equal(result.data, expected)
 
+    def test_extract_surface(self):
+        """Tests to extract the surface from a 3D and 4D cube."""
+        result = extract_surface(self.grid_3d)
+        self.assert_array_equal(result.data, self.grid_2d.data)
+        result2 = extract_surface(self.grid_4d)
+        self.assert_array_equal(self.grid_4d[:, 0, :, :].data, result2.data)
+
 
 if __name__ == '__main__':
     unittest.main()