2

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

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

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

在下面的第一个示例中,所有三角形都被绘制并且Line2D对象没有区别。但试图挽救这个数字会引发IndexError.

在第二个示例中,使用trilot 和手动转换后Line2D来自 triplot 的结果对象transform=projection在应用于其数据的转换结果方面有所不同。所以也许这是组合转换顺序的问题。

注意:我选择 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=projectioncartopy._crs.CRS.transform_points会引发IndexError

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

编辑:

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

编辑2:

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

4

1 回答 1

0

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

"""
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 度的地图上绘图。

这是一个可能的原因吗?

于 2015-09-02T17:50:58.330 回答