基于Python的GIS分析之台风路径可视化（三）

发表于 2017-09-26 分类于闭门造车阅读次数： Disqus：

在上一篇文章中，利用Cartopy和matplotlib实现了台风路径的可视化。为了呈现更加丰富的信息，对单张路径图显示更多的信息，同时制作了单个台风的GIF动图。

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月份，最后呈现的效果，如下所示：
Blue Marble用作Cartopy底图

更丰富的台风信息呈现

台风信息

在CMA提供的台风BST数据中，除了台风的经纬度以外，还有记录的时刻、相应时刻的台风等级、中心最低气压和2分钟平均近中心最大风速。

在作图时，可以将这些信息都呈现出来：

根据台风等级的不同，绘制相应台风路径的颜色；
整个过程中台风中心最低气压的变化情况；
整个过程中台风近中心最大风速的变化情况；
在每个记录时刻，台风的相应信息

根据信息作图

台风等级与路径颜色

热带气旋可以分为弱于热带风暴到超强台风7个等级，此外还有变性，用0~6和9表示，对这9个等级建立颜色的字典，用于后续调用。

1
2

color_dict = {0: '#7fffd4', 1: '#008000', 2: '#0000ff', 3: '#ffff00', 4: '#ffa500', 5: '#ff4500', 6: '#800080', 9: '#808080'}
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参数:

1 2	plt.scatter(lon, lat, color=color_seq, transform=ccrs.PlateCarree()) plt.plot((lon[:-1], lon[1:]), (lat[:-1], lat[1:]), transform=ccrs.PlateCarree())

其他信息子图显示

台风路径所占的经纬度范围是已知的，适当多画一些范围，在底图的右下角分别绘制中心气压和最大风速的底图，调用的函数为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'))

此外，还要在图的左上角显示台风的相关信息，例如台风名称、数据来源、图像对应的记录时刻等：

1
2
3

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')

这样，得到某一时刻的台风路径图如下所示：
台风信息展示

利用pillow绘制台风动图

按照以上步骤，把一个台风对应的所有时刻的台风路径图片绘制出来，按照绘制的顺序进行命名，并保存在一个文件夹下，把这个文件夹下所有的路径按照顺序读取：

1
2
3

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)

最后的效果如下所示(GIF图片较大，加载速度较慢)：
台风发展动图