用 cartopy 和 matplotlib.tri 绘制三角网格

Question

我正在尝试使用 cartopy 和 matplotlib.tri 绘制一个三角形网格。我使用matplotlib.tri.Triangulation 对象并想用matplotlib.pyplot.triplot 绘制它。

当我使用transform=projection 将cartopy 投影作为转换传递给triplot 时，就像我绘制一条线时所做的那样，并非所有三角形都被绘制出来，并且在尝试保存图形时我得到一个IndexError。

另一方面，当我手动转换三角剖分中的所有点并调用 triplot 时，它可以工作。

在下面的第一个示例中，绘制了所有三角形，Line2D 对象没有区别。但是试图保存这个数字会引发IndexError。

在第二个示例中，来自带有transform=projection 的triplot 和来自triplot 的手动转换后生成的Line2D 对象在应用于其数据的转换结果方面有所不同。所以可能是组合变换的顺序问题。

注意：我选择 PlateCarree 是为了让它尽可能简单。我对轴和数据投影的其他组合也有同样的问题。

import matplotlib.pyplot as plt
import numpy as np
from cartopy.crs import PlateCarree
from matplotlib.tri import Triangulation

plt.interactive(False)

# NB. plt.triplot returns a list of two Line2D objects in both cases
# the second of which is empty, therefore only the first is returned

def triplot_with_transform(triangulation):
    """
    triangulation: matplotlib.tri.Triangulation

    """
    plt.figure()
    plt.subplot(projection=PlateCarree())
    lines = plt.triplot(triangulation, transform=PlateCarree())
    return plt.gcf(), lines[0]


def transform_before_triplot(triangulation):
    """
    triangulation: matplotlib.tri.Triangulation

    """
    plt.figure()
    plt.subplot(projection=PlateCarree())
    [x_tr, y_tr, _] = PlateCarree().transform_points(PlateCarree(),
                                                     triangulation.x,
                                                     triangulation.y).T
    triangulation_tr = Triangulation(x_tr, y_tr, triangulation.triangles)
    lines = plt.triplot(triangulation_tr)
    return plt.gcf(), lines[0]


def compare_line2d(line1, line2):
    """
    line1, line2 : matplotlib.lines.Line2D

    """
    data1 = line1.get_xydata()
    transform1 = line1.get_transform()
    data2 = line2.get_xydata()
    transform2 = line2.get_transform()
    data1_ma = np.ma.masked_invalid(data1)
    data2_ma = np.ma.masked_invalid(data2)
    data1_tr = transform1.transform(data1)
    data2_tr = transform2.transform(data2)
    data1_tr_ma = np.ma.masked_invalid(data1_tr)
    data2_tr_ma = np.ma.masked_invalid(data2_tr)

    print 'NaNs in line data match ', np.all(data1_ma.mask == data2_ma.mask)
    print 'Line data mismatch', abs(data1_ma-data2_ma).max()

    print 'NaNs in transformed line data match ', np.all(data1_tr_ma.mask ==
                                                         data2_tr_ma.mask)
    print 'Transformed line data mismatch', abs(data1_tr_ma-data2_tr_ma).max()


def try_so_save(fig):
    try:
        fig.savefig('triangulation.pdf')
    except Exception as e:
        print "Exception while saving figure\n", e.__repr__()


# First example, 100 vertices.
x, y = np.meshgrid(np.arange(0, 50, 5), np.arange(0, 50, 5))
tri = Triangulation(x.ravel(), y.ravel())
fig1, line1 = triplot_with_transform(tri)
fig2, line2 = transform_before_triplot(tri)

# Second example, 2500 vertices.
x, y = np.meshgrid(np.arange(0, 50), np.arange(0, 50))
tri = Triangulation(x.ravel(), y.ravel())
fig3, line3 = triplot_with_transform(tri)
fig4, line4 = transform_before_triplot(tri)

# Compare Line2D objects; since we transform from PlateCarree to
# PlateCarree, the lines should be (almost) equal.
print "Line2D objects from first example"
compare_line2d(line1, line2)
print "Line2D objects from second example"
compare_line2d(line3, line4)

# Try to save the figures.
print "Save fig1 from first example."
try_so_save(fig1)
print "Save fig3 from second example."
try_so_save(fig3) 

# failing plt.savefig command to produce full traceback
fig1.savefig('triangulation.pdf')

输出是

Line2D objects from first example
NaNs in line data match  True
Line data mismatch 3.5527136788e-15
NaNs in transformed line data match  True
Transformed line data mismatch 0.0
Line2D objects from second example
NaNs in line data match  True
Line data mismatch 7.1054273576e-15
NaNs in transformed line data match  True
Transformed line data mismatch 3545.955
Save fig1 from first example.
Exception while saving figure
IndexError('Out of bounds on buffer access (axis 0)',)
Save fig3 from second example.
Exception while saving figure
IndexError('Out of bounds on buffer access (axis 0)',)

当尝试保存使用transform=projection 创建的任何图形时，cartopy._crs.CRS.transform_points 会引发IndexError。

这是 cartopy 还是 matplotlib.tri 的问题？有没有办法避免手动转换并将transform 参数传递给triplot？

编辑：

在pp-mo的回答之后，我发现在triplot命令中添加一个带有标记的格式字符串，比如'o-'，可以让它工作。不过，仅仅改变颜色 ('g-) 并没有帮助。

编辑 2：

在示例代码中添加了fig1.savefig 以产生完整的回溯。

Answer 1

我试过了，机制似乎没有根本问题。 这是一个简单的工作示例...

"""
Example shamelessly poached from http://matplotlib.org/mpl_examples/pylab_examples/triplot_demo.py

We just add a mapping transform to the second figure there...

"""
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import math

import cartopy.crs as ccrs

# Some points defining a triangulation over (roughly) Britain.
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]])
# Make lats + lons
x = xy[:, 0]*180/3.14159
y = xy[:, 1]*180/3.14159

# A selected triangulation of the points.
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]])

# Make a triangulation object
my_tris = tri.Triangulation(x, y, triangles)

# Plot over an OSGB map with a coastline.
crs_pc = ccrs.PlateCarree()
crs_uk = ccrs.OSGB()
plt.figure()
ax = plt.axes(projection=crs_uk)
ax.coastlines(color='red', linewidth=1.5)
plt.triplot(my_tris,
            'go-',
            transform=ccrs.PlateCarree())
plt.show()

输出：

所以...
我只是猜测您的问题可能是当点通过变换时，您的实际三角剖分（哪些点以三角形连接）不再合适。
最明显的原因是三角形穿过投影中的“接缝”。例如，您可能有从 -30 到 +30 的经度，但想要在 0 到 360 度的地图上绘图。

这是可能的原因吗？