用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 number
            continue
        
        #最佳路径数据记录
        newLine = ['numberCN', 'yyyymmddhh', 'latRec', 'lonRec', 'presRec','wndRec','gradeRec']
        #为了数据统一，只提取00，06，12，18时刻的信息
        if line[0][8:10] in ["00","06","12","18"]: 
            numberTy  = numberTy
            yyymmddhh = line[0]
            latRec    = float(line[2]) * 0.1 #unit 1.0 degree
            lonRec    = float(line[3]) * 0.1
            presRec   = line[4]
            wndRec    = line[5]
            gradeRec   = line[1]
            newLine[0] = numberTy
            newLine[1] = yyymmddhh
            newLine[2] = str(latRec)
            newLine[3] = str(lonRec)
            newLine[4] = presRec
            newLine[5] = wndRec
            newLine[6] = gradeRec
            
            #跳过未编号台风
            if numberTy == "0000":
                continue
            
            if path.contains_point((lonRec,latRec),radius=0.01):
            # 计算进入边界范围内的台风
                SCSInfluLine.append(newLine)
                if numberTy not in SCSTC:
                    SCSTC.append(numberTy)
            
            else:
            # 影响范围外的台风，不计入
                continue
                
    fileRead.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"] / 10
        data_path["LONG"] = data_path["LONG"] / 10
        return 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]
                # 积分计算每个台风的pdi
                pdi = trapz(dat.WND**3,dx=6*60*60)
                PDI_allTC.append(pdi)
                pdi_eachyear += pdi 
            except:
                pass
        else:
            continue
    PDI_allyear.append(pdi_eachyear)

大功告成！

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

PDI_allyear表示每年的PDI指数

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

画图的代码略