python 绘制风向频率玫瑰图

研究台风路径和影响，除了直接绘制台风路径，我们还往往想知道研究的台风在海上各个经纬度的频数分布，并直观展示在地图上。

之前和大家交流过如何从台风路径数据集里提取指定条件的台风（用pandas库提取IBTrACS中特定条件的热带气旋最佳路径数据），下面就在这个工作的基础上，进一步来绘制台风路径数据。

这里统计频数的经纬度网格分辨率采用1°×1°，因为这样索引最方便（因为整型可以直接作为数组的索引值），绘制出整个西北太平洋台风的频数图像也不会过于粗糙。

为了不重复统计，每个台风所经过的经纬度网格（1°×1°）只统计一次，即使台风在同个经纬度网格内逗留多个时次或多次进入同个网格，也不会重复统计。

在开始之前，把需要研究的台风编号（JTWC的编号）取出，方便后续用这个编号（也就是每个台风的唯一标识）来进行DataFrame的筛选和提取。

由于本人之前已经进行了相关工作，把筛选下来的台风编号放到了csv文件里了，所以这里代码是直接从csv文件里读取这个列表的，但务必确保这个列表里没有重复的台风，不然就会重复统计了。

有了台风列表和台风路径数据集下面就可以统计每个经纬度网格里台风经过的频数

import pandas as pd
import numpy as np
import os
import csv
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import cartopy.feature as cfeature
import csv
import cmaps as cm

# 读取台风数据集
TCdata = pd.read_csv('IBTrACS-JTWC.csv',index_col=0)

TC_list=[]
# 读取台风编号列表
with open('WNP_TC_1979-2020_USA_ID.csv') as f:
    csv_file = csv.reader(f)
    for row in csv_file :
        TC_list.append(row[0])

# 设定研究海区的经纬度范围
region =[100,180,0,60]   
lon = np.arange(region[0],region[1]+1,1)
lat  =np.arange(region[2],region[3]+1,1)

# 新建数组，用于存放每个格点（1°×1°）的台风频数 
arr_N = np.zeros((len(lon),len(lat)),dtype=int)

# 遍历台风列表
for ID in TC_list:

    data = TCdata[TCdata['USA_ATCF_ID']==ID]
    data = data.copy()

    # 把经纬度变为整型
    data['lon_i']=data['USA_LON'].astype('int')
    data['lat_i']=data['USA_LAT'].astype('int')
   
    # 创建用于存放整型经纬度的df
    data_lonlat=data[['lon_i','lat_i']]

    # 去掉重复值，避免重复统计
    data_lonlat = data_lonlat.drop_duplicates()

    # 去掉不在研究海区范围内的台风
    data_lonlat=  data_lonlat[[ i >= region[0] and  i <= region[1] for i in data_lonlat['lon_i']]]
    data_lonlat=  data_lonlat[[ i >= region[2] and  i <= region[3] for i in data_lonlat['lat_i']]]
    
    # 遍历每行经纬度，即台风经过的点
    for rows in data_lonlat.itertuples():
        # 用经纬度来索引数组，数组的index是从0开始的，所以要减去研究海区的起始经度region[0]和起始维度region[2]
        # 让索引出来的数组的值（即该位置的台风频数）+1
        arr_N[rows[1]-region[0],rows[2]-region[2]] += 1

arr_N就是存放了每个格点台风频数的数组了，第一维是lon，第二维是lat

最后画图看看

# 创建Figure
fig = plt.figure(figsize=(9, 6))

ax1 = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=0))

# 设置绘图区域的经纬度，最好和数据的范围一样
latlon=[100,180,0,60]

# 设置ax1的绘图范围

ax1.set_extent(latlon)
ax1.set_xlim(latlon[:2])
ax1.set_ylim(latlon[2:])

# 为ax1添加地理经纬度标签及刻度

ax1.set_xticks(np.arange(latlon[0],latlon[1]+1,10), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(latlon[2],latlon[3]+1,10), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())

# 为ax1添加填色

cf1 = ax1.contourf(lon,lat,arr_N.T,cmap=cm.WhiteBlueGreenYellowRed,levels=np.arange(0,50,0.1),extend='max')

#为ax1添加海岸线

ax1.coastlines()

# 色标
fig.colorbar(cf1,ax=ax1,ticks=np.arange(0,51,5),fraction=0.035,label="TC Count")

# 标注
ax1.text(102,56,f"TC Count : {len(TC_list)}",fontsize=13,color='k')
ax1.text(155,1,"Data From : IBTrACS-JTWC",fontsize=10,color='k')

# 标题
title='WNP_TC_Count_1979-2020'
plt.title(title,fontsize=18)

 大功告成~~

欢迎交流

祝大家科研顺利~