Repository: ESMValGroup/ESMValCore

Issue:
supplementary variables disappear from cubes after extract_region
After poking around a bit more, I believe that the `cell_measures` and all other supplementary variables actually are disappearing from the cubes upon applying `extract_region`, whereas I would have expected information from spatially-anchored supplementary variables simply to be extracted over the same region. This seems like a major shortcoming to me, not least because to use `area_statistics` (which requires `cell_area`) over any specific region, you need to use `extract_region` first. 

### Discussed in https://github.com/ESMValGroup/ESMValCore/discussions/2161

<div type='discussions-op-text'>

<sup>Originally posted by **rebeccaherman1** August  1, 2023</sup>
I am using [this notebook](https://docs.esmvaltool.org/projects/ESMValCore/en/latest/notebooks/loading-and-processing-data.html) as a guide while I attempt to develop my recipe and diagnostic scripts. Upon loading the initial cube, I believe I see that the supplementary variable `areacella` included in the cube as `cell_area` under `Cell measures`, but the moment I restrict the domain with `extract_region`, `Cell measures` disappears, taking `cell_area` with it. Is areacella no longer part of the cube? If it is, where did it go?</div>

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2166?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) Report
> Merging [#2166](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2166?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (2e8c8b4) into [main](https://app.codecov.io/gh/ESMValGroup/ESMValCore/commit/3222d28a8587e7ce3f5871522c035079b645dd5b?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) (3222d28) will **decrease** coverage by `0.09%`.
> The diff coverage is `48.00%`.

```diff
@@            Coverage Diff             @@
##             main    #2166      +/-   ##
==========================================
- Coverage   93.09%   93.00%   -0.09%     
==========================================
  Files         237      237              
  Lines       12795    12818      +23     
==========================================
+ Hits        11911    11921      +10     
- Misses        884      897      +13     
```


| [Files Changed](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2166?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup) | Coverage Δ | |
|---|---|---|
| [esmvalcore/preprocessor/\_area.py](https://app.codecov.io/gh/ESMValGroup/ESMValCore/pull/2166?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=ESMValGroup#diff-ZXNtdmFsY29yZS9wcmVwcm9jZXNzb3IvX2FyZWEucHk=) | `89.78% <48.00%> (-3.44%)` | :arrow_down: |


valeriupredoi: BTW @ledm I have a feeling that this bit you as well, in your continual battle with fx variables
valeriupredoi: with the latest changes when you try print the cube that's been changed you get this:
```
Traceback (most recent call last):
  File "/home/valeriu/test_intersect.py", line 27, in <module>
    print("ex1", ex1)
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/cube.py", line 2563, in __str__
    return self.summary()
           ^^^^^^^^^^^^^^
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/cube.py", line 2558, in summary
    printer = CubePrinter(self)
              ^^^^^^^^^^^^^^^^^
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/_representation/cube_printout.py", line 165, in __init__
    cube_summary = CubeSummary(cube_or_summary)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/_representation/cube_summary.py", line 364, in __init__
    add_vector_section(
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/_representation/cube_summary.py", line 355, in add_vector_section
    self.vector_sections[title] = VectorSection(
                                  ^^^^^^^^^^^^^^
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/_representation/cube_summary.py", line 203, in __init__
    self.contents = [
                    ^
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/_representation/cube_summary.py", line 204, in <listcomp>
    VectorSummary(cube, vector, iscoord) for vector in vectors
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/valeriu/miniconda3/envs/esmvaltool/lib/python3.11/site-packages/iris/_representation/cube_summary.py", line 141, in __init__
    dims = vector.cube_dims(cube)
           ^^^^^^^^^^^^^^^^
AttributeError: 'Cube' object has no attribute 'cube_dims'
```
valeriupredoi: the test is failing as well bc the metadata of the cube with ancil vars in preserves the metadata from the original cube `AssertionError: assert 'sftlf' == 'psl'`
valeriupredoi: am going back to the the latest working commit 93b085871e9a076156fd9a2fac16408f6b556e59
valeriupredoi: sorry @schlunma a few too many issues with your review suggestions (both at our end and iris - iris issue which I am not prepared to dive into at the moment, see https://github.com/ESMValGroup/ESMValCore/pull/2166#issuecomment-1669576532); if you think we'd benefit enough from those changes, I suggest we merge this, then you can open a PR based off `main` where you apply those changes, and test properly. I'll fix the test not to use the test data (and delete that file), then I'll revert the style changes resulted from `yapf`ing it, then I am ready to merge this :beer:  
valeriupredoi: OK Manu, I have finished poslishing this:
- used your suggestion to allow for irregular grids to be looked at
- fixed the test and removed the now unnecessary test data file
- re-blacked `_area.py` (well, reverted previously added style changes from yapf)

I think we should do this: merge this to have correct functionality for extract_region, then you and me we can go back to the suggestions you made earlier in a follow up PR, unfortunately there are quite a few hitches (not least the iris one that we should look at) to be looked at here; we may as well not even do it all if we'll remove all this when iris fix their intersection, but at least I'd like to understand that print(cube) issue above. What say you? :beer: 
valeriupredoi: Also, as a note to all of us (I am guilty as charged myself) - everytime we suggest major code changes in a PR, in the form of suggested code via the browser/GitHub API, it would be a lot easier for the PR owner if the reviewer tested those changes by either running the test provided in the PR or by running a small independent test - this, to catch serious issues (not style, minor issues with the suggestions, those can be taken care of the PR owner, after including the suggested code changes). I think we should add this to our best practices PR reviewing docs :+1: 
schlunma: Well, if you use `cube.intersection` (and not `extract_region`) to get the cell measure sub cube, this won't work for irregular grids. Hence you need to reinstate the check for `coord('latitude').ndim == 1`.

Also, in this form, this won't work (and will actually crash) for cell measures and ancillary variables that do not exactly span latitude/longitude. We need to take care of this properly here, otherwise it will break recipes.

I am honestly confused why my suggestion didn't work. From the failing test on CircleCI it looks just like a `flake8` [fail](https://app.circleci.com/pipelines/github/ESMValGroup/ESMValCore/9518/workflows/fedbbda2-6e1d-4b49-a08f-1ea8be52b428/jobs/41505)? I will give it a try myself in a local branch.
valeriupredoi: > Well, if you use `cube.intersection` (and not `extract_region`) to get the cell measure sub cube, this won't work for irregular grids. Hence you need to reinstate the check for `coord('latitude').ndim == 1`.
> 
> Also, in this form, this won't work (and will actually crash) for cell measures and ancillary variables that do not exactly span latitude/longitude. We need to take care of this properly here, otherwise it will break recipes.
> 
> I am honestly confused why my suggestion didn't work. From the failing test on CircleCI it looks just like a `flake8` [fail](https://app.circleci.com/pipelines/github/ESMValGroup/ESMValCore/9518/workflows/fedbbda2-6e1d-4b49-a08f-1ea8be52b428/jobs/41505)? I will give it a try myself in a local branch.

sounds like a very good plan! Cheers, Manu! Can grab this branch, edit locally, and when things look OK to you just push changes here - that was my point in https://github.com/ESMValGroup/ESMValCore/pull/2166#issuecomment-1669748583 - reviewer's minimal testing is something we should start doing :beer: 
schlunma: Tests pass fine for me locally (your test and my suggestion). Can I push it? 😁 
valeriupredoi: > Tests pass fine for me locally (your test and my suggestion). Can I push it? grin

Pliz do :beer: 
valeriupredoi: @schlunma let's get this in plz so I can help @rebeccaherman1 with an install at DKRZ :beer: 
valeriupredoi: > Thanks V! Can I be a bother and ask to to test this with the script mentioned [here](https://github.com/ESMValGroup/ESMValCore/issues/2162#issuecomment-1664161639)? I think we need to make sure that this properly works for cubes with time dimensions, too (If you like, you could also add a brief test of this grin).
> 
> I will be in a seminar today with limitied availability, so I pre-approve this. Thanks for the fix!!

Very many thanks, Manu! Yes, that's what I've been using to test the PR (alas, slightly changed, but I'll test with the exact script that Rebecca had initially). Good point, I'll put a time dim in the integration test in this PR :beer: 

Patch:
diff --git a/esmvalcore/preprocessor/_area.py b/esmvalcore/preprocessor/_area.py
index 621a5e47e..947781141 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
 
 


Test patch:
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 000000000..41442c9c0
--- /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