Repository: esmvalgroup/esmvalcore
Issue: extract_region preprocessor no longer works on irregular grids

Hi,

we recently changed the irregular grid section of the `extract_region` preprocessor in `preprocessor/_area.py` from this:

```
    # irregular grids
    lats = cube.coord('latitude').points
    lons = cube.coord('longitude').points
    mask = np.ma.array(cube.data).mask
    mask += np.ma.masked_where(lats < start_latitude, lats).mask
    mask += np.ma.masked_where(lats > end_latitude, lats).mask
    mask += np.ma.masked_where(lons > start_longitude, lons).mask
    mask += np.ma.masked_where(lons > end_longitude, lons).mask
    cube.data = np.ma.masked_where(mask, cube.data)
    return cube
```

to this:
```
    # irregular grids
    lats = cube.coord('latitude').points
    lons = cube.coord('longitude').points
    select_lats = start_latitude < lats < end_latitude
    select_lons = start_longitude < lons < end_longitude
    selection = select_lats & select_lons
    data = da.ma.masked_where(~selection, cube.core_data())
    return cube.copy(data)
```

The new code doesn't work for irregular grids. It produces the error message:
```
  File "/users/modellers/ledm/workspace/ESMValToolTest/ESMValTool/esmvalcore/preprocessor/_area.py", line 82, in extract_region
    select_lats = start_latitude < lats < end_latitude
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

```
Any suggestions on how to fix it and keep it lazy? For now, I've reverted my local code back to the previous version. (but I'm not going to push it!)

The reason that we didn't catch the problem at the time is because the ocean recipes that use this were designed to demonstrate the diagnostics, not test the preprocessor. I'm changing recipe_ocean_example.yml to include an irregular grid CMIP5 model for future testing.



Changed files:
esmvalcore/preprocessor/_area.py, tests/unit/preprocessor/_area/test_area.py

diff --git a/esmvalcore/preprocessor/_area.py b/esmvalcore/preprocessor/_area.py
index bc5ae1a4b..6db338ff5 100644
--- a/esmvalcore/preprocessor/_area.py
+++ b/esmvalcore/preprocessor/_area.py
@@ -46,12 +46,10 @@ def extract_region(cube, start_longitude, end_longitude, start_latitude,
     iris.cube.Cube
         smaller cube.
     """
-    # Converts Negative longitudes to 0 -> 360. standard
-    start_longitude = float(start_longitude)
-    end_longitude = float(end_longitude)
-    start_latitude = float(start_latitude)
-    end_latitude = float(end_latitude)
-
+    if abs(start_latitude) > 90.:
+        raise ValueError(f"Invalid start_latitude: {start_latitude}")
+    if abs(end_latitude) > 90.:
+        raise ValueError(f"Invalid end_latitude: {end_latitude}")
     if cube.coord('latitude').ndim == 1:
         # Iris check if any point of the cell is inside the region
         # To check only the center, ignore_bounds must be set to
@@ -63,14 +61,29 @@ def extract_region(cube, start_longitude, end_longitude, start_latitude,
         )
         region_subset = region_subset.intersection(longitude=(0., 360.))
         return region_subset
-    # irregular grids
+    # Irregular grids
     lats = cube.coord('latitude').points
     lons = cube.coord('longitude').points
-    select_lats = start_latitude <= lats <= end_latitude
-    select_lons = start_longitude <= lons <= end_longitude
+    # Convert longitudes to valid range
+    if start_longitude != 360.:
+        start_longitude %= 360.
+    if end_longitude != 360.:
+        end_longitude %= 360.
+
+    if start_longitude <= end_longitude:
+        select_lons = (lons >= start_longitude) & (lons <= end_longitude)
+    else:
+        select_lons = (lons >= start_longitude) | (lons <= end_longitude)
+
+    if start_latitude <= end_latitude:
+        select_lats = (lats >= start_latitude) & (lats <= end_latitude)
+    else:
+        select_lats = (lats >= start_latitude) | (lats <= end_latitude)
+
     selection = select_lats & select_lons
-    data = da.ma.masked_where(~selection, cube.core_data())
-    return cube.copy(data)
+    selection = da.broadcast_to(selection, cube.shape)
+    cube.data = da.ma.masked_where(~selection, cube.core_data())
+    return cube
 
 
 def zonal_means(cube, coordinate, mean_type):
diff --git a/tests/unit/preprocessor/_area/test_area.py b/tests/unit/preprocessor/_area/test_area.py
index bc4404465..09a4b0345 100644
--- a/tests/unit/preprocessor/_area/test_area.py
+++ b/tests/unit/preprocessor/_area/test_area.py
@@ -154,8 +154,139 @@ class Test(tests.Test):
             extract_named_regions(region_cube, ['region1', 'reg_A'])
 
 
+def create_irregular_grid_cube(data, lons, lats):
+    """Create test cube on irregular grid."""
+    nlat = iris.coords.DimCoord(range(data.shape[0]), var_name='nlat')
+    nlon = iris.coords.DimCoord(range(data.shape[1]), var_name='nlon')
+    lat = iris.coords.AuxCoord(lats,
+                               var_name='lat',
+                               standard_name='latitude',
+                               units='degrees')
+    lon = iris.coords.AuxCoord(lons,
+                               var_name='lon',
+                               standard_name='longitude',
+                               units='degrees')
+    dim_coord_spec = [
+        (nlat, 0),
+        (nlon, 1),
+    ]
+    aux_coord_spec = [
+        (lat, [0, 1]),
+        (lon, [0, 1]),
+    ]
+    cube = iris.cube.Cube(
+        data,
+        var_name='tos',
+        units='K',
+        dim_coords_and_dims=dim_coord_spec,
+        aux_coords_and_dims=aux_coord_spec,
+    )
+    return cube
+
+
+IRREGULAR_EXTRACT_REGION_TESTS = [
+    {
+        'region': (10, 360, 0, 90),
+        'mask':
+        np.array(
+            [
+                [True, True, True],
+                [True, True, False],
+                [True, False, False],
+            ],
+            dtype=bool,
+        ),
+    },
+    {
+        'region': (200, 10, -90, -60),
+        'mask':
+        np.array(
+            [
+                [False, True, False],
+                [True, True, True],
+                [True, True, True],
+            ],
+            dtype=bool,
+        ),
+    },
+    {
+        'region': (-150, 50, 50, -50),
+        'mask':
+        np.array(
+            [
+                [False, True, False],
+                [True, True, True],
+                [False, True, False],
+            ],
+            dtype=bool,
+        ),
+    },
+    {
+        'region': (0, 0, -100, 0),
+        'raises': "Invalid start_latitude: -100."
+    },
+    {
+        'region': (0, 0, 0, 100),
+        'raises': "Invalid end_latitude: -100."
+    },
+]
+
+
+@pytest.fixture
+def irregular_extract_region_cube():
+    """Create a test cube on an irregular grid to test `extract_region`."""
+    data = np.arange(9, dtype=np.float32).reshape((3, 3))
+    lons = np.array(
+        [
+            [0, 120, 240],
+            [0, 120, 240],
+            [0, 120, 240],
+        ],
+        dtype=np.float64,
+    )
+    lats = np.array(
+        [
+            [-60, -61., -60],
+            [0, -1, 0],
+            [60, 60, 60],
+        ],
+        dtype=np.float64,
+    )
+    cube = create_irregular_grid_cube(data, lons, lats)
+    return cube
+
+
+@pytest.mark.parametrize('case', IRREGULAR_EXTRACT_REGION_TESTS)
+def test_extract_region_irregular(irregular_extract_region_cube, case):
+    """Test `extract_region` with data on an irregular grid."""
+    start_lon, end_lon, start_lat, end_lat = case['region']
+    if 'raises' not in case:
+        cube = extract_region(
+            irregular_extract_region_cube,
+            start_longitude=start_lon,
+            end_longitude=end_lon,
+            start_latitude=start_lat,
+            end_latitude=end_lat,
+        )
+
+        data = np.arange(9, dtype=np.float32).reshape((3, 3))
+        np.testing.assert_array_equal(cube.data.mask, case['mask'])
+        np.testing.assert_array_equal(cube.data.data, data)
+    else:
+        with pytest.raises(ValueError) as exc:
+            extract_region(
+                irregular_extract_region_cube,
+                start_longitude=start_lon,
+                end_longitude=end_lon,
+                start_latitude=start_lat,
+                end_latitude=end_lat,
+            )
+            assert exc.value == case['raises']
+
+
 @pytest.fixture
 def make_testcube():
+    """Create a test cube on a Cartesian grid."""
     coord_sys = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
     data = np.ones((5, 5))
     lons = iris.coords.DimCoord(
@@ -174,6 +305,7 @@ def make_testcube():
 
 
 def write_shapefile(shape, path):
+    """Write a shape to a shapefile."""
     # Define a polygon feature geometry with one attribute
     schema = {
         'geometry': 'Polygon',
@@ -239,9 +371,8 @@ def test_extract_shape(make_testcube, square_shape, tmp_path, crop):
 
 
 @pytest.fixture
-def irregular_grid_cube():
-    """Create test cube on irregular grid."""
-    # Define grid and data
+def irreg_extract_shape_cube():
+    """Create a test cube on an irregular grid to test `extract_shape`."""
     data = np.arange(9, dtype=np.float32).reshape((3, 3))
     lats = np.array(
         [
@@ -259,38 +390,12 @@ def irregular_grid_cube():
         ],
         dtype=np.float64,
     )
-
-    # Construct cube
-    nlat = iris.coords.DimCoord(range(data.shape[0]), var_name='nlat')
-    nlon = iris.coords.DimCoord(range(data.shape[1]), var_name='nlon')
-    lat = iris.coords.AuxCoord(lats,
-                               var_name='lat',
-                               standard_name='latitude',
-                               units='degrees')
-    lon = iris.coords.AuxCoord(lons,
-                               var_name='lon',
-                               standard_name='longitude',
-                               units='degrees')
-    dim_coord_spec = [
-        (nlat, 0),
-        (nlon, 1),
-    ]
-    aux_coord_spec = [
-        (lat, [0, 1]),
-        (lon, [0, 1]),
-    ]
-    cube = iris.cube.Cube(
-        data,
-        var_name='tos',
-        units='K',
-        dim_coords_and_dims=dim_coord_spec,
-        aux_coords_and_dims=aux_coord_spec,
-    )
+    cube = create_irregular_grid_cube(data, lons, lats)
     return cube
 
 
 @pytest.mark.parametrize('method', ['contains', 'representative'])
-def test_extract_shape_irregular(irregular_grid_cube, tmp_path, method):
+def test_extract_shape_irregular(irreg_extract_shape_cube, tmp_path, method):
     """Test `extract_shape` with a cube on an irregular grid."""
     # Points are (lon, lat)
     shape = Polygon([
@@ -303,7 +408,7 @@ def test_extract_shape_irregular(irregular_grid_cube, tmp_path, method):
     shapefile = tmp_path / 'shapefile.shp'
     write_shapefile(shape, shapefile)
 
-    cube = extract_shape(irregular_grid_cube, shapefile, method)
+    cube = extract_shape(irreg_extract_shape_cube, shapefile, method)
 
     data = np.arange(9, dtype=np.float32).reshape((3, 3))
     mask = np.array(