rendering issue with geodetic projection

[rsignell-usgs]

here's an interesting one. When we render this triangular mesh using PlateCarree it works fine, but using the WGS84 globe geodetic projection, some of the triangles end up as points instead:

Here's the notebook:
http://nbviewer.ipython.org/gist/rsignell-usgs/d101fec5ba8f6d177333

[ajdawson]

I'm not too surprised about this since cartopy does not implement its own triplot method on the GeoAxes class, so when you do ax.triplot it will use the matplotlib method without any cartopy-specific knowledge. I guess this means the transforms aren't applied properly...

I don't know if this is easy to work around or not, what do you reckon @pelson?

[pelson]

I don't know if this is easy to work around or not, what do you reckon @pelson?

Hmmm, interesting. It is curious that we get triangles for some situations though. At the end of the day, the lines being drawn should be trivial to get working given the machinery in place with cartopy.

I've had a go at reproducing with some canned data with no luck:

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


geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))

xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.045, 0.897], [-0.057, 0.895], [-0.073, 0.900], [-0.087, 0.898],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.073, 0.928], [-0.052, 0.930], [-0.048, 0.942], [-0.062, 0.949],
    [-0.054, 0.958], [-0.069, 0.954], [-0.087, 0.952], [-0.087, 0.959],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.097, 0.975], [-0.092, 0.984], [-0.101, 0.980], [-0.108, 0.980],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.101, 1.007], [-0.090, 1.010], [-0.087, 1.021], [-0.069, 1.021],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.048, 1.005], [-0.031, 1.005], [-0.031, 0.996], [-0.040, 0.987],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [-0.020, 0.954], [-0.006, 0.947], [ 0.003, 0.935], [ 0.006, 0.926],
    [ 0.005, 0.921], [ 0.022, 0.923], [ 0.033, 0.912], [ 0.029, 0.905],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.057, 0.881], [-0.062, 0.876], [-0.078, 0.876], [-0.087, 0.872],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

triangles = np.asarray([
    [67, 66,  1], [65,  2, 66], [ 1, 66,  2], [64,  2, 65], [63,  3, 64],
    [60, 59, 57], [ 2, 64,  3], [ 3, 63,  4], [ 0, 67,  1], [62,  4, 63],
    [57, 59, 56], [59, 58, 56], [61, 60, 69], [57, 69, 60], [ 4, 62, 68],
    [ 6,  5,  9], [61, 68, 62], [69, 68, 61], [ 9,  5, 70], [ 6,  8,  7],
    [ 4, 70,  5], [ 8,  6,  9], [56, 69, 57], [69, 56, 52], [70, 10,  9],
    [54, 53, 55], [56, 55, 53], [68, 70,  4], [52, 56, 53], [11, 10, 12],
    [69, 71, 68], [68, 13, 70], [10, 70, 13], [51, 50, 52], [13, 68, 71],
    [52, 71, 69], [12, 10, 13], [71, 52, 50], [71, 14, 13], [50, 49, 71],
    [49, 48, 71], [14, 16, 15], [14, 71, 48], [17, 19, 18], [17, 20, 19],
    [48, 16, 14], [48, 47, 16], [47, 46, 16], [16, 46, 45], [23, 22, 24],
    [21, 24, 22], [17, 16, 45], [20, 17, 45], [21, 25, 24], [27, 26, 28],
    [20, 72, 21], [25, 21, 72], [45, 72, 20], [25, 28, 26], [44, 73, 45],
    [72, 45, 73], [28, 25, 29], [29, 25, 31], [43, 73, 44], [73, 43, 40],
    [72, 73, 39], [72, 31, 25], [42, 40, 43], [31, 30, 29], [39, 73, 40],
    [42, 41, 40], [72, 33, 31], [32, 31, 33], [39, 38, 72], [33, 72, 38],
    [33, 38, 34], [37, 35, 38], [34, 38, 35], [35, 37, 36]])

ax = plt.axes(projection=ccrs.PlateCarree(0))
ax.coastlines()
plt.triplot(x, y, triangles, transform=geodetic,
            marker='.', linestyle='-', alpha=0.85, color='darkgray', )

plt.show()

I was also changing the projection (central longitude of 180, using Robinson etc.) and had no problem, and have also passed a Triangulation in rather than passing x, y, triangles.

Essentially there is a bug here, we just need to get down to a simple reproducible case. Any chance you'd be willing to dig a little deeper @rsignell-usgs so that we are data and external dependency free?

[rsignell-usgs]

I removed the dependency on ugrid, so it's just cartopy, matplotlib and netCDF4 now. That's okay, right?
http://nbviewer.ipython.org/gist/rsignell-usgs/c410055afd489e139128

[pelson]

I removed the dependency on ugrid, so it's just cartopy, matplotlib and netCDF4 now. That's okay, right?

Well, I was trying to get down to a case which didn't need any data so that we could use it directly as a unit test when we fix it.

[rsignell-usgs]

Okay, this one just generates some random triangles
http://nbviewer.ipython.org/gist/rsignell-usgs/685be5eff19ce15e7340

[rsignell-usgs]

The images produced by this bug are rather pretty, though.

Reminds me of this book I had as a kid.

[pelson]

Cool. Thanks @rsignell-usgs, and thanks for the reference - I agree, there is something beautiful in the result 😄

[pelson]

Looks like this has been introduced between v1.3 and 1.4 of matplotlib. We don't test against v1.4/1.5 at this point, but there is some work in progress to do so.

[QuLogic]

The simplified example appears to be fine now, at least on master. Are you still having trouble with this?

[QuLogic]

This looks like it was fixed by #1146.