使用Matplotlib & Cartopy繪制我國(guó)臺(tái)風(fēng)路徑圖
夏天一到，沿海地區(qū)經(jīng)常會(huì)遭到臺(tái)風(fēng)的襲擾，可謂苦不堪言。
之前在公眾號(hào)做過(guò)一個(gè)關(guān)于我國(guó)1945~2015年歷史臺(tái)風(fēng)統(tǒng)計(jì)的可視化展示，發(fā)現(xiàn)很多有趣的數(shù)據(jù)，比如說(shuō)臺(tái)風(fēng)登陸最多的城市是湛江。
大家可以去翻看歷史文章，附有完整代碼和數(shù)據(jù)，有興趣做些可視化探索。
大數(shù)據(jù)告訴你，臺(tái)風(fēng)最喜歡在我國(guó)哪個(gè)省市登陸
這次的文章不研究臺(tái)風(fēng)數(shù)據(jù)，而是嘗試用Python來(lái)繪制臺(tái)風(fēng)路徑。
主要第三方庫(kù)
用到的主要工具包有pandas、numpy、matplotlib、cartopy、shapely，前三個(gè)庫(kù)大家可能都熟悉，下面介紹下后兩個(gè)庫(kù)的使用場(chǎng)景。
??
cartopy：基于matplotlib的python地理數(shù)據(jù)處理和可視化庫(kù)，本文會(huì)用來(lái)展示地圖shapely: 是一個(gè)對(duì)幾何對(duì)象進(jìn)行操作和分析的Python庫(kù)，本文用來(lái)處理點(diǎn)線數(shù)據(jù)
cartopy文檔：https://scitools.org.uk/cartopy/docs/latest/ shapely文檔：https://shapely.readthedocs.io/en/stable/
臺(tái)風(fēng)路徑數(shù)據(jù)
本文用到的數(shù)據(jù)是我國(guó)2017年所有臺(tái)風(fēng)路徑，包含了時(shí)間、經(jīng)緯度、強(qiáng)度等關(guān)鍵信息。
由于數(shù)據(jù)來(lái)源網(wǎng)絡(luò)，沒(méi)法追溯真實(shí)性，僅供練習(xí)。
原始數(shù)據(jù)比較亂，我重新處理了方便使用：
可以看到共有7個(gè)字段：
?臺(tái)風(fēng)編號(hào)：我國(guó)熱帶氣旋編號(hào)
日期：具體時(shí)間
強(qiáng)度：0~9
緯度：?jiǎn)挝?.1度
經(jīng)度：?jiǎn)挝?.1度
中心氣壓：hPa
中心最大風(fēng)速：m/s
?
繪制地圖
臺(tái)風(fēng)路徑需要在地圖上展示，那么如何獲取地圖呢？
方式有很多種，既可以用離線的GeoJson數(shù)據(jù)，也可以用JPG圖片，或者第三方庫(kù)提供的地圖。
我這里用的是cartopy內(nèi)置的地圖數(shù)據(jù)，可以很方便的修改配置屬性。
首先導(dǎo)入本次會(huì)用到的所有庫(kù)：
# cartopy：用來(lái)獲取地圖
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# matplotlib：用來(lái)繪制圖表
import matplotlib.pyplot as plt
# shapely：用來(lái)處理點(diǎn)線數(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
獲取我國(guó)沿海區(qū)域地圖：
# 通過(guò)cartopy獲取底圖
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
# 用經(jīng)緯度對(duì)地圖區(qū)域進(jìn)行截取，這里只展示我國(guó)沿海區(qū)域
ax.set_extent([85,170,-20,60], crs=ccrs.PlateCarree())
# 設(shè)置名稱
ax.set_title('2017年臺(tái)風(fēng)路徑圖',fontsize=16)
# 設(shè)置地圖屬性，比如加載河流、海洋
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')
這個(gè)數(shù)據(jù)集包含了30個(gè)臺(tái)風(fēng)路徑，所以后面要分別進(jìn)行可視化。
再對(duì)數(shù)據(jù)進(jìn)行處理，依次提取單個(gè)臺(tái)風(fēng)路徑及其經(jīng)緯度。
# 先對(duì)臺(tái)風(fēng)編號(hào)進(jìn)行循環(huán)，提取單個(gè)臺(tái)風(fēng)數(shù)據(jù)
for typhoonNumber in typhoonData['臺(tái)風(fēng)編號(hào)'].unique():
typhoon = typhoonData[typhoonData['臺(tái)風(fēng)編號(hào)']==typhoonNumber]
# 再對(duì)單個(gè)臺(tái)風(fēng)數(shù)據(jù)進(jìn)行處理，提取經(jīng)緯度
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()
能看到所有臺(tái)風(fēng)路徑都被描繪出來(lái)了。
但這里沒(méi)有區(qū)別顯示臺(tái)風(fēng)強(qiáng)度，一般是在.add_geometries()方法中添加參數(shù)調(diào)整。
有兩種方式：
用顏色區(qū)別：不同顏色代表不同強(qiáng)度，參數(shù)-edgecolor 用線條粗細(xì)區(qū)別：越粗則強(qiáng)度越高，參數(shù)-linewidth
顏色區(qū)分
# 按強(qiáng)度區(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
# 先對(duì)臺(tái)風(fēng)編號(hào)進(jìn)行循環(huán)，提取單個(gè)臺(tái)風(fēng)數(shù)據(jù)
for typhoonNumber in typhoonData['臺(tái)風(fēng)編號(hào)'].unique():
typhoon = typhoonData[typhoonData['臺(tái)風(fēng)編號(hào)']==typhoonNumber]
# 再對(duì)單個(gè)臺(tái)風(fēng)數(shù)據(jù)進(jìn)行處理，提取經(jīng)緯度
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()
線條粗細(xì)區(qū)分
# 先對(duì)臺(tái)風(fēng)編號(hào)進(jìn)行循環(huán)，提取單個(gè)臺(tái)風(fēng)數(shù)據(jù)
for typhoonNumber in typhoonData['臺(tái)風(fēng)編號(hào)'].unique():
typhoon = typhoonData[typhoonData['臺(tái)風(fēng)編號(hào)']==typhoonNumber]
# 再對(duì)單個(gè)臺(tái)風(fēng)數(shù)據(jù)進(jìn)行處理，提取經(jīng)緯度
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()
最后
上文用比較簡(jiǎn)單的方式繪制了臺(tái)風(fēng)路徑圖，大家可以嘗試換個(gè)三維地圖，或者用動(dòng)態(tài)顯示臺(tái)風(fēng)走勢(shì)...
玩法挺多的，趕緊嘗試嘗試吧。
?后臺(tái)回復(fù)：t，獲取數(shù)據(jù)集和完整代碼
?
加入知識(shí)星球【我們談?wù)摂?shù)據(jù)科學(xué)】
400+小伙伴一起學(xué)習(xí)！
· 推薦閱讀 ·
數(shù)據(jù)分析最有用的25個(gè)Matplotlib圖表