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

上一篇文章中,利用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添加底图的过程还是有些麻烦,主要包括:

  1. 下载底图,可以从NASA VisibleEarth上下载;

  2. 设定Cartopy底图的环境变量

     os.environ["CARTOPY_USER_BACKGROUNDS"]="存放底图的文件路径"
    
  3. 编写底图的json文件,需要的信息包括

    • "__comment__"
    • "__source__"
    • "__projection__"
    • "__source__"
  4. 绘制底图,例如选择中等分辨率的蓝色弹珠底图(其实Blue Marble提供的分辨率相当高,可以根据自己的需要重新调整):

    ax.background_img(name='blue_marble', resolution='medium')

我选择了北半球绿色最多的7月份,最后呈现的效果,如下所示:
Blue Marble用作Cartpy底图

更丰富的台风信息呈现

台风信息

在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(),例如中心气压的绘制代码如下:

1
2
3
4
5
# 显示中心气压的子图
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()函数可以在图中添加注释。

对每个路径点,在路径点旁边注释该位置的台风等级:

1
2
3
4
5
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的帧率等,代码非常简单:

1
2
3
4
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图片较大,加载速度较慢):
台风发展动图

未完待续

本来其实没有打算做这份工作,直接开始做些数据统计分析的,但是在网上看到了类似的工作,觉得效果不错,我也就做了一下重现,并且简化了其中关于路径点和线的绘制。


参考文献

  1. https://earthobservatory.nasa.gov/Features/BlueMarble/ Blue Marble: Next Generation
  2. http://www.cnblogs.com/vamei/archive/2012/11/07/2758006.html 飓风“桑迪”路径图的制作
  3. http://earthpy.org/cartopy_backgroung.html Background image in cartopy
  4. https://pillow.readthedocs.io/en/4.2.x/ Pillow documentation
张da统帅 wechat
扫码订阅我的微信公众号