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()
结果如下图:
以上内容得益于 stackoverflow。
注:http://stackoverflow.com/questions/23785408/3d-cartopy-similar-to-matplotlib-basemap/23914781#23914781