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指数的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。