SciTools / cartopy

Cartopy - a cartographic python library with matplotlib support
https://scitools.org.uk/cartopy/docs/latest
BSD 3-Clause "New" or "Revised" License
1.41k stars 359 forks source link

Discontinuity in Interupted Goode Homolosine #319

Closed mikewalker92 closed 10 years ago

mikewalker92 commented 11 years ago

The following code produces a plot containing a vertical discontinuity which appears to run through 0 longitude.

import iris import cartopy import matplotlib.pyplot as plt import iris.quickplot as qplt import cartopy.crs as ccrs

cube = iris.load_cube(iris.sample_data_path('A1B.2098.pp'))

plt.axes(projection=ccrs.InterruptedGoodeHomolosine())

qplt.pcolormesh(cube)

plt.gca().coastlines()

plt.show()

discontinuity

The plot also appears not to be plotting the correct data to the right of the discontinuity, as the uniform horizontal stripes do not look correct.

ajdawson commented 10 years ago

A bit of preliminary diagnosis, expert assistance is required! Using a dummy data set to reproduce without iris and adding a matplotlib plot:

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

x = np.arange(0, 360, 2.5)
y = np.arange(-89, 90, 2)
X, Y = np.meshgrid(*[np.deg2rad(c) for c in (x, y)])
Z = np.cos(Y) + 0.375 * np.sin(2. * X)
Z[:, 1:3] = 2.
Z[:, -4:-2] = -2.

plt.pcolormesh(x, y, Z)

plt.figure()
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())
ax.coastlines()
ax.pcolormesh(x, y, Z, transform=ccrs.PlateCarree())

data_0-360 goode_0-360

So as @mikewalker92 said there is a discontinuity, but if I had my data defined on the longitude interval -180 to 180:

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

x = np.arange(-180, 180, 2.5)
y = np.arange(-89, 90, 2)
X, Y = np.meshgrid(*[np.deg2rad(c) for c in (x, y)])
Z = np.cos(Y) + 0.375 * np.sin(2. * X)
Z[:, 1:3] = 2.
Z[:, -4:-2] = -2.

plt.pcolormesh(x, y, Z)

plt.figure()
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())
ax.coastlines()
ax.pcolormesh(x, y, Z, transform=ccrs.PlateCarree())

data_-180-180 goode_-180-180

Things look better (aside from data going outside the map boundary, which is likely a separate issue). Problem is I don't know how to fix it! But hopefully this might be helpful at least.

pelson commented 10 years ago

Very helpful. There is code in pcolormesh which looks for wrapping polygons (https://github.com/SciTools/cartopy/blob/master/lib/cartopy/mpl/geoaxes.py#L1109)[here]. Looks like it either isn't including Interrupted, or there is a case that it isn't handling properly.

ajdawson commented 10 years ago

Interestingly pcolor seems to work as expected. I guess they are rather different implementations under the hood though.

ajdawson commented 10 years ago

OK so I added ccrs.InterruptedGoodeHomolosine to warp_proj_types and the discontinuity is gone, which is good. However I can't figure out how to get it to stop drawing outside the map boundary. Any ideas on that one @pelson?

ajdawson commented 10 years ago

I managed to get a correct looking plot but I fear it is just a hack rather than a fix! As I said before I added ccrs.InterruptedGoodeHomolosine to warp_proj_types in the _pcolormesh_patched() method, and then just set the clip path of the resulting QuadMesh to the projection boundary with collection.set_clip_path(self.outline_patch) at the end of _pcolormesh_patched().

Is that sensible? It feels like a hack but I couldn't figure out what was really wrong. If it sounds reasonable I'll submit a PR.

goode_good

(blue and red bands are there on purpose to make sure the data values are in the correct location, correspond to the first plot in my first comment)

pelson commented 10 years ago

Interestingly pcolor seems to work as expected.

Yep. Each quadrilateral is transformed using cartopy's intelligent polygon projection algorithm with pcolor - that can make it rather slow though!

I have no issue with your proposed solution. The clip path is reasonable, though I would have to check it behaves sensibly when you drag the Axes around the window.

Good job @ajdawson.

ajdawson commented 10 years ago

though I would have to check it behaves sensibly when you drag the Axes around the window

What exactly do you mean by drag the Axes around the window?

pelson commented 10 years ago

What exactly do you mean by drag the Axes around the window?

In a TkAgg GUI (or another backend would be fine), just drag the resulting Axes so that it gets clipped by the Axes bounding box in the same way that coastlines would. If you put up a PR, I would be happy to check it out if that helps.

Cheers!

ajdawson commented 10 years ago

Ah OK you mean click the move button and then drag them round, I forgot you could do that! It behaves identically to the same plot drawn with contourf so I don't think it is a problem if that is what you meant. I'll submit a PR.

ajdawson commented 10 years ago

You may still experience problems with your Iris use case @mikewalker92, but I think that is covered by #325.