气象绘图——复杂的三维图
本节提要：关于在三维绘图环境下，进行复杂的图形描绘。本节文章很长。
修改普通的三维图固定设置
在普通的matplotlib的三维投影中，我们似乎并不能获得我们想要的结果，尤其是视觉上的，虽然倾斜了图形，但是文字等标注仍然是二维的，例如下面这张图片：
明显，无论是刻度标记，还是文字标签，都呈现出与三维图像格格不入的问题，反而像是二维数据窜入了三维世界，所以为了绘制更加漂亮的图片，我们需要将这些更改为具有合适视觉效果的图形。这里主要依靠axes3D自身的相关命令，与matplotlib的demo——Draw flat objects in 3D plot。
关闭默认标签与网格
由于三维图默认的设置不美观，我们必须将其全部删除，使用下面这些语句完成这个需求：
ax.grid(False)
ax.xaxis.pane.fill=False
ax.yaxis.pane.fill=False
ax.zaxis.pane.fill=False
ax.xaxis.pane.set_edgecolor('none')
ax.yaxis.pane.set_edgecolor('none')
ax.zaxis.pane.set_edgecolor('none')
ax.xaxis._axinfo['tick']['outward_factor']=0
ax.xaxis._axinfo['tick']['inward_factor']=0.25
ax.yaxis._axinfo['tick']['outward_factor']=0
ax.yaxis._axinfo['tick']['inward_factor']=0.25
ax.zaxis._axinfo['tick']['outward_factor']=0
ax.zaxis._axinfo['tick']['inward_factor']=0.25
从上往下依次是关闭grid网格、关闭三轴格式与颜色、重新设置刻度系统。这样，我们可以得到一个没有灰色背景与网格的纯色三维图，如果没有更高的要求，到这里其实已经很素净了，不过我们的要求不止这些。
随后，我们调用Draw flat objects in 3D plot里面的代码段，将我们需要的打印的文字与标签变为3D影像投影在平面上的样式：
def text3d(ax, xyz, s, zdir="z", size=None, angle=0, usetex=False, **kwargs):
x, y, z = xyz
if zdir == "y":
xy1, z1 = (x, z), y
elif zdir == "x":
xy1, z1 = (y, z), x
else:
xy1, z1 = (x, y), z
text_path = TextPath((0, 0), s, size=size, usetex=usetex)
trans = Affine2D().rotate(angle).translate(xy1[0], xy1[1])
p1 = PathPatch(trans.transform_path(text_path), **kwargs)
ax.add_patch(p1)
art3d.pathpatch_2d_to_3d(p1, z=z1, zdir=zdir)
text3d(ax, (105,2,1000),'Longitude', zdir="z",size=3,usetex=False,
ec="none", fc="k")
text3d(ax, (140,22,1000),'Latitude', zdir="z",size=3,angle=np.pi*0.5,usetex=False,
ec="none", fc="k")
for m in np.arange(90,140,10):
text3d(ax, (m-2,7,1000), "{}".format(m)+'°E', zdir="z",size=2,usetex=False,
ec="none", fc="k")
ax.plot([m,m],[9,10],[1000,1000],color='k',lw=1,zorder=5)
for n in np.arange(10,60,10):
text3d(ax, (133,n,1000), "{}".format(n)+'°N', zdir="z",size=2,usetex=False,
ec="none", fc="k")
ax.plot([130,131],[n,n],[1000,1000],color='k',lw=1,zorder=5)
ax.plot([90,130],[10,10],[1000,1000],color='k',lw=1,zorder=5)
ax.plot([130,130],[10,50],[1000,1000],color='k',lw=1,zorder=5)
ax.axis('off')
是不是具有三维视感而不是伪劣的窜入效果了呢。
在三维图中导入地图
三维图中插入地图 利用左边这个链接里的相关介绍，我们给三维图加上地图：
proj_ax=plt.figure().add_subplot(111,projection=ccrs.PlateCarree())
proj_ax.set_xlim(ax.get_xlim())
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='lightgrey',alpha=1,closed=False)
ax.add_collection3d(lc,1000)
proj_ax.spines['geo'].set_visible(False)
这样，具有真三维效果的图形就做好了。
在三维图中实现栅格可视化
在之前的推文三维图形迁移中，我们已经介绍了如何使用收集collection的办法，来实现贴瓷砖式的数据可视化，这里我们仍然使用这种办法。假定使用FNL的再分析资料，精度为1×1。取出相对湿度的值进行剖面与平面图的绘制，并裁剪数据的轮廓。
RH=xr.open_dataset(file,engine='cfgrib',backend_kwargs={'filter_by_keys': {'typeOfLevel': 'isobaricInhPa','cfVarName': 'r'}})
lon=RH['longitude'].loc[90:130]
lat=RH['latitude'].loc[50:10]
rh=RH['r'].loc[1000:300,50:10,90:130]
lev=RH['isobaricInhPa'].loc[1000:300]
cmap=mpl.cm.plasma_r
colorlevel=np.arange(70,101,1)
norm=mcolors.BoundaryNorm(colorlevel,cmap.N)
ax2=fig.add_axes([0,0,0.5,0.5])
X,Y=np.meshgrid(lat,lev)
a1=ax2.pcolormesh(X,Y,rh.loc[:,:,90],shading='auto',cmap=cmap,norm=norm)
rh_flat=np.array(rh.loc[:,:,90]).flatten()
paths1=a1.get_paths()
polys1=concat(path.to_polygons() for path in paths1)
lc01=PolyCollection(polys1,cmap=cmap,edgecolor='none',closed=False,array=rh_flat)
ax.add_collection3d(lc01,zdir='x',zs=89)
##################################
X,Y=np.meshgrid(lon,lev)
a2=ax2.pcolormesh(X,Y,rh.loc[:,50,:],shading='auto',cmap=cmap,norm=norm)
rh_flat=np.array(rh.loc[:,50,:]).flatten()
paths2=a2.get_paths()
polys2=concat(path.to_polygons() for path in paths2)
lc02=PolyCollection(polys2,cmap=cmap,edgecolor='none',closed=False,array=rh_flat)
ax.add_collection3d(lc02,zdir='y',zs=51)
ax2.set_visible(False)
上述代码表示将数据分为两轮添加到地图上，第一次表示提取出沿东经90进行剖面，相对湿度沿纬向的分布；第二次表示提取沿北纬50进行剖面，相对湿度在经向上的分布。
从上图不难看出，在东经90°剖面上，相对湿度在北纬25°以北突然降低到非常小的范围；在北纬50°剖面上，相对湿度在东经105°-115°有一个深厚的湿区。由于我们使用的是pcolormesh函数，所以所有的栅格类数据都可以这样进行剖面可视化，经过与平面出图对比，应该是没有多大问题的。
使用plot_surface命令栅格化
在当前的三维投影中，暂时没有axes3D.pcolormesh这个平面图中常用的栅格化绘图函数，但是，我们可以使用plot_surface命令替代这个效果。plot_surface是一个通过拼接polygon来实现立体可视化效果，具体可见李开元老师绘制的一张假相当位温图：
具体来说，与contourf函数类似，x,y负责经纬定位，z值不仅负责垂直定位，还负责给曲面上色。不过此处有一个冲突，即我们需要平面，不需要z变化起伏。查阅plot_surface函数可知，默认情况下，z确定曲面上色，但是，提供了facecolors参量来超越z值上色。这时不妨一个想法，如果我将z全部设为一个值，则曲面变为平面，不妨新建一个z数组：
LON,LAT=np.meshgrid(lon,lat)
cmap=plt.get_cmap('viridis_r')
rh_850=rh.loc[850].values
z=np.full((LON.shape[0],LAT.shape[1]),850)
ax.plot_surface(LON,LAT,z,cmap=cmap)
我们提取850百帕上的相对湿度，使z全部填充为850,然后绘图：
可以看出，plot_surface实现了pcolormesh效果，颜色完全一样，是因为全部z都为850，默认使用z为颜色映射依据，则处处相等，故有此现象。下面，我们来新建一个facecolor数据，来给每个小方框上色：
facecolor_array=np.empty((LON.shape[0],LAT.shape[1]),dtype=tuple)
for i in range(LON.shape[0]):
for j in range(LON.shape[1]):
facecolor_array[i,j]=cmap((rh_850[i,j]/100))
z=np.full((LON.shape[0],LAT.shape[1]),850)
ax.plot_surface(LON,LAT,z,facecolors=facecolor_array,cmap=cmap)
效果还是可以的。上述两种栅格化，具体有什么用，目前视觉效果最好的就是这一种：
在三维图中实现contourf可视化
我们之前曾经推送过如何进行contourf的三维可视化，但是有一定的问题，小值区的色块总是会被遮盖，幸运的是，matplotlib官网新上线了一个demo——3D box surface plot，可以解决我们很多问题。 首先导入我们需要的数据：
RH=xr.open_dataset(file,engine='cfgrib',backend_kwargs={'filter_by_keys': {'typeOfLevel': 'isobaricInhPa','cfVarName': 'r'}})
lon=RH['longitude'].loc[90:130]
lat=RH['latitude'].loc[50:10]
rh=RH['r'].loc[1000:300,50:10,90:130]
lev=RH['isobaricInhPa'].loc[1000:300]
随后，类似于二维平面的网格化，对三维坐标的lon,lat,level进行网格化：
Y,Z,X=np.meshgrid(lat,lev,lon)
因为这段程序使用的是我当时学习的原始程序，所以网格化顺序严格与demo相同，后期可以不使用这种网格化顺序。随后使用contourf的投影参数:
kw={'vmin':rh.min(),
'vmax':rh.max(),
'levels':np.arange(0,110,10)}
cross_lon_90=ax.contourf(
rh[:,:,0].values, Y[:,:,0], Z[:,:,0],cmap='Blues',
zdir='x', offset=89.9, **kw)
cross_lat_50=ax.contourf(
X[:,0,:], rh[:,0, :].values, Z[:,0, :],cmap='Blues',
zdir='y', offset=50.1, **kw)
在三维图中使用数据白化
接下来，我们不妨考究这一张图片，以这一张图片看看其主要元素怎么使用单纯使用matplotlib进行仿制。首先看底层，左侧为青藏高原地形下的数据，且仅含有青藏高原数据。右侧为依靠经线为分界线，使东经100°左侧无数据；中层为两个链接柱；上层为箭头风场。其他小细节，上面已有论述，不必絮叨，但是怎么实现白化数据挖空，需要另外指出。对于contourf函数来说，使用levels参数，就可以控制数据的上下界限，例如：
kw={'levels':np.arange(70,105,5)}
cross_lon_90=ax.contourf(
rh[:,:,0].values, Y[:,:,0], Z[:,:,0],cmap='Blues',
zdir='x', offset=89.9, **kw)
cross_lat_50=ax.contourf(
X[:,0,:], rh[:,0, :].values, Z[:,0, :],cmap='Blues',
zdir='y', offset=50.1, **kw)
但是对于pcolormesh和plot_surface的栅格，这个命令就没用了，因为这两个没有levels参数。即便是指定上下限，也会出现额外填色，我们在pcolormesh章节已经谈到过了。上面这个青藏高原挖空，明显需要对原始数据进行处理。我们使用polygon的contains语句判断数据是否落在青藏高原内，然后将没有落在里面的变成np.nan，当然也可以使用regionmask。
from shapely.geometry import Point
from matplotlib.path import Path
shp=shpreader.Reader(r'C:\Users\lenovo\Desktop\DBATP_Polygon.shp')
geo_list=list(shp.geometries())
poly=geo_list[0]
上面将青藏高原作为polygon取出：
LON,LAT=np.meshgrid(lon,lat)
rh_850=rh.loc[850].values
for i in range(LON.shape[0]):
for j in range(LAT.shape[1]):
plon=LON[i][j]
plat=LAT[i][j]
if not poly.contains(Point(plon,plat)):
rh_850[i][j]=np.nan
else:
pass
然后将筛选后的数据重新绘图：
ax2=fig.add_axes([0,0,0.5,0.5])
a1=ax2.pcolormesh(LON,LAT,rh_850,shading='auto',cmap=cmap,norm=norm)
rh_flat=np.array(rh_850.flatten())
paths1=a1.get_paths()
polys1=concat(path.to_polygons() for path in paths1)
lc01=PolyCollection(polys1,cmap=cmap,edgecolor='none',closed=False,array=rh_flat)
ax.add_collection3d(lc01,zdir='z',zs=850)
text3d(ax, (105,50,100),'雲台書使', zdir="z",size=3,usetex=False,
ec="none", fc="k")
ax2.set_visible(False)
这样青藏高原以外的数据便消失了。contourf命令如法炮制便可。
腾讯云开发者
