Repository: esmvalgroup/esmvalcore
Issue: Cubes lose derived coordinates during concatenation

**Describe the bug**
Due to a bug in `iris` (SciTools/iris#2478) , cubes with a derived coordinate (`AuxFactory`) lose this during the concatenation step of our preprocessor.

I will provide a PR for this.

Changed files:
esmvalcore/_recipe.py, tests/integration/preprocessor/_io/test_concatenate.py, esmvalcore/preprocessor/_io.py

diff --git a/esmvalcore/_recipe.py b/esmvalcore/_recipe.py
index ac8cbcbee..ac68be8c4 100644
--- a/esmvalcore/_recipe.py
+++ b/esmvalcore/_recipe.py
@@ -283,7 +283,7 @@ def _get_default_settings(variable, config_user, derive=False):
     settings['load'] = {
         'callback': concatenate_callback,
     }
-    # Configure merge
+    # Configure concatenation
     settings['concatenate'] = {}
 
     # Configure fixes
diff --git a/esmvalcore/preprocessor/_io.py b/esmvalcore/preprocessor/_io.py
index a43824c48..51fa0c201 100644
--- a/esmvalcore/preprocessor/_io.py
+++ b/esmvalcore/preprocessor/_io.py
@@ -13,6 +13,7 @@ import numpy as np
 import yaml
 
 from .._task import write_ncl_settings
+from ..cmor._fixes.shared import AtmosphereSigmaFactory
 from ._time import extract_time
 
 logger = logging.getLogger(__name__)
@@ -28,6 +29,55 @@ VARIABLE_KEYS = {
 }
 
 
+def _fix_aux_factories(cube):
+    """Fix :class:`iris.aux_factory.AuxCoordFactory` after concatenation.
+
+    Necessary because of bug in :mod:`iris` (see issue #2478).
+
+    """
+    coord_names = [coord.name() for coord in cube.coords()]
+
+    # Hybrid sigma pressure coordinate
+    # TODO possibly add support for other hybrid coordinates
+    if 'atmosphere_hybrid_sigma_pressure_coordinate' in coord_names:
+        new_aux_factory = iris.aux_factory.HybridPressureFactory(
+            delta=cube.coord(var_name='ap'),
+            sigma=cube.coord(var_name='b'),
+            surface_air_pressure=cube.coord(var_name='ps'),
+        )
+        for aux_factory in cube.aux_factories:
+            if isinstance(aux_factory, iris.aux_factory.HybridPressureFactory):
+                break
+        else:
+            cube.add_aux_factory(new_aux_factory)
+
+    # Hybrid sigma height coordinate
+    if 'atmosphere_hybrid_height_coordinate' in coord_names:
+        new_aux_factory = iris.aux_factory.HybridHeightFactory(
+            delta=cube.coord(var_name='lev'),
+            sigma=cube.coord(var_name='b'),
+            orography=cube.coord(var_name='orog'),
+        )
+        for aux_factory in cube.aux_factories:
+            if isinstance(aux_factory, iris.aux_factory.HybridHeightFactory):
+                break
+        else:
+            cube.add_aux_factory(new_aux_factory)
+
+    # Atmosphere sigma coordinate
+    if 'atmosphere_sigma_coordinate' in coord_names:
+        new_aux_factory = AtmosphereSigmaFactory(
+            pressure_at_top=cube.coord(var_name='ptop'),
+            sigma=cube.coord(var_name='lev'),
+            surface_air_pressure=cube.coord(var_name='ps'),
+        )
+        for aux_factory in cube.aux_factories:
+            if isinstance(aux_factory, AtmosphereSigmaFactory):
+                break
+        else:
+            cube.add_aux_factory(new_aux_factory)
+
+
 def _get_attr_from_field_coord(ncfield, coord_name, attr):
     if coord_name is not None:
         attrs = ncfield.cf_group[coord_name].cf_attrs()
@@ -107,9 +157,16 @@ def concatenate(cubes):
             concatenated = _concatenate_overlapping_cubes(concatenated)
 
     if len(concatenated) == 1:
-        return concatenated[0]
-
-    logger.error('Can not concatenate cubes into a single one.')
+        cube = concatenated[0]
+        _fix_aux_factories(cube)
+        return cube
+
+    # Concatenation not successful -> retrieve exact error message
+    try:
+        iris.cube.CubeList(cubes).concatenate_cube()
+    except iris.exceptions.ConcatenateError as exc:
+        msg = str(exc)
+    logger.error('Can not concatenate cubes into a single one: %s', msg)
     logger.error('Resulting cubes:')
     for cube in concatenated:
         logger.error(cube)
@@ -119,7 +176,7 @@ def concatenate(cubes):
             pass
         else:
             logger.error('From %s to %s', time.cell(0), time.cell(-1))
-    raise ValueError('Can not concatenate cubes.')
+    raise ValueError(f'Can not concatenate cubes: {msg}')
 
 
 def save(cubes, filename, optimize_access='', compress=False, **kwargs):
diff --git a/tests/integration/preprocessor/_io/test_concatenate.py b/tests/integration/preprocessor/_io/test_concatenate.py
index 2ee24bae0..1fbace4eb 100644
--- a/tests/integration/preprocessor/_io/test_concatenate.py
+++ b/tests/integration/preprocessor/_io/test_concatenate.py
@@ -1,18 +1,223 @@
 """Integration tests for :func:`esmvalcore.preprocessor._io.concatenate`."""
 
+import warnings
 import unittest
+from unittest.mock import call
 
 import numpy as np
+import pytest
 from cf_units import Unit
-from iris.coords import DimCoord
-from iris.cube import Cube
+from iris.aux_factory import HybridHeightFactory, HybridPressureFactory
+from iris.coords import AuxCoord, DimCoord
+from iris.cube import Cube, CubeList
 from iris.exceptions import ConcatenateError
 
 from esmvalcore.preprocessor import _io
 
 
+def get_hybrid_pressure_cube():
+    """Return cube with hybrid pressure coordinate."""
+    ap_coord = AuxCoord([1.0], bounds=[[0.0, 2.0]], var_name='ap', units='Pa')
+    b_coord = AuxCoord([0.0], bounds=[[-0.5, 1.5]], var_name='b')
+    ps_coord = AuxCoord([[[100000]]], var_name='ps', units='Pa')
+    x_coord = AuxCoord(
+        0.0,
+        var_name='x',
+        standard_name='atmosphere_hybrid_sigma_pressure_coordinate',
+    )
+    cube = Cube([[[[0.0]]]], var_name='x',
+                aux_coords_and_dims=[(ap_coord, 1), (b_coord, 1),
+                                     (ps_coord, (0, 2, 3)), (x_coord, ())])
+    return cube
+
+
+def get_hybrid_pressure_cube_list():
+    """Return list of cubes including hybrid pressure coordinate."""
+    cube_0 = get_hybrid_pressure_cube()
+    cube_1 = get_hybrid_pressure_cube()
+    cube_0.add_dim_coord(get_time_coord(0), 0)
+    cube_1.add_dim_coord(get_time_coord(1), 0)
+    cubes = CubeList([cube_0, cube_1])
+    for cube in cubes:
+        aux_factory = HybridPressureFactory(
+            delta=cube.coord(var_name='ap'),
+            sigma=cube.coord(var_name='b'),
+            surface_air_pressure=cube.coord(var_name='ps'),
+        )
+        cube.add_aux_factory(aux_factory)
+    return cubes
+
+
+def get_time_coord(time_point):
+    """Time coordinate."""
+    return DimCoord([time_point], var_name='time', standard_name='time',
+                    units='days since 6453-2-1')
+
+
+@pytest.fixture
+def mock_empty_cube():
+    """Return mocked cube with irrelevant coordinates."""
+    cube = unittest.mock.create_autospec(Cube, spec_set=True, instance=True)
+    a_coord = AuxCoord(0.0, var_name='a')
+    b_coord = AuxCoord(0.0, var_name='b')
+    cube.coords.return_value = [a_coord, b_coord]
+    return cube
+
+
+@pytest.fixture
+def mock_hybrid_height_cube():
+    """Return mocked cube with hybrid height coordinate."""
+    cube = unittest.mock.create_autospec(Cube, spec_set=True, instance=True)
+    lev_coord = AuxCoord([1.0], bounds=[[0.0, 2.0]], var_name='lev', units='m')
+    b_coord = AuxCoord([0.0], bounds=[[-0.5, 1.5]], var_name='b')
+    orog_coord = AuxCoord([[[100000]]], var_name='orog', units='m')
+    cube.coord.side_effect = [lev_coord, b_coord, orog_coord,
+                              lev_coord, b_coord, orog_coord]
+    cube.coords.return_value = [
+        lev_coord,
+        b_coord,
+        orog_coord,
+        AuxCoord(0.0, standard_name='atmosphere_hybrid_height_coordinate'),
+    ]
+    aux_factory = HybridHeightFactory(
+        delta=lev_coord,
+        sigma=b_coord,
+        orography=orog_coord,
+    )
+    cube.aux_factories = ['dummy', aux_factory]
+    return cube
+
+
+@pytest.fixture
+def mock_hybrid_pressure_cube():
+    """Return mocked cube with hybrid pressure coordinate."""
+    cube = unittest.mock.create_autospec(Cube, spec_set=True, instance=True)
+    ap_coord = AuxCoord([1.0], bounds=[[0.0, 2.0]], var_name='ap', units='Pa')
+    b_coord = AuxCoord([0.0], bounds=[[-0.5, 1.5]], var_name='b')
+    ps_coord = AuxCoord([[[100000]]], var_name='ps', units='Pa')
+    cube.coord.side_effect = [ap_coord, b_coord, ps_coord,
+                              ap_coord, b_coord, ps_coord]
+    cube.coords.return_value = [
+        ap_coord,
+        b_coord,
+        ps_coord,
+        AuxCoord(0.0,
+                 standard_name='atmosphere_hybrid_sigma_pressure_coordinate'),
+    ]
+    aux_factory = HybridPressureFactory(
+        delta=ap_coord,
+        sigma=b_coord,
+        surface_air_pressure=ps_coord,
+    )
+    cube.aux_factories = ['dummy', aux_factory]
+    return cube
+
+
+@pytest.fixture
+def real_hybrid_pressure_cube():
+    """Return real cube with hybrid pressure coordinate."""
+    return get_hybrid_pressure_cube()
+
+
+@pytest.fixture
+def real_hybrid_pressure_cube_list():
+    """Return real list of cubes with hybrid pressure coordinate."""
+    return get_hybrid_pressure_cube_list()
+
+
+def check_if_fix_aux_factories_is_necessary():
+    """Check if _fix_aux_factories() is necessary (i.e. iris bug is fixed)."""
+    cubes = get_hybrid_pressure_cube_list()
+    cube = cubes.concatenate_cube()
+    coords = [coord.name() for coord in cube.coords()]
+    msg = ("Apparently concatenation of cubes that have a derived variable "
+           "is now possible in iris (i.e. issue #2478 has been fixed). Thus, "
+           "this test and ALL appearances of the function "
+           "'_fix_aux_factories' can safely be removed!")
+    if 'air_pressure' in coords:
+        warnings.warn(msg)
+
+
+def test_fix_aux_factories_empty_cube(mock_empty_cube):
+    """Test fixing with empty cube."""
+    check_if_fix_aux_factories_is_necessary()
+    _io._fix_aux_factories(mock_empty_cube)
+    assert mock_empty_cube.mock_calls == [call.coords()]
+
+
+def test_fix_aux_factories_hybrid_height(mock_hybrid_height_cube):
+    """Test fixing of hybrid height coordinate."""
+    check_if_fix_aux_factories_is_necessary()
+
+    # Test with aux_factory object
+    _io._fix_aux_factories(mock_hybrid_height_cube)
+    mock_hybrid_height_cube.coords.assert_called_once_with()
+    mock_hybrid_height_cube.coord.assert_has_calls([call(var_name='lev'),
+                                                    call(var_name='b'),
+                                                    call(var_name='orog')])
+    mock_hybrid_height_cube.add_aux_factory.assert_not_called()
+
+    # Test without aux_factory object
+    mock_hybrid_height_cube.reset_mock()
+    mock_hybrid_height_cube.aux_factories = ['dummy']
+    _io._fix_aux_factories(mock_hybrid_height_cube)
+    mock_hybrid_height_cube.coords.assert_called_once_with()
+    mock_hybrid_height_cube.coord.assert_has_calls([call(var_name='lev'),
+                                                    call(var_name='b'),
+                                                    call(var_name='orog')])
+    mock_hybrid_height_cube.add_aux_factory.assert_called_once()
+
+
+def test_fix_aux_factories_hybrid_pressure(mock_hybrid_pressure_cube):
+    """Test fixing of hybrid pressure coordinate."""
+    check_if_fix_aux_factories_is_necessary()
+
+    # Test with aux_factory object
+    _io._fix_aux_factories(mock_hybrid_pressure_cube)
+    mock_hybrid_pressure_cube.coords.assert_called_once_with()
+    mock_hybrid_pressure_cube.coord.assert_has_calls([call(var_name='ap'),
+                                                      call(var_name='b'),
+                                                      call(var_name='ps')])
+    mock_hybrid_pressure_cube.add_aux_factory.assert_not_called()
+
+    # Test without aux_factory object
+    mock_hybrid_pressure_cube.reset_mock()
+    mock_hybrid_pressure_cube.aux_factories = ['dummy']
+    _io._fix_aux_factories(mock_hybrid_pressure_cube)
+    mock_hybrid_pressure_cube.coords.assert_called_once_with()
+    mock_hybrid_pressure_cube.coord.assert_has_calls([call(var_name='ap'),
+                                                      call(var_name='b'),
+                                                      call(var_name='ps')])
+    mock_hybrid_pressure_cube.add_aux_factory.assert_called_once()
+
+
+def test_fix_aux_factories_real_cube(real_hybrid_pressure_cube):
+    """Test fixing of hybrid pressure coordinate on real cube."""
+    check_if_fix_aux_factories_is_necessary()
+    assert not real_hybrid_pressure_cube.coords('air_pressure')
+    _io._fix_aux_factories(real_hybrid_pressure_cube)
+    air_pressure_coord = real_hybrid_pressure_cube.coord('air_pressure')
+    expected_coord = AuxCoord([[[[1.0]]]], bounds=[[[[[-50000., 150002.]]]]],
+                              standard_name='air_pressure', units='Pa')
+    assert air_pressure_coord == expected_coord
+
+
+def test_concatenation_with_aux_factory(real_hybrid_pressure_cube_list):
+    """Test actual concatenation of a cube with a derived coordinate."""
+    concatenated = _io.concatenate(real_hybrid_pressure_cube_list)
+    air_pressure_coord = concatenated.coord('air_pressure')
+    expected_coord = AuxCoord(
+        [[[[1.0]]], [[[1.0]]]],
+        bounds=[[[[[-50000.0, 150002.0]]]], [[[[-50000.0, 150002.0]]]]],
+        standard_name='air_pressure',
+        units='Pa',
+    )
+    assert air_pressure_coord == expected_coord
+
+
 class TestConcatenate(unittest.TestCase):
     """Tests for :func:`esmvalcore.preprocessor._io.concatenate`."""
+
     def setUp(self):
         """Start tests."""
         self._model_coord = DimCoord([1., 2.],
@@ -37,7 +242,7 @@ class TestConcatenate(unittest.TestCase):
             concatenated.coord('time').points, np.array([1, 2, 3, 4, 5, 6]))
 
     def test_concatenate_with_overlap(self):
-        """Test concatenation of time overalapping cubes"""
+        """Test concatenation of time overalapping cubes."""
         self._add_cube([6.5, 7.5], [6., 7.])
         concatenated = _io.concatenate(self.raw_cubes)
         np.testing.assert_array_equal(