Fixed clip_timerange if only a single time point is extracted

esmvalcore/preprocessor/_time.py

@@ -157,6 +157,21 @@ def _duration_to_date(duration, reference, sign):
     return date
 
 
+def _restore_time_coord_position(cube, original_time_index):
+    """Restore original ordering of coordinates."""
+    # Coordinates before time
+    new_order = list(np.arange(original_time_index) + 1)
+
+    # Time coordinate
+    new_order.append(0)
+
+    # Coordinates after time
+    new_order = new_order + list(range(original_time_index + 1, cube.ndim))
+
+    # Transpose cube in-place
+    cube.transpose(new_order)
+
+
 def _extract_datetime(cube, start_datetime, end_datetime):
     """Extract a time range from a cube.
 
@@ -213,6 +228,17 @@ def _extract_datetime(cube, start_datetime, end_datetime):
             f"to {end_datetime.strftime('%Y-%m-%d')} is outside "
             f"cube time bounds {time_coord.cell(0)} to {time_coord.cell(-1)}.")
 
+    # If only a single point in time is extracted, the new time coordinate of
+    # cube_slice is a scalar coordinate. Convert this back to a regular
+    # dimensional coordinate with length 1. Note that iris.util.new_axis always
+    # puts the new axis at index 0, so we need to reorder the coordinates in
+    # case the original time coordinate was not at index 0.
+    if cube_slice.ndim < cube.ndim:
+        cube_slice = iris.util.new_axis(cube_slice, 'time')
+        original_time_index = cube.coord_dims(time_coord)[0]
+        if original_time_index != 0:
+            _restore_time_coord_position(cube_slice, original_time_index)
+
     return cube_slice
 
 

tests/unit/preprocessor/_time/test_time.py

@@ -325,6 +325,119 @@ def test_clip_timerange_30_day(self):
         assert_array_equal(
                 sliced_cube.coord('time').points, expected_time)
 
+    def test_clip_timerange_single_year_1d(self):
+        """Test that single year stays dimensional coordinate."""
+        cube = self._create_cube([0.0], [150.0], [[0.0, 365.0]], 'standard')
+        sliced_cube = clip_timerange(cube, '1950/1950')
+
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert_array_equal(sliced_cube.coord('time').bounds, [[0.0, 365.0]])
+        assert cube.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+
+        # Repeat test without bounds
+        cube.coord('time').bounds = None
+        sliced_cube = clip_timerange(cube, '1950/1950')
+
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert sliced_cube.coord('time').bounds is None
+        assert cube.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+
+    def test_clip_timerange_single_year_2d(self):
+        """Test that single year stays dimensional coordinate."""
+        cube = self._create_cube([[0.0, 1.0]], [150.0], [[0.0, 365.0]],
+                                 'standard')
+        lat_coord = iris.coords.DimCoord([10.0, 20.0],
+                                         standard_name='latitude')
+        cube.add_dim_coord(lat_coord, 1)
+        sliced_cube = clip_timerange(cube, '1950/1950')
+
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert_array_equal(sliced_cube.coord('time').bounds, [[0.0, 365.0]])
+        assert cube.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+
+        # Repeat test without bounds
+        cube.coord('time').bounds = None
+        sliced_cube = clip_timerange(cube, '1950/1950')
+
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert sliced_cube.coord('time').bounds is None
+        assert cube.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+
+    def test_clip_timerange_single_year_4d(self):
+        """Test time is not scalar even when time is not first coordinate."""
+        cube = self._create_cube([[[[0.0, 1.0]]]], [150.0], [[0.0, 365.0]],
+                                 'standard')
+        plev_coord = iris.coords.DimCoord([1013.0],
+                                          standard_name='air_pressure')
+        lat_coord = iris.coords.DimCoord([10.0], standard_name='latitude')
+        lon_coord = iris.coords.DimCoord([0.0, 1.0], standard_name='longitude')
+        cube.add_dim_coord(plev_coord, 1)
+        cube.add_dim_coord(lat_coord, 2)
+        cube.add_dim_coord(lon_coord, 3)
+
+        # Order: plev, time, lat, lon
+        cube_1 = cube.copy()
+        cube_1.transpose([1, 0, 2, 3])
+        assert cube_1.shape == (1, 1, 1, 2)
+        sliced_cube = clip_timerange(cube_1, '1950/1950')
+
+        assert sliced_cube is not cube_1
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert_array_equal(sliced_cube.coord('time').bounds, [[0.0, 365.0]])
+        assert cube_1.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+        for coord_name in [c.name() for c in cube_1.coords()]:
+            assert (sliced_cube.coord_dims(coord_name) ==
+                    cube_1.coord_dims(coord_name))
+
+        # Order: lat, lon, time, plev
+        cube_2 = cube.copy()
+        cube_2.transpose([2, 3, 0, 1])
+        assert cube_2.shape == (1, 2, 1, 1)
+        sliced_cube = clip_timerange(cube_2, '1950/1950')
+
+        assert sliced_cube is not cube_2
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert_array_equal(sliced_cube.coord('time').bounds, [[0.0, 365.0]])
+        assert cube_2.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+        for coord_name in [c.name() for c in cube_2.coords()]:
+            assert (sliced_cube.coord_dims(coord_name) ==
+                    cube_2.coord_dims(coord_name))
+
+        # Order: lon, lat, plev, time
+        cube_3 = cube.copy()
+        cube_3.transpose([3, 2, 1, 0])
+        assert cube_3.shape == (2, 1, 1, 1)
+        sliced_cube = clip_timerange(cube_3, '1950/1950')
+
+        assert sliced_cube is not cube_3
+        assert sliced_cube.coord('time').units == Unit(
+            'days since 1950-01-01', calendar='standard')
+        assert_array_equal(sliced_cube.coord('time').points, [150.0])
+        assert_array_equal(sliced_cube.coord('time').bounds, [[0.0, 365.0]])
+        assert cube_3.shape == sliced_cube.shape
+        assert sliced_cube.coord('time', dim_coords=True)
+        for coord_name in [c.name() for c in cube_3.coords()]:
+            assert (sliced_cube.coord_dims(coord_name) ==
+                    cube_3.coord_dims(coord_name))
+
 
 class TestExtractSeason(tests.Test):
     """Tests for extract_season."""