BUG: Broken cartopy img_transform masks

[cpelley]

Putting this PR up here as I have not managed to finish the fixes for the bugs I found before Christmas.

Relates to changes made in #370
(specifically commit 02e764c)

Two bugs:
1. Mask array operations which are unsuitable for the case when dat.mask = False
TypeError: 'numpy.bool_' object does not support item assignment

Should instead be indexing the array itself and assigning a np.ma.masked value to these indexes.
i.e. dat[index] = np.ma.masked to mask these indexes.

2. The rectangle specified is unsuitable for the general case of determining which target points fall outside the source domain.

outside_source_domain = ((target_in_source_x < source_x_coords.min()) |
    (target_in_source_x > source_x_coords.max()) |
    (target_in_source_y < source_y_coords.min()) |
    (target_in_source_y > source_y_coords.max()))

This issue reveals itself when running one of the iris tests which has data on the range [0, 360] rather than [-180, 180].

These bugs currently don't allow iris to work with cartopy master.

[ajdawson]

1. Shouldn't we be creating a mask instead of setting data elements, or does that do the same thing?
2. Can we just transform the source points into the expected domain first?

What is broken in Iris due to this bug?

[cpelley]

Apologies for the delay, looking back at this now.

Shouldn't we be creating a mask instead of setting data elements, or does that do the same thing?

If assigning np.ma.masked then this is the same thing.

Can we just transform the source points into the expected domain first?

Looking into creating a shapely geometry in the source CS and transforming that into the resulting CS.

What is broken in Iris due to this bug?

Any data on the [0, 360] lon range:

...
filename = iris.tests.get_data_path(['PP', 'aPPglob1', 'global.pp'])
cube = iris.load_cube(filename)[::15, ::15]

target_proj = ccrs.PlateCarree()

fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=target_proj)
# Transform cube to target projection
new_cube, extent = iris.analysis.cartography.project(cube, target_proj,
                                                     nx=10, ny=10)
plt.pcolor(new_cube.coord('projection_x_coordinate').points,
           new_cube.coord('projection_y_coordinate').points,
           new_cube.data)
iplt.pcolor(new_cube)
ax.coastlines()

[ajdawson]

Any data on the [0, 360] lon range:

I did not realise you were talking about the iris.analysis.cartography.project function, I didn't know this used the cartopy.img_transform code.

[pelson]

Hmmm. Given we'd like to get a solid release out of cartopy, I'm minded to simply comment out the change that was added until we have a fix. @ajdawson how do you feel about that?

Cheers,

[ajdawson]

Not great but if it has to be done then so be it. You'll need to reopen the relevant bug(s?) and I'll have to retract my comments on StackOverflow saying the image projection bug is fixed in v0.10.0.

[pelson]

@cpelley - can you provide a SSCCE reproducing the problem - I'd like to investigate some ideas before moving forward with disabling the functionality.

Thanks.

[rhattersley]

What is broken in Iris due to this bug?

Any data on the [0, 360] lon range:

How about making a rectangle geometry from the source coordinates, and then passing it through cartopy's projection pipeline (to normalise the coordinate range and deal with splitting into two pieces). The result should be one or two rectangles whose extents can be plugged into the same sort of vectorised bounds checking as before.

[cpelley]

@pelson here you go. I welcome any ideas you might have, as I'm afraid I have exhausted what I think I can do within time constraints.

import iris.analysis.cartography
import cartopy.crs as ccrs
import iris
import iris.plot as iplt
import iris.tests as tests
import matplotlib.pyplot as plt
plt.switch_backend('tkagg')

filename = tests.get_data_path(['PP', 'aPPglob1', 'global.pp'])

#filename = '/data/local/cpelley/git/iris-test-data/test_data/PP/aPPglob1/global.pp'
cube = iris.load_cube(filename)[::25, ::25]

#target_proj = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)
target_proj = ccrs.PlateCarree()

# Set up figure
fig = plt.figure(figsize=(10, 10))

ax = plt.axes(projection=target_proj)
# Transform cube to target projection
new_cube, extent = iris.analysis.cartography.project(cube, target_proj,
                                                     nx=10, ny=10)
iplt.pcolor(new_cube)
# Add coastlines
ax.coastlines()
plt.show()

[cpelley]

I have tried two ways. In both, I to define a rectangle in the source coordinates (shapely geometry) and project that to the target projection. To determine determine which points fall inside/outside this shape is where the problem lies. There is no way using shapely to do this except on a point-by-point call. I attempted to do this with matplotlib instead using paths (contains_points method of the Path class). There are a few bugs with this however and after using matplotlib master (many of these bugs appear fixed) I still stumble at the final hurdle, see the following issue:

Here is the issue that got me stuck with matplotlib (possibly a bug):
matplotlib/matplotlib#2708

[cpelley]

How about making a rectangle geometry from the source coordinates, and then passing it through cartopy's projection pipeline (to normalise the coordinate range and deal with splitting into two pieces). The result should be one or two rectangles whose extents can be plugged into the same sort of vectorised bounds checking as before.

Thanks for your comments @rhattersley, can you clarify if what you suggest is different to what I have already done?

[rhattersley]

I to define a rectangle in the source coordinates (shapely geometry) and project that to the target projection.

I was suggesting projecting the source geometry into the source "projection". But in general that specific method won't necessarily exist because the source coordinate system might not be a projection.

To determine determine which points fall inside/outside this shape is where the problem lies. There is no way using shapely to do this except on a point-by-point call.

We definitely don't want to be doing point-by-point geometry comparisons!

To keep it quick and efficient we have to be doing comparisons in the source coordinate system as that is the only coordinate system where the bounds are rectangular. So all that's really needed is a way to normalise the source coordinate range to form one or more rectangles. The target_in_source_x and target_in_source_y values can then be quickly and efficiently checked against those rectangles.

[cpelley]

The result should be one or two rectangles whose extents can be plugged into the same sort of vectorised bounds checking as before.

So I transform a rectangle to the target projection and obtain bounds in the source projection (which could be two rectangles due to wrap around I'm guessing? Then check that points lie within these limits of set of bounds.
Oh I see, will that work for the general case?

[cpelley]

ignore my last comment (I posted before the page refreshed)

[pelson]

@pelson here you go

Thank you. I was hoping to remove Iris from the loop, but don't think you've got to that point yet so I will do that.

Cheers,

[cpelley]

@pelson I have re-added the unittest I wrote before (doesn't use cartopy).
Looking at a simple solution in line with what @rhattersley describes

[pelson]

Ok. So the underlying problem is that it is actually not possible to identify what is "in" and what is "out" with a spherical point-wise interpolation. If some assumptions are made about the nature of the input x and y coordinates (namely they are recti-linear) then it is possible, given some more information about whether the data "wraps" the globe or not, to mask out "target" pixels which were not within the range of the source coordinates.

In short, the stripy behaviour seen in http://stackoverflow.com/questions/19127697/difficulty-reprojecting-an-image-with-cartopy isn't incorrect on the sphere, though I do accept that it also isn't the desired behaviour. In order to get around this problem in the short term, I think we could add an experimental "recilinear_extrapolation" keyword which toggles whether the new masking takes place. Given the biggest usecase is for reprojecting images, which are rectilinear in their source coordinate system, I'd also suggest that the keyword is False by default but that GeoAxes.imshow always sets the value to True.

The only remaining complication, then, is that passing images to GeoAxes.imshow which are in the range 0-360 will result in an image which is only half of the longitudinal extent. The correct answer to this IMHO is that it is not reasonable to expect this functionality, and instead a PlateCarree projection with an appropriate central longitude should be sought.

I think this leaves us with the best of both world at the cost of a boolean keyword argument (sorry @rhattersley). Anyone care to comment, or even submit the PR? Otherwise I'll take a look at this early this afternoon.

[rhattersley]

BUG: Fix masks for specific cases.

This is very specific in terms of what it will fix, (and I'm not madly fond of all the somewhat opaque indexing going on, e.g. [0][0]), but it's a step in the right direction. I'm tempted to take it as it is.

[pelson]

I'm tempted to take it as it is.

Please see my (fairly lengthy) comment. We currenty have no drive to "fix" the masking as far as I can tell - it either needs to be on for image resampling (via imshow) or off for those calling regrid directly (via iris most likely). Anything else is probably speculative at best - I think we can do the simple thing here and get just as much benefit.

[cpelley]

getattr(source_cs, 'is_geodetic', False) would be equivalent

But is_geodetic is a property, sorry if I'm missing something obvious.

[cpelley]

Have pushed up some changes, thanks for looking at this @rhattersley and @pelson

[pelson]

I'm not sure I've got the name of this keyword correct. Care to make another suggestion?

[cpelley]

rectilinear on its own?

[pelson]

What part of the behaviour is being toggled?

[cpelley]

The masking of points falling outside the source domain, but I'm not sure a single word would do it. The keyword has to have rectilinear in it. How about extrapolate_rectilinear_limits?

[pelson]

Perhaps mask_extrapolated and add that this only makes sense for rectilinear source data in the docstring?

[cpelley]

OK

[cpelley]

Updated and squashed

[pelson]

The source coordinate is rectilinear in this function.

[pelson]

Ok. This is looking mergeable now. Anyone else care to comment? Otherwise I'll get on and merge this tomorrow.

Cheers,

[ajdawson]

Only comment is that the copyright dates probably should be updated before a release (I forgot to do that previously too).

[cpelley]

Ah yes, thanks @ajdawson

[cpelley]

All done + squashed