Cube coord bounds make cube intersection misbehave

[lukasbrunner]

Issue in a nutshell

• operation: _area.py/extract_region()
• task at hand: func needs to perfrom cube intersection on 0-360 lons
• comment: extract_region works well when run on a cube with negative longitude points and run as a stand alone preprocessor step, see comment but the intersection operation misbehaves when a regrid step has been run before, see comment

Original Post

Hi all!

I am currently testing the new recipe ClimWIP, which, among other things, calls the pre-processor to extract the following region:
https://github.com/ESMValGroup/ESMValTool/blob/414d7d2e84a62a777dfada79297c0ed400b6bbd3/esmvaltool/recipes/recipe_climwip.yml#L52-L56

For me the resulting longitude values look like this:
lon = -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75, 16.25, 18.75, 21.25, 23.75, 26.25, 28.75, 31.25, 33.75, 36.25, 38.75, 351.25, 353.75, 356.25 ;

Can anyone recreate this behaviour? Is it even expected? To me this looks wrong. My understanding is that the longitude should either be in [-180, 180) OR in [0, 360) but neither case is true here.

Many thanks,
Lukas

main_log_debug.txt

[valeriupredoi]

Hello @lukasbrunner cheers for raising this! This should indeed not happen since the longitude coordinate is CMOR standardized and has 0-360 bounds, in fact, the preprocessor function extract_region subsets the region on 0-360 by:

region_subset = region_subset.intersection(longitude=(0., 360.))

I will have a look and try reproduce the issue. In the meantime I suggest @Peter9192 and/or @stefsmeets have a look at the recipe since you guys are recipe maintainers 👍 🍺

[valeriupredoi]

I ran the extract_region preprocessor on (most of) the data in that recipe (a few files are missing off the CEDA ESGF node on JASMIN), and I can not reproduce the behavior of negative longitude points; could you please tell me which dataset gives you those negative points by running this little script, on data that has been preprocessed only with the extract_region preprocessor:

import iris
import glob

files = glob.glob("$YOUR_PREPROC_DIR/tas/*nc")
for fn in files:
    c = iris.load_cube(fn)
    long_points = c.coord("longitude").points
    if long_points[0] < 0.:
        print(fn)
        print(long_points[0])

where $YOUR_PREPROC_DIR is your preproc dir, of course. Cheers 🍺

[lukasbrunner]

Hey @valeriupredoi!

Thanks for looking into this! Attached is the output from your script. Its wired that you don't get that behaviour, I had a college check on her setup and she got the same issue (though it is the same data - the ETH archive). Is there anything else I can provide or check? Not sure how much of a problem that can become for us.

esmval_check_lon.log

[Peter9192]

perhaps you can run with
esmvaltool run ... ... --save_intermediary_cubes=True
Then you can see exactly after which step these errors are introduced. But it produces a lot of data, so you may want to remove them manually afterwards.

[Peter9192]

I managed to reproduce the problem, and indeed it first appears after the extract region step.

[Peter9192]

Further pinned down to exactly this line you mentioned @valeriupredoi

region_subset = region_subset.intersection(longitude=(0., 360.))

this doesn't seem to work as expected.

[valeriupredoi]

good stuff @Peter9192 - could you please tell me what datasets produce this behavior and also tell me what iris version you using? 🍺

[Peter9192]

{dataset: ACCESS1-0, project: CMIP5, exp: [historical, rcp85], ensemble: r1i1p1, start_year: 1995, end_year: 2014, mip: Amon, short_name: tas}
iris 2.4.0

[valeriupredoi]

cool, cheers for the info @Peter9192 - there is something rather fishy going on here, I am printing out the lon points righte before and right after the cube intersection on the [0, 360] interval:

• before: [-9.375, -7.5 , -5.625, -3.75 , -1.875, 0. , 1.875, 3.75 ,..., 37.5]
• after: [ 0. , 1.875, 3.75... 37.5 , 350.625, 352.5 , 354.375, 356.25 , 358.125]

So as you can see the selection works and rotates the negative values on the 0-360 base; can you do the same for yours please (print the coord right before and after l.64 in _area.py from esmvalcore/master). Cheers 🍺

[Peter9192]

Interesting! My values are different...

Before: [-8.75, -6.25, -3.75, -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, ..., 36.25, 38.75]
After: [ -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75, 16.25, ..., 38.75, 351.25, 353.75, 356.25]

[Peter9192]

do you have bounds? So for me the extracted region (the outer bounds) is -10 to 40 rather than -10 to 39.

[valeriupredoi]

OK I found the cause - this thing gets spookier the closer you look - the points I have in the comment are from running just the extract_region prproc step in the preprocessor; now, if you add the regrid step right before that I can reproduce the problem ad literam (before/after that cube intersection)

• before: DimCoord(array([-8.75, -6.25, -3.75, -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75, 16.25, 18.75, 21.25, 23.75, 26.25, 28.75, 31.25, 33.75, 36.25, 38.75])
• after: DimCoord(array([ -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75, 16.25, 18.75, 21.25, 23.75, 26.25, 28.75, 31.25, 33.75, 36.25, 38.75, 351.25, 353.75, 356.25])

so it looks like the intersection works fine, but then it stops working fine after regrid - I have to investigate more. I am, however, transferring this issue to esmvalcore since the problem is there 🍺 Oh and I did change from -10 to 40 here, but this is not the problem (I changed it and ran it without regrid and it worked well, no neg points)

[Peter9192]

Good! At least we get the same results now, that helps.

[valeriupredoi]

it's the bounds that fool intersection, here's a minimal example:

import iris
import numpy as np

lon_points = np.array([-8.75, -6.25, -3.75, -1.25,  1.25,  3.75,  6.25,  8.75, 11.25,
                       13.75, 16.25, 18.75, 21.25, 23.75, 26.25, 28.75, 31.25, 33.75,
                       36.25, 38.75])
lon_bounds = np.array([
    [ -2.5,   0. ],
    [  0. ,   2.5],
    [  2.5,   5. ],
    [  5. ,   7.5],
    [  7.5,  10. ],
    [ 10. ,  12.5],
    [ 12.5,  15. ],
    [ 15. ,  17.5],
    [ 17.5,  20. ],
    [ 20. ,  22.5],
    [ 22.5,  25. ],
    [ 25. ,  27.5],
    [ 27.5,  30. ],
    [ 30. ,  32.5],
    [ 32.5,  35. ],
    [ 35. ,  37.5],
    [ 37.5,  40. ],
    [350. , 352.5],
    [352.5, 355. ],
    [355. , 357.5]])
lons = iris.coords.DimCoord(lon_points, standard_name='longitude', bounds=lon_bounds,
                            units='degrees_east')
cube = iris.cube.Cube(np.zeros((20,)), dim_coords_and_dims=[(lons, 0)])
isect_with_bounds = cube.intersection(longitude=(0, 360))
print(isect_with_bounds.coord("longitude"))
cube.coord("longitude").bounds = None
isect_without_bounds = cube.intersection(longitude=(0, 360))
print(isect_without_bounds.coord("longitude"))

[Peter9192]

Seems to be this issue then: SciTools/iris#3391

[valeriupredoi]

Indeed, good find @Peter9192 - I posted the minimal code so the iris folk can see the behavior 👍

[Peter9192]

Okay let's wait for their reply then. It would be better to solve this in iris then for us to hack around it.

[valeriupredoi]

too late - I already created a temporary fix for us 😁

[valeriupredoi]

this has been fixed in iris=3.0.2 and tested by meself - you can view the original Iris SciTools issue at SciTools/iris#3391 and the PR that fixed it at SciTools/iris#4059 (note that the unwanted behavior still persists for iris=3.0.1)