Cartopy自定义底图
Cartopy自定义底图
蓝色弹珠
前几天微信的启动页面做了短时间的“变脸”：从NASA的蓝色弹珠(Blue Marble) 换成了风云四号拍摄的高清图像。
蓝色弹珠可谓是妇孺皆知的地球照片，是阿波罗17号太空船船员于1982年12月7日所拍摄。2000年，NASA利用卫星遥感数据以“蓝色弹珠”为名发布了全新的地球全景影像。2005年，NASA更是发布了“Blue Marble: Next Generation ”,包含总共十二个月份拍摄的地球卫星影响，可以显示地球在这十二个月内不同季节产生的变化，包含冰层的扩张与缩减、植被的疏密的改变、季节性干旱、南北半球因季节不同而形成的差异等。此后，蓝色弹珠又多次更新，并且可以在NASA的网站上选择合适的版本下载。
2005年的“蓝色弹珠：下一代”
“蓝色弹珠”用作底图
从前面的文章可以看到，Cartopy的底图实在是乏善可陈，默认的底图只有stock_img()中的一张。幸运的是，从Cartopy 0.15开始，可以利用background_img()自定义底图，而Blue Marble就是一个合适的选择。
不过，Cartopy添加底图的过程还是有些麻烦，主要包括：
-
下载底图，可以从NASA VisibleEarth上下载；
-
设定Cartopy底图的环境变量
os.environ["CARTOPY_USER_BACKGROUNDS"]="存放底图的文件路径" -
编写底图的json文件，需要的信息包括
"__comment__""__source__""__projection__""__source__"
-
绘制底图，例如选择中等分辨率的蓝色弹珠底图（其实Blue Marble提供的分辨率相当高，可以根据自己的需要重新调整）：
ax.background_img(name='blue_marble', resolution='medium')
我选择了北半球绿色最多的7月份，最后呈现的效果，如下所示：
更丰富的台风信息呈现
台风信息
在CMA提供的台风BST数据中，除了台风的经纬度以外，还有记录的时刻、相应时刻的台风等级、中心最低气压和2分钟平均近中心最大风速。
在作图时，可以将这些信息都呈现出来：
- 根据台风等级的不同，绘制相应台风路径的颜色；
- 整个过程中台风中心最低气压的变化情况；
- 整个过程中台风近中心最大风速的变化情况；
- 在每个记录时刻，台风的相应信息
根据信息作图
台风等级与路径颜色
热带气旋可以分为弱于热带风暴到超强台风7个等级，此外还有变性，用0~6和9表示，对这9个等级建立颜色的字典，用于后续调用。
level_dict = {0: '弱于热带风暴', 1: '热带低压', 2: '热带风暴', 3: '强热带风暴', 4: '台风', 5: '强台风', 6: '超强台风', 9: '变性'}
路径点&线作图
对于k时刻的台风，需要将该时刻及以前的k个路径点和（k-1）条线段按照强度对应的颜色画出来。
首先，根据k个点的强度值，设定颜色的循环顺序：
color_seq = [color_dict[lev] for lev in level] plt.rc('axes', prop_cycle=(cycler('color', color_seq)))
对于折线图，可以按照设定好的循环直接画图，但是对于散点图，如果利用scatter()一次性画出所有路径点，颜色循环是无效的，需要直接传入color参数:
其他信息子图显示
台风路径所占的经纬度范围是已知的，适当多画一些范围，在底图的右下角分别绘制中心气压和最大风速的底图，调用的函数为add_axes()，例如中心气压的绘制代码如下:
显示中心气压的子图
ax_pressure = fig.add_axes([0.7, 0.45, 0.15, 0.1])
ax_pressure.set_title('中心最低气压(hPa)', fontsize=10, fontproperties=font)
ax_pressure.plot(range(1, N + 1), single_track.pres)
ax_pressure.axvline(i, color='red')
文字注释
matplotlib中的annotate()函数可以在图中添加注释。
对每个路径点，在路径点旁边注释该位置的台风等级:
ax.annotate(level_dict[temp.level], xy=(temp.lons, temp.lats), xycoords= transform,
xytext=(-5, -30), textcoords='offset points', ha='right', va='bottom',
bbox=dict(boxstyle='round,pad=0.5', fc=color_dict[temp.level], alpha=0.8),
fontproperties=font,
arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0'))
此外，还要在图的左上角显示台风的相关信息，例如台风名称、数据来源、图像对应的记录时刻等：
ax.annotate('台风名称: %s\n\n时间：%s\n数据来源：CMA\n张da统帅制作' % (temp['name'], temp.time),
xy=(0, 1), xytext=(12, -12), va='top', ha='left', fontproperties=font,
xycoords='axes fraction', textcoords='offset points')
这样，得到某一时刻的台风路径图如下所示：
image.png
利用pillow绘制台风动图
按照以上步骤，把一个台风对应的所有时刻的台风路径图片绘制出来，按照绘制的顺序进行命名，并保存在一个文件夹下，把这个文件夹下所有的路径按照顺序读取：
os.chdir(path)
imgFiles = [fn for fn in os.listdir('.') if fn.endswith('.png')]
imgFiles.sort(key=lambda x:int(x[:-4]))
图片的读取和GIF图片的制作，可以用Pillow`工具包完成，Anaconda安装的时候已经自带，可以直接是用。在制作gif时，可以设定gif的帧率等，代码非常简单：
images = [Image.open(fn) for fn in imgFiles]
im = images[0]
filename = 'test.gif'
im.save(fp=filename, format='gif', save_all=True, append_images=images[1:], duration=500)
</pre>
最后的效果如下所示(GIF图片较大，加载速度较慢)：
参考文献
- https://earthobservatory.nasa.gov/Features/BlueMarble/ Blue Marble: Next Generation
- http://www.cnblogs.com/vamei/archive/2012/11/07/2758006.html 飓风“桑迪”路径图的制作
- http://earthpy.org/cartopy_backgroung.html Background image in cartopy
- https://pillow.readthedocs.io/en/4.2.x/ Pillow documentation
