本节提要：关于如何利用matplotlib cartopy绘制酷炫的三维地图。
这是我在比较久远之前看到的问题。首先必须明确一点，matplotlib的axes3D这个投影中 ，是不能用add_geometry这个功能来直接将读取到的shp文件添加上去的。add_geometry这个功能是cartopy下的geoaxes才能使用，同理add_feature也不能再3d图中使用。
但是这个功能确实又是比较常用，而且酷炫的。所以我在s站上查到了cartopy库包的开发人员直接给出的回答。我们不直接开讲怎么绘制，而先回忆在不久之前的推文中，我们使用过的一个功能。
在气象绘图加强版（八）——框线的添加这一节中，我们说了有几种添加框线的方式：polygon、plot、add_geometry。我们添加框线时的add_geometry方法是不是读取的shp文件呢，其实不是。而是使用了shapely的一个功能，构造了一个geometry
代码语言：javascript复制import shapely.geometry as sgeom
x2=[105,115,115,105,105]
y2=[25,25,35,35,25]
plot_geometries= sgeom.LineString(zip(x2,y2))
ax.add_geometries([plot_geometries],proj,facecolor='none',
edgecolor='k',lw=0.3,ls=":")通过这么简单的一个问题，引出了这么一个变化过程，即matplotlib中的plot、matplotlib中的polygon、地图geometry（几何图形）是可以相互转化的，他们本质上是横纵坐标下的点线面。
恰巧，matplotlib的axes3D投影中，允许我们使用polygon功能。那么我们是不是可以将shp文件中的geometry读取出来，转变成polygon，然后添加到三维图中呢。
这就是我理解的开发者的想法。下面就是如何将shp文件中的geometry转化成polygon的问题了。cartopy有没有这个功能呢，打开官网文档，可以查到这么一条：
在这一行下，有一个geos_to_path，如果我的理解正确，这个功能全称应该是geometries to path。而在matplotlib中，path是可以转化成patch的，而polygon就是patch的一种。但是显然，上述只是一个简单的解释为什么这几样之间可以相互转化。由于shp文件中的geometry信息众多，所以使用的是polycollection（多边形集合）功能，这个功能能够存储全部polygon，并一次性绘制。
看完以上内容，我再来贴出开发者的原始程序就好理解多了：
代码语言：javascript复制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
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False#负号
fig = plt.figure(figsize=(10,8),dpi=200)
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90])
ax.set_zlim(0,0.5)
##############################################################
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) #geom转path
polys = concat(path.to_polygons() for path in paths) #path转poly
lc = PolyCollection(polys, edgecolor='black',
facecolor='yellow', closed=False)
ax.add_collection3d(lc)
ax.set_xlabel('经度')
ax.set_ylabel('纬度')
ax.set_zlabel('高度')
plt.show()
使用时，只需要修改shp文件的路径，读取本地文件。限制绘图范围即可。这个办法读取的地图已经是polygon格式而不是geomtery格式的了。通言之，在平面图中也是可以读取并放置的：
代码语言：javascript复制ax.add_collection(lc)
但是，这个开发者给出的原始版本答案存在一个问题，即你在收缩3d投影的长宽时polygon不会改变导致地图突出的问题：
代码语言：javascript复制ax = Axes3D(fig, xlim=[-70, 130], ylim=[0, 70])
由于这个polygon的范围是按照依据的shp生成的，并不会按照3d投影中的方式来收缩。于是开发者提出了一个更加完备的封装函数完成收缩功能。
代码语言：javascript复制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 as cf
from cartopy.mpl.patch import geos_to_path
import cartopy.crs as ccrs
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False#负号
fig = plt.figure(figsize=(10,8),dpi=200)
ax = Axes3D(fig, xlim=[70, 130], ylim=[-0, 70])
ax.set_zlim(0,0.5)
####################预先设置地图的参数######################################
proj_ax=plt.figure().add_subplot(111,projection=ccrs.PlateCarree())
proj_ax.set_xlim(ax.get_xlim())#使地图投影获得当前3d投影一样的绘图范围
proj_ax.set_ylim(ax.get_ylim())
concat = lambda iterable: list(itertools.chain.from_iterable(iterable))
target_projection=proj_ax.projection
feature=cf.NaturalEarthFeature('physical', 'land', '50m')
geoms=feature.geometries()
boundary=proj_ax._get_extent_geom()
geoms = [target_projection.project_geometry(geom, feature.crs)
for geom in geoms]
geoms2=[]
for i in range(len(geoms)):
if geoms[i].is_valid:
geoms2.append(geoms[i])
geoms=geoms2
geoms=[boundary.intersection(geom)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='yellow', closed=False)
ax.add_collection3d(lc)
proj_ax.spines['geo'].set_visible(False)#解除掉用于确定地图的子图
plt.show()
这些步骤里，只需替换shp文件即可复制粘贴使用。我用了手头的五六份文件，都没有报错过。
只要前贤造了能用的工具，我这个懒鬼是肯定不会想新方法的。
ax.add_collection3d(lc)这个函数里有个参数zs参数控制添加时的位置，我们修改这个参数可使其放置在不同的高度：
代码语言：javascript复制lc = PolyCollection(polys, edgecolor='black',
facecolor='yellow', closed=False)
lc2 = PolyCollection(polys, edgecolor='black',
facecolor='green', closed=False)
lc3 = PolyCollection(polys, edgecolor='black',
facecolor='white', closed=False)
ax.add_collection3d(lc,zs=0)
ax.add_collection3d(lc2,zs=0.25)
ax.add_collection3d(lc3,zs=0.5)
另外，由于ax.add_collection3d还有一个参数zdir,他使我们可以更改更加酷炫的地图添加方式。例如zdir传入的一个y，表示当前添加的poly平面与y轴垂直，其他可自行摸索：
代码语言：javascript复制lc = PolyCollection(polys, edgecolor='black',
facecolor='yellow', closed=False)
lc3 = PolyCollection(polys, edgecolor='black',
facecolor='white', closed=False)
lc4 = PolyCollection(polys, edgecolor='black',
facecolor='orange', closed=False)
ax.add_collection3d(lc3,zs=-90)
ax.add_collection3d(lc4,zdir='y',zs=180)
ax.add_collection3d(lc,zdir='x',zs=-180)
而且，我们还可以引入cartopy中的经纬度格式化器来格式化3d地图中的经纬度。
代码语言：javascript复制from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xlabel('经度')
ax.set_ylabel('纬度')
ax.set_zlabel('高度')
接下来，我们简单介绍一下如何在3d图中使用contourf函数绘制平面图。
由于3d图的投影结构完全与当前我们的世界相符合，都是三维空间，所以3d图中的contourf与真实世界等值线相同都是立体的，这与我们平时见到的二维等值线图不一致。
代码语言：javascript复制def data():
x=np.arange(80,120)
y=np.arange(0,60)
X,Y=np.meshgrid(x,y)
Z=(X-100)**2/100 (Y-30)**2/10
return X,Y,Z
X,Y,Z=data()
ax.contourf(X,Y,Z)
能够看出来，这是现实中等值线填色图存在的意义，z值即代表其存在的高度。如果用面绘图命令能够更清楚的看出等值线代表的实际意义:
代码语言：javascript复制ax.plot_surface(X,Y,Z,cmap='viridis')
这大概就是天气学上高度场的槽的立体模型。
不过显然我们不需要这样立体的表示这个等值线场。利用在前面提到的zdir和offset命令，我们就可以修改三维图的等值线图变为二维等值线图，zdir通俗理解就是将传入的维度降维，offset表示降维后的等值线填色图放置的层次：
代码语言：javascript复制ax.contourf(X,Y,Z,zdir='z',offset=0)
ax.contourf(X,Y,Z,zdir='z',offset=0)
ax.contourf(X,Y,Z,zdir='z',offset=30)
ax.contourf(X,Y,Z,zdir='z',offset=80)
ax.contourf(X,Y,Z,zdir='z',offset=0)
ax.contourf(X,Y,Z,zdir='y',offset=70)
ax.contourf(X,Y,Z,zdir='x',offset=70)
上述三种降维方式，z降维就是我们平时看到的那种平面等值线填色图。其他两种与我们的习惯不符合，所以很少使用。
