Python气象绘图教程（十五）—Cartopy_5
本节提要：仿制中央气象台气象服务图片、关于cartopy里的投影与转换、cartopy中extent与boundary。
一、仿制中央气象台图片
从鄙人高三填报了南信的志愿开始，就一直持续的关注中央气象台，也算是一个老看客了。期间还经历了比较大型的改版（我感觉改版后变丑了O(∩_∩)O哈哈~），我觉得央台的气象服务图片做得比较好，基本上以填色图和散点图为主，辅以少量雷达图，绝大多数图片都能很好的向各层次的读者展现天气预报与实况。
今天就接着上次谈论的图例和色条，来谈谈py下仿制央台图片流程。
我仿制的图片如下：
A、仿制的高温图
恩施想要这样的高温还是比较困难的（恩施算得上矮高原了，鄂西凉都），所以修改成了地质灾害的预报。
首先涉及到资料的问题，地质灾害不在常规预报里，但是气象局也必须发这项预警。这里只能用事先做好的实验数据，预报在清江两岸有一定的山洪泥石流风险。使用的仍然是不均匀的站点资料，所以先要将站点资料网格化，变为格点资料后再用等值线填色的方式画出危险区域。最后，通过前面几节提到的添加图例的方法，完善图形。关键代码如下：
olon=np.linspace(108,111,90)
olat=np.linspace(29,32,90)
olon,olat=np.meshgrid(olon,olat)
func=Rbf(lon,lat,danger,function='linear')
danger_new=func(olon,olat)
colorlevel=[0,60,80,100]#危险等级
colordict=['white','orange','red']#颜色列表
danger_map=mcolors.ListedColormap(colordict)#产生颜色映射
norm=mcolors.BoundaryNorm(colorlevel,danger_map.N)#生成索引
ax.contourf(olon,olat,danger_new,levels=colorlevel,cmap=danger_map,norm=norm)#填色图
clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000)#白化
ax.set_title('地质灾害风险（sample data）',fontsize=12)
danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")
low_danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="orange")
labels = ['泥石流高风险','内涝积水风险']
ax.legend([danger, low_danger], labels,fontsize=7,frameon=False,title='图例',
loc='lower left', bbox_to_anchor=(-0.01, 0), fancybox=False)
B、降水量图的仿制
使用的是最开始实验自定义colorbar时的那张图的数据，但是当时用的是色条来表示降水量，这次我们用图例的方式表示降水量，前面的步骤和A中的类似。最后添加图例时候，我用的是逐个添加，实际上可以用循环的方式添加。代码如下：
filename=r'C:\Users\lenovo\Desktop\ex1.xlsx'
df=pd.read_excel(filename)
lon=df['lon']#经度
lat=df['lat']#纬度
olon=np.linspace(108,111,90)
olat=np.linspace(29,32,90)
olon,olat=np.meshgrid(olon,olat)
rain=df['rain_new']
func=Rbf(lon,lat,rain,function='linear')
rain_new=func(olon,olat)
colorlevel=[0.1,10.0,25.0,50.0,100.0,250.0,500.0]#雨量等级
colordict=['#A6F28F','#3DBA3D','#61BBFF','#0000FF','#FA00FA','#800040']#颜色列表
rain_map=mcolors.ListedColormap(colordict)#产生颜色映射
norm=mcolors.BoundaryNorm(colorlevel,rain_map.N)#生成索引
cs= ax.contourf(olon,olat,rain_new,levels=colorlevel,cmap=rain_map,norm=norm)
clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000)#白化
############添加图例#######################
larger1 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#A6F28F")
larger2 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#3DBA3D")
larger3 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#61BBFF")
larger4 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#0000FF")
larger5 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#FA00FA")
labels = ['0~10','10~25','25~50','50~100','100~250']
ax.legend([larger1,larger2,larger3,larger4,larger5], labels,fontsize=12,frameon=False,title='图例（mm）',
loc='lower left', bbox_to_anchor=(-0.01, -0.03), fancybox=False)
其中，mpatches.Rectangle这个语句就是添加图例上的色条，rectangle就是矩形的意思，延伸的，你也可以使用其他形状添加图例（包括自定义形状）。如何提供更丰富的标记，只能烦请各位在官网文档上查找了。
二、Cartopy里的投影与转换
我在两个月前经常碰到这个问题，有两个关于投影的—crs、transform。前一个设定投影方式，后一个涉及投影与数据转换。
首先说crs，这个是GeoAxes的基础，只有在projection设置了投影之后，才能添加地图。如果没有一步转化，有一些专属于GeoAxes的命令就无法生效：
'AxesSubplot' object has no attribute 'add_feature'比较典型的如add_feature、add_geometries添加地图和地理特征的命令不能使用。
crs还控制着数据绘图与边界的裁剪，比如set_extent(crs=ccrs.PlateCarree()) ，就使裁剪的方式按照PlateCarree()的方式进行边界的裁剪，一个经典的案例即兰勃脱下的使用extent进行裁剪，那么将会在地图中出现范围之外的地图：
如果主地图投影为默认投影方式（PlateCarree），那么，你在后续的绘制中是不需要进行投影转换的。
但是，在主地图不是默认投影方式时，需要进行转换。比如：
ax=fig.add_subplot(1,1,1,projection=ccrs.LambertConformal(central_longitude=110, central_latitude=35))
ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r')这种情况下，将无法绘制出数据：
这是因为，数据默认的格式为PlateCarree下的投影，需要transform的介入：
ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r',transform=ccrs.PlateCarree())
transform=PlateCarree()的意思是，这个数据原来是PlateCarree()投影下的，转化为主图的格式后（即LambertConformal），再绘制等值线。
主图ax的投影方式 | 数据的投影方式 | 转换命令 |
|---|---|---|
LambertConformal | PlateCarree | transform=PlateCarree |
LambertConformal | Mercator | transform=Mercator |
PlateCarree | Mercator | transform=Mercator |
三、通过set_boundary和set_extent绘制扇形图
cartopy官网实例有一个Custom Boundary Shape，通过该例子可以裁剪边框形状。现在尝试裁剪东亚地区的兰勃脱投影。
import numpy as np
from matplotlib.path import Path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat
vertices = [(65, 10), (65, 60), (145, 60), (145, 10), (65, 10)]#五个点，但是首尾是一样的，以连接为封闭的四边形
boundary = Path(vertices)#边界形状
fig=plt.figure(figsize=(5,5),dpi=500)
ax=plt.axes([0,0,1,1],projection=ccrs.LambertConformal(central_longitude=110, central_latitude=35))
ax.add_feature(cfeat.OCEAN)#添加海洋
ax.add_feature(cfeat.LAND)#添加陆地
ax.set_boundary(boundary, transform=ccrs.PlateCarree())#裁剪边界
ax.set_extent([65,145,0,60])#裁剪地图范围
这时，可能因为投影转换的原因，出现了几个明显的折痕，于是，进一步提高path上点的数量：
latmin = 10
latmax = 60
lonmin = 70
lonmax = 140
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
[(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())
ax.set_extent(extent)
此时，因为vertices中的边界点的数量增多，边界连接更加圆滑。请注意，在这种方式下裁剪了边界，一定要使用set_extent来裁剪地图大小，不然就会出现地图蜷缩在图里这种情况：
在cartopy=0.17中，不能使用draw_labels=True来为除PlateCarree、Mercator之外的投影添加经纬标签，不过据说在0.18版本中已经优化，读者可以试试。通过前面提到的数据转换，绘制了一张气温图：
ax.contourf(data_lon,data_lat,temp,levels=np.arange(-40,40),cmap='Spectral_r',transform=ccrs.PlateCarree())
腾讯云开发者
