diff --git a/esmvalcore/preprocessor/_area.py b/esmvalcore/preprocessor/_area.py index 621a5e47e1..9477811415 100644 --- a/esmvalcore/preprocessor/_area.py +++ b/esmvalcore/preprocessor/_area.py @@ -65,6 +65,10 @@ def extract_region(cube, start_longitude, end_longitude, start_latitude, iris.cube.Cube smaller cube. """ + # first examine if any cell_measures are present + cell_measures = cube.cell_measures() + ancil_vars = cube.ancillary_variables() + if abs(start_latitude) > 90.: raise ValueError(f"Invalid start_latitude: {start_latitude}") if abs(end_latitude) > 90.: @@ -87,6 +91,56 @@ def extract_region(cube, start_longitude, end_longitude, start_latitude, start_latitude, end_latitude, ) + + # put back cell measures and ancillary_variables; + # iris.Cube.cube.intersection removes them both. + # This is a workaround resulting from opening upstream + # https://github.com/SciTools/iris/issues/5413 + # When removing this block after iris have a fix, make sure to remove the + # test too tests/integration/preprocessor/_extract_region/ + + def _extract_region_from_dim_metadata(dim_metadata, dim_metadata_dims): + """Extract region from dimensional metadata.""" + idx = tuple(( + slice(None) if d in dim_metadata_dims else 0 + for d in range(cube.ndim) + )) + subcube = cube[idx].copy(dim_metadata.core_data()) + for sub_cm in subcube.cell_measures(): + subcube.remove_cell_measure(sub_cm) + for sub_av in subcube.ancillary_variables(): + subcube.remove_ancillary_variable(sub_av) + subcube = extract_region( + subcube, + start_longitude, + end_longitude, + start_latitude, + end_latitude, + ) + return dim_metadata.copy(subcube.core_data()) + + # Step 1: cell measures + if cell_measures and not region_subset.cell_measures(): + for cell_measure in cell_measures: + cell_measure_dims = cube.cell_measure_dims(cell_measure) + cell_measure_subset = _extract_region_from_dim_metadata( + cell_measure, cell_measure_dims + ) + region_subset.add_cell_measure( + cell_measure_subset, cell_measure_dims + ) + + # Step 2: ancillary variables + if ancil_vars and not region_subset.ancillary_variables(): + for ancil_var in ancil_vars: + ancil_var_dims = cube.ancillary_variable_dims(ancil_var) + ancil_var_subset = _extract_region_from_dim_metadata( + ancil_var, ancil_var_dims + ) + region_subset.add_ancillary_variable( + ancil_var_subset, ancil_var_dims + ) + return region_subset diff --git a/tests/integration/preprocessor/_extract_region/test_intersect.py b/tests/integration/preprocessor/_extract_region/test_intersect.py new file mode 100644 index 0000000000..41442c9c01 --- /dev/null +++ b/tests/integration/preprocessor/_extract_region/test_intersect.py @@ -0,0 +1,98 @@ +""" +Temporary test. + +Remove this test and test file after iris fixes this +https://github.com/SciTools/iris/issues/5413 . +""" +import iris +import numpy as np +from cf_units import Unit +from iris.fileformats.pp import EARTH_RADIUS + +from esmvalcore.preprocessor._area import extract_region +from esmvalcore.preprocessor._shared import guess_bounds + + +def make_cube(): + """Make a realistic cube with cell measure and ancil var.""" + coord_sys = iris.coord_systems.GeogCS(EARTH_RADIUS) + data = np.ones((10, 192, 288), dtype=np.float32) + time = iris.coords.DimCoord( + np.arange(0., 10., 1.), + standard_name='time', + units=Unit('days since 1950-01-01 00:00:00', + calendar='360_day')) + lons = iris.coords.DimCoord( + [i + .5 for i in range(288)], + standard_name='longitude', + bounds=[[i, i + 1.] for i in range(288)], + units='degrees_east', + coord_system=coord_sys) + lats = iris.coords.DimCoord( + [i + .5 for i in range(192)], + standard_name='latitude', + bounds=[[i, i + 1.] for i in range(192)], + units='degrees_north', + coord_system=coord_sys, + ) + coords_spec = [(time, 0), (lats, 1), (lons, 2)] + simple_cube = iris.cube.Cube(data, dim_coords_and_dims=coords_spec) + + # add a cell measure + simple_cube = guess_bounds(simple_cube, ['longitude', 'latitude']) + grid_areas = iris.analysis.cartography.area_weights(simple_cube) + measure = iris.coords.CellMeasure( + grid_areas, + standard_name='cell_area', + units='m2', + measure='area') + simple_cube.add_cell_measure(measure, range(0, measure.ndim)) + + # add ancillary variable + ancillary_var = iris.coords.AncillaryVariable( + simple_cube.data, + standard_name='land_ice_area_fraction', + var_name='sftgif', + units='%') + simple_cube.add_ancillary_variable(ancillary_var, + range(0, simple_cube.ndim)) + + return simple_cube + + +def test_extract_region_cell_ancil(): + """Test re-adding cell measures ancil variables after extract region.""" + cube = make_cube() + + # intersection cube loses cellmeas/ancillary variables + # under normal (unpatched) conditions of extract_region + ex1 = extract_region(cube, + start_longitude=-90, + end_longitude=40, + start_latitude=20, + end_latitude=80) + + # intersection cube doesn't lose cellmeas/ancillary variables + # under normal (unpatched) conditions of extract_region + # so duplication must be avoided + ex2 = extract_region(cube, + start_longitude=160, + end_longitude=280, + start_latitude=-5, + end_latitude=5) + + expected_cm = cube.cell_measures()[0] + result_cm = ex1.cell_measures() + assert result_cm + assert expected_cm.measure == result_cm[0].measure + assert expected_cm.var_name == result_cm[0].var_name + np.testing.assert_array_equal(result_cm[0].shape, (10, 60, 58)) + assert expected_cm.standard_name == result_cm[0].standard_name + expected_ancil = cube.ancillary_variables()[0] + result_ancil = ex1.ancillary_variables() + assert result_ancil + assert expected_ancil.var_name == result_ancil[0].var_name + assert expected_ancil.standard_name == result_ancil[0].standard_name + np.testing.assert_array_equal(result_ancil[0].shape, (10, 60, 58)) + assert len(ex2.cell_measures()) == 1 + assert len(ex2.ancillary_variables()) == 1