Rolling cubes

[ghost]

When working with global fields it is occasionally useful to extract data across the longitude boundary, for example when attempting to extract the region 30W-30E. At the moment it is not straightforward to extract such a region in iris (for plotting or analysis), but Carwyn (cpelley) has provided a solution using the following code

import iris
import numpy as np

extb = iris.load_strict(fname)

longco = extb.coord('longitude')
if longco.is_monotonic():
    extb.data = np.roll(extb.data, len(longco.points)//2)
    longco.points -= 180.
    if longco.has_bounds():
        longco.bounds = longco.bounds - 180.
else:
    print 'ERROR: the points are not monotonic, this workaround will not work'
    sys.exit(0)

after which a suitably phrased constraint can be used to extract the data from the cube extb.

The numpy method to "roll" an array along an index could be directly implemented within the cube object allowing this task to be performed more easily. Such a method could also be very useful in preparing global fields to be plotted with the most important region centred.

[pelson]

Hey Matthew, welcome to github!

Thanks for reporting/suggesting this. Obviously the code you have provided has some significant limitations, but is clearly useful in some cases.

FYI There is already a discussion going on about bounding box extractions, for which this would make a good usecase (at #77).

I agree that there is a need for this functionality, especially with regards to analysis. It should be noted that I am currently working on a tool which should make visualisation of data much easier (on a map) in its native coordinate system, avoiding the need to transform the coordinate system of the underlying data (which the example code is doing).

All the best,

[bblay]

Welcome, @msmizielinski !

[cpelley]

I have at least two other users on previous occasions who have also reported a need for this so I think it is definitely a common requirement and a necessary enhancement

[rhattersley]

No one wants to modify a Cube like this as an end in itself - it's a means to an end. As @pelson suggests, there may be other changes which would remove the need to do this, so it's worth considering why this is being done so we can best address the real need.

[cpelley]

The underlying issue here is that a safe extraction is required irrespective of whether it over the meridian or not. This is what needs to be addressed. In addition to this, it would be useful to be able to modify the plotting region so that you don't have data on either side with nothing at the centre.

There is another user who has contacted Ian in regards to this issue. This is now 4+ users who require to extract a region over a longitudinal boundary and adjust the plot so that it is not split.

[rhattersley]

Sounds like the 4+ users might be best served by:

1. A way to extract a region, which ensures the extracted cube has a simple, contiguous domain. e.g. Source longitudes in range [-180 to 180] combined with an extraction longitude range of [170 to 190] would give result longitudes in the range [170 to 190].
2. The default map for a non-global cube would be centred on the cube's longitude range. e.g. A cube of [170 to 190] would give a map centred on 180.

[isedwards]

+1
(or perhaps +4!). A simple extract (specifying a region using the same coords as the data) with the above behaviour is exactly the functionality that is being requested. It is clear that the users do not just want to visualise the data, but instead require a contiguous dataset to be returned for analysis.

[pelson]

It is clear that the users do not just want to visualise the data, but instead require a contiguous dataset to be returned for analysis.

I don't disagree that this is desirable. But I do want to stress that the "contiguous dataset" is not a requirement to do a lot of analyses (masking would work equally well in a large number of cases. Where it wouldn't work might be in gradients and other "connectivity" related analyses).

[isedwards]

I believe I need the "contiguous dataset" solution for outputting the data array in a more primitive file format (geotiff #224). The data must be arranged sequentially from -180:180.

[bblay]

It's time to nail this down into a concrete plan.
The purpose of this post is to present a solution for critical analysis and/or mockery.

I suggest:

• Upgrade _iris.constraints.Constraint.extract():
• Look for a coord constraint on a circular coord where the constraint goes beyond the points range.
  • This can be done by taking the left and right adjacent coords (i.e with + and - mod) and checking if any points pass the test in those, but not the original.
  • Possibly ugly, but I think this works? Does floating point accuracy get in the way?
• If found, apply an updated version of Carwyn's code (above).
  • The update will be to the shift amount, which will apply an "appropriate shift", not just halve it.
  • An "appropriate shift" is the number of points from the first point to the first extraction point.

So. Is this solution acceptable? Is it flawed? Are there any nicer alternatves?

[bblay]

Moving to the discussion forum https://groups.google.com/forum/#!topic/scitools-iris-dev/puXL8oTjNlw

[pelson]

Updating @cpelley's function to work with iris v1.1/v1.2:

import numpy as np

import iris
import iris.tests.stock as stock


def roll_cube(cube):
    """Takes a cube which goes longitude 0-360 back to -180-180."""

    lon = cube.coord('longitude')

    cube.data = np.roll(cube.data, len(lon.points) // 2)
    lon.points -= 180.
    if lon.bounds is not None:
        lon.bounds -= 180.
    return cube



cube = stock.simple_pp()

print cube
print cube.coord('longitude').points.min()
print roll_cube(cube).coord('longitude').points.min()

[bblay]

Also see this topic for a comment on the axis keyword.

[bblay]

The default map for a non-global cube would be centred on the cube's longitude range. e.g. A cube of [170 to 190] would give a map centred on 180.

Also see this ticket for a related issue.

[rhattersley]

I'd like to submit a PR to address this issue so I've posted to the corresponding iris-dev discussion topic.