类似于Matplotlib-Basemap的3D CartoPy

Question

我刚开始学习Python，想问一下关于Cartopy能否用于3D绘图的问题。下面是一个使用matplotlibBasemap的例子。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.basemap import Basemap

m = Basemap(projection='merc',
            llcrnrlat=52.0,urcrnrlat=58.0,
            llcrnrlon=19.0,urcrnrlon=40.0,
            rsphere=6371200.,resolution='h',area_thresh=10)

fig = plt.figure()
ax = Axes3D(fig)
ax.add_collection3d(m.drawcoastlines(linewidth=0.25))
ax.add_collection3d(m.drawcountries(linewidth=0.35))
ax.add_collection3d(m.drawrivers(color='blue'))

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Height')

fig.show()

这个例子创建了一个3D坐标系中的地图，这样你就可以在地图表面上绘制物体。但是使用Cartopy时，它返回的是matplotlib.axes.GeoAxesSubplot。我不太清楚如何将这个结果添加到像上面那样的3D图形/坐标系中，使用matplotlib-basemap。

所以，有人能给我一些关于如何用Cartopy做类似3D绘图的建议吗？

Answer 1

basemap mpl3d 是一个很有趣的技巧，但它并不是为了这样使用而设计的。因此，目前你只能用它来处理简单的海岸线。比如说，填充的大陆就不太行，至少我觉得是这样。

不过，使用 cartopy 时有一个类似的技巧可以用。因为我们可以通用地访问 shapefile 信息，所以这个方法应该适用于任何多边形线的 shapefile，比如海岸线。

第一步是获取 shapefile 以及相关的几何信息：

feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m')
geoms = feature.geometries()

接下来，我们可以把这些信息转换成想要的投影：

target_projection = ccrs.PlateCarree()
geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]

因为这些是 shapely 几何体，所以我们接下来想把它们转换成 matplotlib 的路径，使用：

from cartopy.mpl.patch import geos_to_path
import itertools

paths = list(itertools.chain.from_iterable(geos_to_path(geom)
                                             for geom in geoms))

有了路径后，我们应该能在 matplotlib 中创建一个 PathCollection，并把它添加到坐标轴上，但可惜的是，Axes3D 似乎不支持 PathCollection 实例，所以我们需要通过构建一个 LineCollection 来解决这个问题（就像 basemap 那样）。可惜的是，LineCollections 不接受路径，而是接受线段，我们可以通过以下方式计算这些线段：

segments = []
for path in paths:
    vertices = [vertex for vertex, _ in path.iter_segments()]
    vertices = np.asarray(vertices)
    segments.append(vertices)

把这些都整合在一起，我们就能得到一个和你代码生成的 basemap 图类似的结果：

import itertools

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import numpy as np

import cartopy.feature
from cartopy.mpl.patch import geos_to_path
import cartopy.crs as ccrs


fig = plt.figure()
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90])
ax.set_zlim(bottom=0)


target_projection = ccrs.PlateCarree()

feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m')
geoms = feature.geometries()

geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]

paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms))

# At this point, we start working around mpl3d's slightly broken interfaces.
# So we produce a LineCollection rather than a PathCollection.
segments = []
for path in paths:
    vertices = [vertex for vertex, _ in path.iter_segments()]
    vertices = np.asarray(vertices)
    segments.append(vertices)

lc = LineCollection(segments, color='black')

ax.add_collection3d(lc)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Height')

plt.show()

另外，mpl3d 似乎对 PolyCollection 的处理也不错，这可能是我会考虑用来处理填充几何体的方式，比如土地轮廓（而不是海岸线，海岸线只是一个轮廓）。

关键的一步是把路径转换成多边形，并在 PolyCollection 对象中使用这些多边形：

concat = lambda iterable: list(itertools.chain.from_iterable(iterable))

polys = concat(path.to_polygons() for path in paths)
lc = PolyCollection(polys, edgecolor='black',
                    facecolor='green', closed=False)

这个案例的完整代码大概是这样的：

import itertools

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection, PolyCollection
import numpy as np

import cartopy.feature
from cartopy.mpl.patch import geos_to_path
import cartopy.crs as ccrs


fig = plt.figure()
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90])
ax.set_zlim(bottom=0)


concat = lambda iterable: list(itertools.chain.from_iterable(iterable))

target_projection = ccrs.PlateCarree()

feature = cartopy.feature.NaturalEarthFeature('physical', 'land', '110m')
geoms = feature.geometries()

geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]

paths = concat(geos_to_path(geom) for geom in geoms)

polys = concat(path.to_polygons() for path in paths)

lc = PolyCollection(polys, edgecolor='black',
                    facecolor='green', closed=False)

ax.add_collection3d(lc)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Height')

plt.show()

最终效果是：