IRIS longitude constraints...

[Nick Dunstone]

Hi,

I'm just starting the transition from the comfortable world of IDL to the brave new one of Python/IRIS. This is my first post to this group.

Sometimes different observation datasets (or even models) define the longitude coordinate differently, sometimes 0->360, other times -180->180.
In IDL if I want to extract a region, the plotting code would be clever enough to sort this out behind the scenes, for example:
when I ask for longitudes -70 to 35 [70E to 35W] it would know that is equivalent to [270 to 0 combined with 0 to 35], if the file uses 0->360.

In Python/IRIS, I need to do this:
            if obs_name == 'hadisst':
                loncon = iris.Constraint(longitude=lambda v: 290 < v < 360 or 0 < v < 35)
            else:
                loncon = iris.Constraint(longitude=lambda v: -70 < v < 35)
which seems a wee bit clumsy to me. Is there a special method for selecting longitude ranges that I'm unaware of, or is this additional functionality that could be added?

Thanks,

Nick

[Richard Hattersley]

Hi Nick,

Thanks for giving Iris a try and for posting your feedback/question. I completely agree that having to use different Constraints is clumsy, and quite few other people do besides! So in the development version of Iris (which will be published as Iris v1.7 in the next few weeks) we've added a Cube method which takes care of this for you:

region_cube = global_cube.intersection(longitude=(-70, 35))

The result will always be a cube with longitudes in the range -70 to 35, irrespective of the original longitudes. NB. At the moment there can still be some problems using this method to select other coordinate ranges, so for latitude it's best to stick with constraints for now.

If you're using Iris within the Met Office you should be able to use the "testing" release which will give you early access to this method. If you have any problems with that then please just contact AVD and they'll be very happy to get you up and running.

Regards,
Richard

[Nick Dunstone]

Hi Richard,

Thanks for your reply. It's great to see that a solution has been implemented.

This seems to work, but when plotted it still includes the whole (global) longitude range:
iplt.pcolormesh(danom_regrid[i].intersection(longitude=(-70, 35),latitude=(25,75)), vmin=vmin_model,vmax=vmax_model,cmap=brewer_cmap)

I've included an attached image to show this.

Is there a way around this? Else the "intersection" fix is of limited use I guess.

Cheers,
Nick

[Andrew Dawson]

You can manually change the limits of the plot, which should suffice for the short term:

 

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
...

iplt.pcolormesh(danom_regrid[i].intersection(longitude=(-70, 35),latitude=(25,75)), vmin=vmin_model,vmax=vmax_model,cmap=brewer_cmap)

ax = plt.gca()
ax.set_extent([-70, 35, 25, 75], crs=ccrs.PlateCarree())
...

[Richard Hattersley]

That looks like a bug.

In the result of the intersection, the longitude coordinate is being erroneously marked as "circular". This means the last column is regarded as connected to the first column. Hence the weird stretched out lines on your plot and its global longitude extent.

As a workaround for now I would suggest:

region = danom_regrid[i].intersection(longitude=(-70, 35),latitude=(25,75))
region.coord('longitude').circular = False

Then when plotted it will default to the chosen area.

[Richard Hattersley]

I've raised a fix as https://github.com/SciTools/iris/pull/1169

[Nick Dunstone]

Thanks both Richard and Andrew for solutions to this.

I tried yours Richard, it sounds like the eventual solution, but removing the "circular" attribute still produced the same plot (with the global longitude domain) for some reason.

Andrew, your resetting of the axis worked well for me - a good short term fix.

Thanks again for the help - I've learned lots.

Nick

[Richard Hattersley]

On Tuesday, 10 June 2014 13:40:44 UTC+1, Nick Dunstone wrote:

... removing the "circular" attribute still produced the same plot (with the global longitude domain) for some reason.

Have you cleared the circular flag before doing the plot? For example:

region = danom_regrid[i].intersection(longitude=(-70, 35),latitude=(25,75))
region.coord('longitude').circular = False


iplt.pcolormesh(region, vmin=vmin_model,vmax=vmax_model,cmap=brewer_cmap)

Richard

[Nick Dunstone]

Apologies Richard, that does indeed work! I'm not sure what I did wrong on my earlier test - probably what you suggested!

I like your solution, look forward to it being implemented in the .intersection routine itself.

Thanks,
Nick

[Jamie Kettleborough]

Hello,

I've just hit a similar problem and came across this post which covered the same issue.  Can I ask for some clarification?  Will constraints in longitude ever support periodic coordinates in the way expected - or is the advice to move to cube.intersection?  (If that's the case - will Constraints be deprecated at some stage?)

Sorry if this has been discussed somewhere else and I've missed it.

Jamie

On Tuesday, 10 June 2014 11:07:54 UTC+1, Richard Hattersley wrote:

[Nick Dunstone]

Hello,

I've also just hit another related problem.

I found this using two different cubes now...

I'm trying to use .intersection like so:
obs_reg_tmp = obs.intersection(longitude=[-170.,180.],latitude=[0.,90.])
*** ValueError: The bounds array must be strictly monotonic.

Through trial-and-error I've identified the "longitude" coordinate as the problem.
However, to my mind the bounds array for the coordinate is perfectly "monotonic":
print obs.coord('longitude')
DimCoord(array([-177.5, -172.5, -167.5, -162.5, -157.5, -152.5, -147.5, -142.5,
       -137.5, -132.5, -127.5, -122.5, -117.5, -112.5, -107.5, -102.5,
        -97.5,  -92.5,  -87.5,  -82.5,  -77.5,  -72.5,  -67.5,  -62.5,
        -57.5,  -52.5,  -47.5,  -42.5,  -37.5,  -32.5,  -27.5,  -22.5,
        -17.5,  -12.5,   -7.5,   -2.5,    2.5,    7.5,   12.5,   17.5,
         22.5,   27.5,   32.5,   37.5,   42.5,   47.5,   52.5,   57.5,
         62.5,   67.5,   72.5,   77.5,   82.5,   87.5,   92.5,   97.5,
        102.5,  107.5,  112.5,  117.5,  122.5,  127.5,  132.5,  137.5,
        142.5,  147.5,  152.5,  157.5,  162.5,  167.5,  172.5,  177.5], dtype=float32), bounds=array([[-180., -175.],
       [-175., -170.],
       [-170., -165.],
       [-165., -160.],
       [-160., -155.],
       [-155., -150.],
       [-150., -145.],
       [-145., -140.],
       [-140., -135.],
       [-135., -130.],
       [-130., -125.],
       [-125., -120.],
       [-120., -115.],
       [-115., -110.],
       [-110., -105.],
       [-105., -100.],
       [-100.,  -95.],
       [ -95.,  -90.],
       [ -90.,  -85.],
       [ -85.,  -80.],
       [ -80.,  -75.],
       [ -75.,  -70.],
       [ -70.,  -65.],
       [ -65.,  -60.],
       [ -60.,  -55.],
       [ -55.,  -50.],
       [ -50.,  -45.],
       [ -45.,  -40.],
       [ -40.,  -35.],
       [ -35.,  -30.],
       [ -30.,  -25.],
       [ -25.,  -20.],
       [ -20.,  -15.],
       [ -15.,  -10.],
       [ -10.,   -5.],
       [  -5.,    0.],
       [   0.,    5.],
       [   5.,   10.],
       [  10.,   15.],
       [  15.,   20.],
       [  20.,   25.],
       [  25.,   30.],
       [  30.,   35.],
       [  35.,   40.],
       [  40.,   45.],
       [  45.,   50.],
       [  50.,   55.],
       [  55.,   60.],
       [  60.,   65.],
       [  65.,   70.],
       [  70.,   75.],
       [  75.,   80.],
       [  80.,   85.],
       [  85.,   90.],
       [  90.,   95.],
       [  95.,  100.],
       [ 100.,  105.],
       [ 105.,  110.],
       [ 110.,  115.],
       [ 115.,  120.],
       [ 120.,  125.],
       [ 125.,  130.],
       [ 130.,  135.],
       [ 135.,  140.],
       [ 140.,  145.],
       [ 145.,  150.],
       [ 150.,  155.],
       [ 155.,  160.],
       [ 160.,  165.],
       [ 165.,  170.],
       [ 170.,  175.],
       [ 175.,  180.]], dtype=float32), standard_name='longitude', units=Unit('degrees'), coord_system=GeogCS(6371229.0), circular=True)

I hope that someone can tell me what I'm doing wrong, or whether this is a bug in the .intersection code?

Thanks,
Nick

[Niall Robinson]

Is it anything to do with this special case bug?
https://github.com/SciTools/iris/issues/1173

if so, then I think its fixed in the incoming v1.7.0

[Niall Robinson]

Actually its different - but looks as though it may be part of the same problem - things failing the monotonicity test that shouldn't