Python提取WRF模拟的台风路径

WRF模式广泛应用于台风的模拟和预报，但是其并不能输出台风路径信息（除非编译的时候nesting选择vortex following选项，运行时会输出ATCF数据，不过这个不常用）。

因此使用WRF模拟台风，后处理需要解决如何从wrfout数据提取出台风路径的问题。NCL官网提供的脚本中，是以每个时刻海平面气压(SLP)最小处作为台风的中心位置，简单易行。但是在一些特定情况存在问题，比如双台风的时候。

为了避免上述的问题，提取台风在T时的位置时必须要参考T-1时刻的位置。即以T-1时刻台风的位置为圆心，给定半径radius确定一个圆，在此范围内寻找T时刻的台风位置。以下为具体的操作步骤：

（1）第1时刻，找到min(SLP)的的位置索引，然后提取台风中心经纬度(lonTC, latTC)。或者由观测或者最佳路径数据集给定模拟初始时刻的台风中心位置。

（2）第it个时刻，根据it-1时刻的台风中心位置(lonTC, latTC)计算it-1时刻的台风中心位置索引(ic, jc)。北纬30°以南，台风一般移动速度spd不超过0.5deg/hour，北纬30°以北一般不超过2.0deg/hour，wrfout输出时间间隔为history_interval(hours)。根据移速和输出间隔，计算台风最大移动半径radius=spd*history_interval。最大移动半径除以网格分辨率，得到移动半径所对应的索引半径indexRadius。根据索引半径确定台风所在的范围，此范围内的SLP最小处即是it时刻的台风中心位置。

T时刻的台风中心，应该在T-1时刻台风中心的附近，根据T时刻的SLP场在这个范围内的最小值位置，确定T时刻的台风中心

（3）遍历所有时刻，得到完整的台风路径，并且每一步可以保存最低气压Pmin，最大风速Vmax等信息。

以1713号台风天鸽(Hato)为例，进行模拟。WRF模式分辨率为15km，模拟起止时间：2017-08-21_00 至 2017-08-24_00 UTC。

从此次模拟结果看，基本模拟出了台风"天鸽"的大致走向，只是模拟的台风移速偏慢。

模拟的台风路径与观测对比

该算法需要使用海平面气压(SLP)，但是WRF的输出数据中是没有海平面气压的，因此需要根据已有的气压和温度场等信息诊断得到SLP。使用wrf-python库的getvar函数可以直接实现这个功能，直接得到SLP。

PS1: 实际应用发现一个问题，WRF初始时刻诊断得到的SLP，台风中心特别不明显（甚至不存在）。如果模拟区域存在高山，比如台湾的山脉，则往往山脉的SLP最低，因此根据SLP最小确定台风的位置则容易出错，即导致台风位置错误的出现在高山上。因此初始时刻最好给定台风的大致位置，并且最大移动半径设置较小(<0.5°)。如果无法给出初始时刻的台风位置，可以跳过初始时刻，从第二个时刻开始。

初始时刻SLP诊断场，无明显的台风中心，但实际上此时台风已经比较成熟。

初始时刻基于SLP最小确定台风中心位置，SLP场的结果较差导致台风位置出现在山脉。

PS2: 实际使用，如果双台风距离太近，台风之间的距离小于最大移动半径时，也可能会发生误判，此类情况可以缩小最大移动半径。

以下为全部代码：

import numpy as npimport netCDF4 as ncimport datetimefrom wrf import getvar, ALL_TIMES, to_npimport matplotlib.pyplot as pltimport cartopy.crs as ccrsfrom cartopy.io.shapereader import Readerimport cartopy.feature as cfeaturefrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterimport pandas as pd
def nearest_position( stn_lat, stn_lon, lat2d, lon2d):    """获取最临近格点坐标索引    stn_lat  : 站点纬度    stn_lon  : 站点经度    lat2d    : numpy.ndarray网格二维经度坐标    lon2d    : numpy.ndarray网格二维纬度坐标    Return: y_index, x_index    """    difflat = stn_lat - lat2d;    difflon = 


    
stn_lon - lon2d;    rad = np.multiply(difflat,difflat)+np.multiply(difflon , difflon)    aa=np.where(rad==np.min(rad))    ind=np.squeeze(np.array(aa))    return tuple(ind)
if __name__  == "__main__":    # setting    wrfout_dir  = "wrf_wnp15km/wrfout_d01_2017-08-21_00:00:00"    lonTC0 = 125.0  #台风初始时刻位置    latTC0 = 20.3
    wrfout      = nc.Dataset(wrfout_dir, mode="r")    lat2D       = to_np(getvar(wrfout, "lat"  ))  # units: decimal degrees    lon2D       = to_np(getvar(wrfout, "lon"  ))  # units: decimal degrees    times       = to_np(getvar(wrfout, "times", timeidx=ALL_TIMES))  #    nt = len(times)    ny, nx = np.shape(lat2D)
    date0 = datetime.datetime.strptime(str(times[0])[0:19], "%Y-%m-%dT%H:%M:%S")    date1 = datetime.datetime.strptime(str(times[1])[0:19], "%Y-%m-%dT%H:%M:%S")    history_interval = (date1 - date0).total_seconds()/3600  #hours
    lonMax = np.max(lon2D)    lonMin = np.min(lon2D)    latMax = np.max(lat2D)    latMin = np.min(lat2D)
    date = [] # 2020-07-20T00    lons = [] # degree    lats = [] # degree    pmin = [] # hPa    vmax = [] # m/s
    lonTC = lonTC0    latTC = latTC0
    for it in range(nt):        slp         = to_np(getvar(wrfout ,"slp"         ,units="hPa"   ,timeidx=it))        wspd_wdir10 = to_np(getvar(wrfout ,"wspd_wdir10" ,units="m s-1" ,timeidx=it))        wspd10      = wspd_wdir10[0]
        if it==0 :            # 如果给定了初始时刻的TC位置，则使用给定的TC位置            # 未给定初始时刻TC位置，则取全场slp最小的位置为TC中心            if latTC0 > -60.0 :                latTC = latTC0                lonTC = lonTC0                #lons.append(round(lonTC,2))                #lats.append(round(latTC,2))                #continue            else:                slpMin = np.min(slp[:,:])                indexMin = np.argwhere(slp[:,:] == slpMin)                jTC = indexMin[0][0]                iTC = indexMin[0][1]                lonTC = lon2D[jTC,iTC]                latTC = lat2D[jTC,iTC]
        ### 1 找到TC中心(lonTC, latTC)的索引(iTC, jTC)        indexTC = nearest_position(latTC, lonTC, lat2D, lon2D)        jTC = indexTC[0]        iTC = indexTC[1]        # 避免台风中心选在边界点        jTC = np.max((1,jTC))    # jTC [1,ny-2]        jTC = np.min((jTC,ny-2))


    
        iTC = np.max((1,iTC))    # iTC [1,nx-2]        iTC = np.min((iTC,nx-2))
        ### 2 计算TC center附近的网格分辨率dAvg，        dLat = lat2D[jTC+1,iTC] - lat2D[jTC,iTC]        dLon = lon2D[jTC,iTC+1] - lon2D[jTC,iTC]        dAvg = (dLat + dLon)/2.0
        ### 3 根据移速计算台风中心最大可能半径，根据这个        if latTC < 30.0 : # 纬度30°以南，台风移速每小时0.5°           radius = 0.5*history_interval  # 0.5 degree/hour        else:             # 纬度30°以北，台风移速每小时1.0°           radius = 1.0*history_interval  # 1.0 degree/hour           radius = 0.5*history_interval  # 0.5 degree/hour        radius = 1.0*history_interval  # 1.0 degree/hour        if it==0 :           radius = 0.5        indexRadius = int(radius/dAvg) + 1
        ### 找到最大可能半径内，slp最小值及其位置索引        iStart = iTC - indexRadius        iEnd   = iTC + indexRadius        jStart = jTC - indexRadius        jEnd   = jTC + indexRadius        jStart = np.max((1,jStart))        jEnd   = np.min((jEnd,ny-2))        iStart = np.max((1,iStart))        iEnd   = np.min((iEnd,nx-2))
        slpMin = np.min(slp[jStart:jEnd,iStart:iEnd])        w10Max = np.max(wspd10[jStart:jEnd,iStart:iEnd])        indexMin = np.argwhere(slp[jStart:jEnd,iStart:iEnd] == slpMin)        jTC = indexMin[0][0] + jStart        iTC = indexMin[0][1] + iStart        lonTC = lon2D[jTC,iTC]        latTC = lat2D[jTC,iTC]        print("date:", str(times[it])[0:19],"TC center:",round(lonTC,2), round(latTC,2)," p min:",round(slpMin,2), " vmax:",round(w10Max,2))        date.append(str(times[it])[0:19])        lons.append(round(lonTC,2))        lats.append(round(latTC,2))        pmin.append(round(slpMin,2))        vmax.append(round(w10Max,2))
    #read CMA-STI best track    f_Con = 'cma_bst_hato.txt'    col_names =['date','grade','lat', 'lon', 'pres', 'vmax']    widths = [10,2,4,5,4,3]    df = pd.read_fwf(f_Con,usecols=[0,1,2,3,4,5],widths=widths,names=col_names)    latObs  = df['lat'].values[:]    lonObs  = df['lon'].values[:]    latObs  = np.array(latObs)/10    lonObs  = np.array(lonObs)/10
    ### plot track    fig = plt.figure()    ax = plt.axes(projection=ccrs.PlateCarree())    ax.set_extent([lonMin, lonMax, latMin, latMax],crs=ccrs.PlateCarree())    xticks = np.linspace(lonMin, lonMax, 6, endpoint=True)    yticks = np.linspace(latMin, latMax, 6, endpoint=True)    ax.set_xticks(xticks, crs=ccrs.PlateCarree())    ax.set_yticks(yticks, crs=ccrs.PlateCarree())    lon_formatter = LongitudeFormatter(zero_direction_label=False)    


    
lat_formatter = LatitudeFormatter()    ax.xaxis.set_major_formatter(lon_formatter)    ax.yaxis.set_major_formatter(lat_formatter)    ax.gridlines()    ax.coastlines()    for it in range(1,len(lons)):        if it == 1:           plt.plot((lons[it-1],lons[it]), (lats[it-1],lats[it]),color='red',linewidth=1.5,transform=ccrs.PlateCarree(), label="wrf")        else:           plt.plot((lons[it-1],lons[it]), (lats[it-1],lats[it]),color='red',linewidth=1.5,transform=ccrs.PlateCarree())    for it in range(1,len(lonObs)):        if it == 1:            plt.plot((lonObs[it-1],lonObs[it]), (latObs[it-1],latObs[it]),color='black',linewidth=2,transform=ccrs.PlateCarree(), label="obs")        else:            plt.plot((lonObs[it-1],lonObs[it]), (latObs[it-1],latObs[it]),color='black',linewidth=2,transform=ccrs.PlateCarree())    plt.legend()    plt.show()

关注下方公众号，获取更多大气海洋领域的内容。

往期推荐

构建适合大气与海洋应用的Anaconda环境

WRF后处理：降雨量的说明以及降雨的绘制

python爬取中央气象台台风网当前台风实况和预报数据