使用matplotlib更高效地绘制多边形

Question

我有一个数据集，大约包含60000个形状（每个形状都有各自的经度/纬度坐标），我想使用matplotlib和basemap在地图上绘制这些形状。

目前我是这样做的：

for ii in range(len(data)):
    lons = np.array([data['lon1'][ii],data['lon3'][ii],data['lon4'][ii],data['lon2'][ii]],'f2')
    lats = np.array([data['lat1'][ii],data['lat3'][ii],data['lat4'][ii],data['lat2'][ii]],'f2')
    x,y = m(lons,lats)
    poly = Polygon(zip(x,y),facecolor=colorval[ii],edgecolor='none')
    plt.gca().add_patch(poly)

然而，在我的计算机上这需要大约1.5分钟，我在思考是否有可能稍微加快一些速度。有没有更有效的方法来绘制多边形并将它们添加到地图上？

- HyperCube

Answer 1

你可以考虑创建多边形的集合，而不是单个的多边形。

相关文档可在这里找到：http://matplotlib.org/api/collections_api.html ，其中有一个值得分析的例子：http://matplotlib.org/examples/api/collections_demo.html 例如：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
import matplotlib as mpl

# Generate data. In this case, we'll make a bunch of center-points and generate
# verticies by subtracting random offsets from those center-points
numpoly, numverts = 100, 4
centers = 100 * (np.random.random((numpoly,2)) - 0.5)
offsets = 10 * (np.random.random((numverts,numpoly,2)) - 0.5)
verts = centers + offsets
verts = np.swapaxes(verts, 0, 1)

# In your case, "verts" might be something like:
# verts = zip(zip(lon1, lat1), zip(lon2, lat2), ...)
# If "data" in your case is a numpy array, there are cleaner ways to reorder
# things to suit.

# Color scalar...
# If you have rgb values in your "colorval" array, you could just pass them
# in as "facecolors=colorval" when you create the PolyCollection
z = np.random.random(numpoly) * 500

fig, ax = plt.subplots()

# Make the collection and add it to the plot.
coll = PolyCollection(verts, array=z, cmap=mpl.cm.jet, edgecolors='none')
ax.add_collection(coll)
ax.autoscale_view()

# Add a colorbar for the PolyCollection
fig.colorbar(coll, ax=ax)
plt.show()

希望有所帮助，

- pelson

---

谢谢，很好的例子！我认为PolyCollection是关键。然而，我不知道如何将我的经纬度转换成多边形。在你的例子中是“verts”。 - HyperCube

@JoeKington：非常好的补充。不幸的是，我将得到你所有辛勤工作的功劳... - pelson

同样的限制...绘制27000个多边形需要15秒。不算太糟糕，但有没有什么方法可以改进呢？ 另外，一旦多边形被绘制出来，是否有一种方法只更新值而不是重新绘制？ - PBrockmann

Answer 2

我调整了我的代码，现在它完美无缺地工作着：）

下面是可行的示例：

lons = np.array([data['lon1'],data['lon3'],data['lon4'],data['lon2']])
lats = np.array([data['lat1'],data['lat3'],data['lat4'],data['lat2']])
x,y = m(lons,lats)
pols = zip(x,y)
pols = np.swapaxes(pols,0,2)
pols = np.swapaxes(pols,1,2)
coll = PolyCollection(pols,facecolor=colorval,cmap=jet,edgecolor='none',zorder=2)
plt.gca().add_collection(coll)

- HyperCube

---

完美无缺并且快速吗？你比原来的1.5分钟节省了多少时间？ - pelson

现在只需要32秒，所以它真的加快了事情的进展！ - HyperCube

m是什么？添加导入/定义会更好。 - Skylar Saveland

@SkylarSaveland m是mpl_toolkits.basemap模块中的Basemap对象。它可以将经纬度转换为绘制地图的轴单位。 - Syrtis Major