Cartopy 系列教程：3D 地图绘制

bugsuse

发布于 2020-04-21 17:47:38

2.4K0

发布于 2020-04-21 17:47:38

文章被收录于专栏：气象杂货铺气象杂货铺

Basemap 的 mpl3d 绘制3D地图时非常强大，但目前仍然存在一些小问题，比如在填充陆地时有时会出现问题。

虽然 Cartopy 中没有提供直接绘制 3D 地图的方法，但是使用 Cartopy 同样可以绘制 3D 地图。下面就逐步看一下如何绘制：

首先，获取 shp 文件及相应的几何图形（geometries）

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]

因为这些都是均匀的几何图形，我们需要将其转换为 matplotlib path

from cartopy.mpl.patch import geos_to_path
import itertools

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

完成这些后，只是创建了 PathCollection，然后要添加这些到 axes，但是 Axes3D 似乎并不支持 PathCollection 实例。因此，需要构建 LineCollection 实例，这也是 basemap 所使用的方法。因为 LineCollection 不支持 paths，所以需要使用 segments：

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 似乎能很好的解决 PathCollection 问题，后面我们可以试一下填充 polygon:

最重要的一点是转换 paths 为 polygon，然后传递给 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()

结果如下图：