用CMA热带气旋最佳路径数据集计算南海台风PDI指数

PDI（power dissipation index），中文直译总功率耗散指数，是表征热带气旋破坏潜力的指标，由麻省理工学院的Kerry Emanuel教授在2005年提出，综合考虑了热带气旋的数量、强度和持续时间，是一个比较综合全面的指数。（Increasing destructiveness of tropical cyclones over the past 30 years | Nature）

PDI的计算公式简单直观，物理意义明确，但网上搜不到计算台风PDI指数的Python代码，那就当作笔记在这里发一个吧。献丑了。

计算分为两个部分：第一部分提取出历史上的南海台风编号，第二部分计算PDI指数。

#需要用到的库 from typing import List from typing import Union from typing import Tuple import pandas as pd import numpy as np from pathlib import Path from scipy.integrate import trapz,simps import csv import matplotlib.pyplot as plt import matplotlib.path as mpath import datetime

第一部分：提取历史南海台风编号

先定义好最佳路径数据集的文件名放到列表里备用

#先定义需要计算的年份 begYear = 1949 endYear = 2020 totalYear = endYear-begYear+1# 获取指定年份的最佳路径数据的文件名 fileName = [] for i in range(begYear,endYear+1):iFile = 'CH'+str(i)+'BST.txt'fileName.append(iFile)

定义南海的边界范围

#南海边界经纬度（或者用更精细的地图文件更好） limlonlat = (105,120,5,20) #圈定南海边界 boundary = [(limlonlat[0], limlonlat[2]), (limlonlat[1], limlonlat[2]), (limlonlat[1], limlonlat[3]), (limlonlat[0], limlonlat[3]),(0,0)] #定义好一个南海经纬度范围的path类 path = mpath.Path(boundary, closed=True)

 读取最佳路径数据集，保存影响南海的台风编号和影响时段的台风信息

SCSTC = [] # 保存影响的台风编号 SCSInfluLine = [] #保存影响时段的台风信息 #逐年读取最佳路径数据 for i in range(0,totalYear):filePath = fileName[i]fileRead = open(filePath) while True:line = fileRead.readline().split()if not line:break#头记录if line[0] == "66666": numberTy = line[4] # typhoon numbercontinue#最佳路径数据记录newLine = ['numberCN', 'yyyymmddhh', 'latRec', 'lonRec', 'presRec','wndRec','gradeRec']#为了数据统一，只提取00，06，12，18时刻的信息if line[0][8:10] in ["00","06","12","18"]: numberTy = numberTyyyymmddhh = line[0]latRec = float(line[2]) * 0.1 #unit 1.0 degreelonRec = float(line[3]) * 0.1presRec = line[4]wndRec = line[5]gradeRec = line[1]newLine[0] = numberTynewLine[1] = yyymmddhhnewLine[2] = str(latRec)newLine[3] = str(lonRec)newLine[4] = presRecnewLine[5] = wndRecnewLine[6] = gradeRec#跳过未编号台风if numberTy == "0000":continueif path.contains_point((lonRec,latRec),radius=0.01):# 计算进入边界范围内的台风SCSInfluLine.append(newLine)if numberTy not in SCSTC:SCSTC.append(numberTy)else:# 影响范围外的台风，不计入continuefileRead.close

至此，已提取出历史南海台风的编号，放在名为SCSTC的列表中

第二部分：计算PDI指数

为了方便逐个台风地计算，定义一个reader函数（参考了大神的代码），用于读取指定文件和指定台风的最佳路径数据

#读取CMA热带气旋最佳路径数据集 def reader(typhoon_txt: os.PathLike, code: Union[str, int]) -> Tuple[List[str], pd.DataFrame]:typhoon_txt = Path(typhoon_txt)if isinstance(code, int):code = "{:04}".format(code)with open(typhoon_txt, "r") as txt_handle:while True:header = txt_handle.readline().split()if not header:raise ValueError(f"没有在文件里找到编号为{code}的台风")if header[4].strip() == code:break[txt_handle.readline() for _ in range(int(header[2]))]data_path = pd.read_table(txt_handle,sep=r"\s+",header=None,names=["TIME", "I", "LAT", "LONG", "PRES", "WND", "OWD"],nrows=int(header[2]),dtype={"I": np.int8 ,"LAT": np.float32,"LONG": np.float32,"PRES": np.float32,"WND": np.float32,"OWD": np.float32,},parse_dates=True,date_parser=lambda x: pd.to_datetime(x, format=f"%Y%m%d%H"),index_col="TIME",)data_path["LAT"] = data_path["LAT"] / 10data_path["LONG"] = data_path["LONG"] / 10return header, data_path

 万事俱备，正式计算PDI

PDI_allyear=[] #保存逐年的PDI PDI_allTC=[] #保存每个台风的PDI# 获取指定年份的最佳路径数据 for year in range(begYear,endYear+1):File = 'CH'+str(year)+'BST.txt'pdi_eachyear=0# 遍历影响南海的台风列表for TCnum in SCSTC:if str(year)[2:] == TCnum[:2]:try:head, dat = reader(File,TCnum)# #为了数据统一，只提取00，06，12，18时刻的信息dat['every6']=[str(i)[11:13] in ['00','06','12','18'] for i in dat.index]dat=dat[dat['every6']==True]# 积分计算每个台风的pdipdi = trapz(dat.WND**3,dx=6*60*60)PDI_allTC.append(pdi)pdi_eachyear += pdi except:passelse:continuePDI_allyear.append(pdi_eachyear)

大功告成！

PDI_allTC表示的是每个台风的PDI指数

PDI_allyear表示每年的PDI指数

如果有更简练的代码也欢迎交流~

画图的代码略

总结

以上是生活随笔为你收集整理的用CMA热带气旋最佳路径数据集计算南海台风PDI指数的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。