使用Matplotlib & Cartopy繪制我國臺風路徑圖

夏天一到，沿海地區(qū)經常會遭到臺風的襲擾，可謂苦不堪言。

之前在公眾號做過一個關于我國1945~2015年歷史臺風統(tǒng)計的可視化展示，發(fā)現(xiàn)很多有趣的數(shù)據(jù)，比如說臺風登陸最多的城市是湛江。

大家可以去翻看歷史文章，附有完整代碼和數(shù)據(jù)，有興趣做些可視化探索。

大數(shù)據(jù)告訴你，臺風最喜歡在我國哪個省市登陸

這次的文章不研究臺風數(shù)據(jù)，而是嘗試用Python來繪制臺風路徑。

主要第三方庫

用到的主要工具包有pandas、numpy、matplotlib、cartopy、shapely，前三個庫大家可能都熟悉，下面介紹下后兩個庫的使用場景。

?

cartopy：基于matplotlib的python地理數(shù)據(jù)處理和可視化庫，本文會用來展示地圖
shapely: 是一個對幾何對象進行操作和分析的Python庫，本文用來處理點線數(shù)據(jù)

?

• cartopy文檔：https://scitools.org.uk/cartopy/docs/latest/
  
• shapely文檔：https://shapely.readthedocs.io/en/stable/

臺風路徑數(shù)據(jù)

本文用到的數(shù)據(jù)是我國2017年所有臺風路徑，包含了時間、經緯度、強度等關鍵信息。

由于數(shù)據(jù)來源網絡，沒法追溯真實性，僅供練習。

原始數(shù)據(jù)比較亂，我重新處理了方便使用：

可以看到共有7個字段：

?

臺風編號：我國熱帶氣旋編號 

日期：具體時間 

強度：0~9

緯度：單位0.1度 

經度：單位0.1度 

中心氣壓：hPa

中心最大風速：m/s

?

繪制地圖

臺風路徑需要在地圖上展示，那么如何獲取地圖呢？

方式有很多種，既可以用離線的GeoJson數(shù)據(jù)，也可以用JPG圖片，或者第三方庫提供的地圖。

我這里用的是cartopy內置的地圖數(shù)據(jù)，可以很方便的修改配置屬性。

首先導入本次會用到的所有庫：

# cartopy：用來獲取地圖
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# matplotlib：用來繪制圖表
import matplotlib.pyplot as plt
# shapely：用來處理點線數(shù)據(jù)
import shapely.geometry as sgeom
import warnings
import re
import numpy as np
import pandas as pd
warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = [u'SimHei']
plt.rcParams['axes.unicode_minus'] = False

獲取我國沿海區(qū)域地圖：

# 通過cartopy獲取底圖
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# 用經緯度對地圖區(qū)域進行截取，這里只展示我國沿海區(qū)域
ax.set_extent([85,170,-20,60], crs=ccrs.PlateCarree())

# 設置名稱
ax.set_title('2017年臺風路徑圖',fontsize=16)

# 設置地圖屬性，比如加載河流、海洋
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.LAKES, alpha=0.5)

# 展示地圖
plt.show()

處理數(shù)據(jù)并可視化

先加載數(shù)據(jù)：

typhoonData = pd.read_csv('typhoonData.csv')

這個數(shù)據(jù)集包含了30個臺風路徑，所以后面要分別進行可視化。

再對數(shù)據(jù)進行處理，依次提取單個臺風路徑及其經緯度。

# 先對臺風編號進行循環(huán)，提取單個臺風數(shù)據(jù)
for typhoonNumber in typhoonData['臺風編號'].unique():
    typhoon = typhoonData[typhoonData['臺風編號']==typhoonNumber]
    # 再對單個臺風數(shù)據(jù)進行處理，提取經緯度
    for typhoonPoint in np.arange(len(typhoon)-1):
        lat_1 = typhoon.iloc[typhoonPoint,3]/10
        lon_1 = typhoon.iloc[typhoonPoint,4]/10
        lat_2 = typhoon.iloc[typhoonPoint+1,3]/10
        lon_2 = typhoon.iloc[typhoonPoint+1,4]/10
        point_1 = lon_1,lat_1
        point_2 = lon_2,lat_2
        # 最后可視化
        ax.add_geometries([sgeom.LineString([point_1,point_2])],crs=ccrs.PlateCarree(),edgecolor='red')
        
    
# 展示圖像
plt.show()

能看到所有臺風路徑都被描繪出來了。

但這里沒有區(qū)別顯示臺風強度，一般是在.add_geometries()方法中添加參數(shù)調整。

有兩種方式：

• 用顏色區(qū)別：不同顏色代表不同強度，參數(shù)-edgecolor
• 用線條粗細區(qū)別：越粗則強度越高，參數(shù)-linewidth

1. 顏色區(qū)分

# 按強度區(qū)分顏色
def get_color(level):
    if level in (0,1):
        color='#ff00ff'
    elif level in (2,3):
        color='#ff00cc'
    elif level in (4,5):
        color='#ff0066'
    elif level in (6,7):
        color='#ff0033'
    elif level in (8,9):
        color='#ccff00'
    return color

# 先對臺風編號進行循環(huán)，提取單個臺風數(shù)據(jù)
for typhoonNumber in typhoonData['臺風編號'].unique():
    typhoon = typhoonData[typhoonData['臺風編號']==typhoonNumber]
    # 再對單個臺風數(shù)據(jù)進行處理，提取經緯度
    for typhoonPoint in np.arange(len(typhoon)-1):
        lat_1 = typhoon.iloc[typhoonPoint,3]/10
        lon_1 = typhoon.iloc[typhoonPoint,4]/10
        lat_2 = typhoon.iloc[typhoonPoint+1,3]/10
        lon_2 = typhoon.iloc[typhoonPoint+1,4]/10
        point_1 = lon_1,lat_1
        point_2 = lon_2,lat_2
        # 最后可視化，添加顏色參數(shù)
        ax.add_geometries([sgeom.LineString([point_1,point_2])],crs=ccrs.PlateCarree(),edgecolor=get_color(typhoon.iloc[typhoonPoint,2]))
        
# 展示圖像
plt.show()

2. 線條粗細區(qū)分

# 先對臺風編號進行循環(huán)，提取單個臺風數(shù)據(jù)
for typhoonNumber in typhoonData['臺風編號'].unique():
    typhoon = typhoonData[typhoonData['臺風編號']==typhoonNumber]
    # 再對單個臺風數(shù)據(jù)進行處理，提取經緯度
    for typhoonPoint in np.arange(len(typhoon)-1):
        lat_1 = typhoon.iloc[typhoonPoint,3]/10
        lon_1 = typhoon.iloc[typhoonPoint,4]/10
        lat_2 = typhoon.iloc[typhoonPoint+1,3]/10
        lon_2 = typhoon.iloc[typhoonPoint+1,4]/10
        point_1 = lon_1,lat_1
        point_2 = lon_2,lat_2
        # 最后可視化
        ax.add_geometries([sgeom.LineString([point_1,point_2])],crs=ccrs.PlateCarree(),linewidth = typhoon.iloc[typhoonPoint,2],edgecolor='red')
        
# 展示圖像
plt.show()

最后

上文用比較簡單的方式繪制了臺風路徑圖，大家可以嘗試換個三維地圖，或者用動態(tài)顯示臺風走勢...

玩法挺多的，趕緊嘗試嘗試吧。

?

后臺回復：t，獲取數(shù)據(jù)集和完整代碼

?