supplementary variables disappear from cubes after extract_region

[rebeccaherman1]

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 #2161

^{Originally posted by rebeccaherman1 August 1, 2023}
I am using this notebook 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?

[valeriupredoi]

hi @rebeccaherman1 could be able to provide me with a minimal recipe that I can use to replicate this behaviour please? 🍺

[rebeccaherman1]

@valeriupredoi Here's a notebook:
extract_region_supplementary_example.zip

Would you prefer to have it in a different format? I could potentially dump my cube into a pkl file and shorten the notebook.

Thank you!

[valeriupredoi]

cheers @rebeccaherman1 - could you pls give me the recipe you are running, or, if you are not running a recipe, a minimal reproducible Python code that I can run, sorry, Notebooks are too fancy for my 1990s ways of reproducing issues 😆

[rebeccaherman1]

@valeriupredoi no worries :)
I've been finding it easier to debug in the notebook; I'll turn it into a recipe when it works the way I want it to.
In the meantime, here's a python script version of it:
extract_region_supplementary.zip
(sorry I didn't see your response right away)

[rebeccaherman1]

A rather confusing (at least to me) update:
This behavior seems to be dependent on the region that I extract!

cell_area and land_area_fraction disappear when I use the region [-90,40]x[20,80], but remain when I use the region [160,280]x[-5,5]. I have no idea what to make of this.

This region-dependent behavior holds for at least two variables: I've tested psl (used in the python script) and ts.

[valeriupredoi]

no worries @rebeccaherman1 - I think I can make do with your Python script you attached above - I'll have a looksee at what/how it does things and I'll investigate the fairly bizarre cases you are seeing 👍

[valeriupredoi]

@rebeccaherman1 I have tried to reproduce this issue, and I used a pytest-like Python test script:

import os
import iris
import esmvalcore.preprocessor as pp

from esmvalcore.dataset import Dataset
from esmvalcore.config import CFG
from esmvalcore.esgf import download, ESGFFile


root_fil = "/home/valeriu/climate_data/CMIP6/CMIP/NCAR/CESM2/piControl/r1i1p1f1/Amon/psl/gn/v20190320"
CESM2 = Dataset(
    mip='Amon',
    project='CMIP6',
    dataset='CESM2',
    ensemble='*',
    exp='piControl',
    timerange='P501Y/*',
    grid='gn'
)
CFG['search_esgf'] = 'when_missing'

def add_facet(D, fname, fvalue):
    R = D.copy()
    R.facets[fname]=fvalue
    R.augment_facets()
    return R


def get_files(D):
    D.find_files()
    files = list(D.files)
    for supplementary_ds in D.supplementaries:
        files.extend(supplementary_ds.files)
    files = [file for file in files if isinstance(file, ESGFFile)]
    if(len(files)>0):
        print("Downloading Files:")
        print(files)
        download(files, CFG['download_dir'])
    D.find_files()
    return D


def bug_extract_region():
    """Test bug of reg extract."""
    NAO_xlim = [-90, 40]
    NAO_ylim = [20, 80]
    files = os.listdir(root_fil)

    # individual files
    for fil in files:
        print(fil)
        psl = iris.load_cube(os.path.join(root_fil, fil))
        print(f"Initial cell measures: {psl.cell_methods}")
        test = pp.extract_region(psl,
                                 start_longitude=NAO_xlim[0],
                                 end_longitude=NAO_xlim[1],
                                 start_latitude=NAO_ylim[0],
                                 end_latitude=NAO_ylim[1])
        print(f"Final cell measures: {test.cell_methods}")

    # cubelist concatenated into one cube
    cube_list = iris.cube.CubeList(
        [iris.load_cube(os.path.join(root_fil, fil)) for fil in files]
    )
    cubes = pp.concatenate(cube_list)
    test = pp.extract_region(cubes,
                             start_longitude=NAO_xlim[0],
                             end_longitude=NAO_xlim[1],
                             start_latitude=NAO_ylim[0],
                             end_latitude=NAO_ylim[1])
    print(test)
    print(f"Final cell measure concat cubes: {test.cell_methods}")

    # use Dataset object
    psl = get_files(add_facet(CESM2, 'short_name', 'psl'))
    ps = psl.load()
    ds_test = pp.extract_region(ps,
                                start_longitude=NAO_xlim[0],
                                end_longitude=NAO_xlim[1],
                                start_latitude=NAO_ylim[0],
                                end_latitude=NAO_ylim[1])
    print(ds_test)
    print(f"DATASET Final cell measures: {ds_test.cell_methods}")


if __name__ == '__main__':
    bug_extract_region()

You can run it yourself by simply changing the root_fil to point to where your data is stored (run it like any bogstandard Python script) - I don't see any loss of CellMethod() in any of the cubes, after looking at each of the files individually, nor after concatenating them into a single cube - all this using the pure iris load, nor after invoking the Dataset class and doing everything using it (see the third test case inside the function. Moreover I looked at the two cubes (the concatenated cube loaded by pure iris and the one obtained from filtering the data through Dataset class) and they are (almost) identical - almost here since Dataset fixes time axes, attributes etc. I have run this with esmvalcore=2.9.0. Could you please maybe grab this code, run it with the esmvalcore version I ran it with, and point to me what you think it's wrong/missing? I am attaching it as a txt file as well, for ease of download from here, just rename it 🍺
bug_extract_region.py.txt

[rebeccaherman1]

I think you actually removed the lines of code where I added the supplementary variables that are disappearing:

for s in ['areacella', 'sftlf']:
    CESM2.add_supplementary(short_name=s, mip='fx')

Then, though you mention cell measures in your print statements, it looks to me like you are examining cell_methods rather than cell_measures and ancillary_variables, which is where I believe the supplementary variable information is stored. To be honest, I'm not yet entirely sure what CellMethod() represents; it's the cell measures that I've been concerned about.

Could you try again to replicate my error after adding the supplementary variables, and while looking at cell_measures and ancillary_variables instead of cell_methods ?

[rebeccaherman1]

try this? and maybe repeat with ancillary_variables instead of cell_measures?

import os
import iris
import esmvalcore.preprocessor as pp

from esmvalcore.dataset import Dataset
from esmvalcore.config import CFG
from esmvalcore.esgf import download, ESGFFile


root_fil = "/home/valeriu/climate_data/CMIP6/CMIP/NCAR/CESM2/piControl/r1i1p1f1/Amon/psl/gn/v20190320"
CESM2 = Dataset(
    mip='Amon',
    project='CMIP6',
    dataset='CESM2',
    ensemble='*',
    exp='piControl',
    timerange='P501Y/*',
    grid='gn'
)
for s in ['areacella', 'sftlf']:
    CESM2.add_supplementary(short_name=s, mip='fx')     #ADDED CODE HERE

CFG['search_esgf'] = 'when_missing'

def add_facet(D, fname, fvalue):
    R = D.copy()
    R.facets[fname]=fvalue
    R.augment_facets()
    return R


def get_files(D):
    D.find_files()
    files = list(D.files)
    for supplementary_ds in D.supplementaries:
        files.extend(supplementary_ds.files)
    files = [file for file in files if isinstance(file, ESGFFile)]
    if(len(files)>0):
        print("Downloading Files:")
        print(files)
        download(files, CFG['download_dir'])
    D.find_files()
    return D


def bug_extract_region():
    """Test bug of reg extract."""
    NAO_xlim = [-90, 40]
    NAO_ylim = [20, 80]
    files = os.listdir(root_fil)

    # individual files
    for fil in files:
        print(fil)
        psl = iris.load_cube(os.path.join(root_fil, fil))
        print(f"Initial cell measures: {psl.cell_measures}")       #CHANGED 
        test = pp.extract_region(psl,
                                 start_longitude=NAO_xlim[0],
                                 end_longitude=NAO_xlim[1],
                                 start_latitude=NAO_ylim[0],
                                 end_latitude=NAO_ylim[1])
        print(f"Final cell measures: {test.cell_measures}")              #CHANGED

    # cubelist concatenated into one cube
    cube_list = iris.cube.CubeList(
        [iris.load_cube(os.path.join(root_fil, fil)) for fil in files]
    )
    cubes = pp.concatenate(cube_list)
    test = pp.extract_region(cubes,
                             start_longitude=NAO_xlim[0],
                             end_longitude=NAO_xlim[1],
                             start_latitude=NAO_ylim[0],
                             end_latitude=NAO_ylim[1])
    print(test)
    print(f"Final cell measure concat cubes: {test.cell_measures}")          #CHANGED

    # use Dataset object
    psl = get_files(add_facet(CESM2, 'short_name', 'psl'))
    ps = psl.load()
    ds_test = pp.extract_region(ps,
                                start_longitude=NAO_xlim[0],
                                end_longitude=NAO_xlim[1],
                                start_latitude=NAO_ylim[0],
                                end_latitude=NAO_ylim[1])
    print(ds_test)
    print(f"DATASET Final cell measures: {ds_test.cell_measures}")           #CHANGED


if __name__ == '__main__':
    bug_extract_region()

bug_extract_region_edited.py.txt

[valeriupredoi]

@rebeccaherman1 I think you are misinterpreting the functionality of add_supplementary() convenience function:

from docs:

        This is a convenience function that will create a copy of the current
        dataset, update its facets with the values specified in ``**facets``,
        and append it to :obj:`Dataset.supplementaries`. For more control
        over the creation of the supplementary dataset, first create a new
        :class:`Dataset` describing the supplementary dataset and then append
        it to :obj:`Dataset.supplementaries`.

I don't think there is a guarantee that the supplementary variables will persist through the preprocessing chain as if they were hard-added to the data cube; the way cell_measure variables are added as actual cube members is done via esmvalcore/preprocessor/_supplementary_vars.py and that invokes a specific iris function cube.add_cell_measure() that changes the cube in accordance and saves the information in the cube metadata. That's why I didn't reuse that bit of your code - in all fairness, I am not entirely sure why @bouweandela has created that convenience function, and not made it private

[valeriupredoi]

also, as an advice - please run your things using a recipe - that way you are guaranteed to make full use of the power of the preprocessor since the preprocessing chain is executed correctly and uses its internal functionalities in the correct manner 😃