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

发表于 分类于 阅读次数： Disqus：

在上一篇文章中，我们实现了台风最佳路径数据的读取和预处理。本篇我们将利用Cartopy和matplotlib显示地图和台风路径，在这里我们重点参考了Cartopy的例子。

Python库介绍

matplotlib绘制图表

matplotlib可以说是Python中最著名的绘图库，看这个名字就可以知道matplotlib提供了一套与Matlab中非常类似的绘图函数，可以非常方便地绘制各种精美的图表。此外，还有许多基于matplotlib的库可供使用

Basemap是matplotlib的一个自包，顾名思义可知是一个用于绘制地图的工具，Basemap提供了许多绘制地图相关的功能。但是，如果是采用Anaconda安装的PPython，需要自己安装Basemap。我用的是Python 3.6，结果直接采用conda安装Basemap失败，折腾了好久（这就是另外一个故事了……），最后采用离线安装的方式解决了，简而言之就是：

  1. 根据Basemap安装说明安装’Basemap’依赖的GEOS、PROJ4等包；
  2. 从unofficial python Extension Packages下载Basemap的whl文件；
  3. 采用pip在本地安装Basemap。

其实我在自己电脑上没有这么复杂也安装成功了，实在是不可思议。

言归正传，Basemap的功能虽然强大，但是应用起来略显繁琐。对于我这种一知半解的人，用一些更加方便的GIS包也就是情理之中的事情了。

Cartopy GIS包

Cartopy功能强大，方便快捷，实属居家旅行之必备，啧啧……
采用conda install -c scitools cartopy瞬间安装好，然后就可以用了。

Cartopy 绘制地图

底图显示

Cartopy本身只有一个背景底图，采用stock_img()即可绘制：

1
2
3
4
5
6
7
8
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

f = plt.figure(figsize=(16,9))
ax = plt.axes(projection=ccrs.Robinson())
ax.stock_img()

plt.show()

即可显示世界地图：
世界地图

此外还可以画出世界的海岸线轮廓：

1
2
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()

显示效果如下：costlines

但是！这些还不能满足我们的需求，我们需要显示台风经过的行政区划。这就需要shapefile了。

shapefile

简介&来源选择

shapefile，全称是ESRI Shapefile，是由全世界最大的GIS厂商、大名鼎鼎的ArcGIS的出品方ESRI开发的空间数据格式，已经称为GIS的开放标准。Shapefile用于描述几何体对象，是一种矢量图形格式。
关于中国行政区划的shapefile，有多种渠道可以获得：

  1. Cartopy自带的 shpreader.natural_earth()函数可以自动从Natural Earth上下载不同经度、不同类型的shapefile,当然也可以自己手动下载然后在本地加载；
  2. Global Administative Aream的GADM数据库提供全世界各个国家的数据；
  3. 百度网盘里可以找到原来中国国家地理信息中心发布的中国的shapefile，信息相对较老；
  4. 国家地理信息中心现在已经不再提供公开下载，但是如果需要可以国家基础地理信息中心信息服务行政许可受理大厅购买。

总结一下：方式1,2,3都是公开免费的，方式4需要花钱；1，2由于是国外提供的，在国家领土和边界处理上存在问题。综上，我选择了方式3。

shapefile导入与显示

利用Cartopy的cartopy.io.shapereader可以读取shapefile。在利用coastlines()已经画好的轮廓图上叠加shapefile的图层，就可以得到中国的行政区划。我们采用的是精确到县级的行政区划，加载的时间比较长，显示的范围为东经100-160度，赤道到北纬50度。
在Jupyter Notebook中，最好用%matplotlib qt5将图显示到GUI中，方便放缩。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import cartopy.io.shapereader as shpreader
file_name = 'F:/Git/ZGis/chinashp/县级bou4_4m/BOUNT_poly.shp'
reader = shpreader.Reader(file_name)

adm1_shapes = list(reader.geometries())

ax = plt.axes(projection=ccrs.PlateCarree())

ax.coastlines(resolution='10m')

ax.add_geometries(adm1_shapes, ccrs.PlateCarree(),
                  edgecolor='black', facecolor='gray', alpha=0.5)

ax.set_extent([100, 160, 0, 50], ccrs.PlateCarree())

plt.show()

显示效果如下：
中国行政区划（县市级）

绘制台风路径和影响县市

在上一篇文章中，可以利用set_index()将台风编号设置为索引，然后选择某一个编号对应的BST数据，根据其经纬度坐标，利用plot()即可画出台风路径。为了描述其影响区域，我们利用shapely.geometry的buffer()函数可以画出一个包络。但是需要注意的是，整理的包络实际上是shapely通过调整“画笔”的粗细得出的，并不是实际的距离，因此更多的是作为示意。

1
2
3
4
track = sgeom.LineString(zip(lons, lats))
track_buffer = track.buffer(1)
plt.plot(lons, lats, color='black', transform=ccrs.PlateCarree())
ax.add_geometries([track_buffer], ccrs.PlateCarree(), facecolor=facecolor, alpha=0.5)

继续使用Cartopy的geometry()函数提取县市的边界，用intersects()判断其与台风路径或者影响区域的相交情况,如果与台风路径相交则设置为红色，与影响区域相交则为橙色。

1
2
3
4
5
6
for county in adm1_shapes:
    facecolor = [0.9375, 0.9375, 0.859375]
    if county.intersects(track):
        facecolor = 'red'
    elif county.intersects(track_buffer):
        facecolor = '#FF7E00'

我们选择1508号台风为例，最后显示结果如下：
1508影响范围

插曲：中文图例的显示

matplotlib的中文图例默认会显示为乱码，需要进行一定的设置才可以正常显示。我用了网上的两种方法，都可以：

  1. 利用Fontmanger临时设置字体，有点像matlab的做法
1
2
3
font = FontProperties(fname=r"c:\windows\fonts\NotoSansCJKsc-Regular.otf", size=14)
plt.legend([direct_hit, within_2_deg], labels,
           loc='lower left', bbox_to_anchor=(0.025, -0.1), fancybox=True, prop=font)
  1. 在matplotlib的字体里增加支持中文的字体文件（同时要修改配置文件），路径为C:\Anaconda3\Lib\site-packages\matplotlib\mpl-data\

未完待续

到这里，我只是完成了选择一个台风编号，显示其影响的区域和县市。这只是非常初步的工作，接下来：

  • 考虑将人口和GDP的数据加进来，判断其影响的人口和经济，然后根据实际造成的损失，估计一个模型，用于未来损失的预测；
  • 统计历史上台风的强度、寿命、时间等并进行可视化;
  • 对了，这里面有没有什么幂律之类的存在吗？

一点感受

第一次写略微有点“技术含量”的博客，感觉还是收获颇丰，不管是Jupyter Notebook、Markdown的使用，还是思路的整理、语言的组织，当然，更让我看到了进步的空间。

参考文献

  1. http://scitools.org.uk/cartopy/docs/latest/examples/hurricane_katrina.html Cartopy hurricane_katrina example
  2. http://www.naturalearthdata.com Natural Earth
  3. http://desktop.arcgis.com/zh-cn/arcmap/10.3/manage-data/shapefiles/what-is-a-shapefile.htm 什么是shapefile？
  4. https://matplotlib.org/users/customizing.html Customizing matplotlib