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 np
import netCDF4 as nc
import datetime
from wrf import getvar, ALL_TIMES, to_np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import 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), "%Y-%m-%dT%H:%M:%S")
date1 = datetime.datetime.strptime(str(times), "%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
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
iTC = indexMin
lonTC = lon2D
latTC = lat2D
### 1 找到TC中心(lonTC, latTC)的索引(iTC, jTC)
indexTC = nearest_position(latTC, lonTC, lat2D, lon2D)
jTC = indexTC
iTC = indexTC
# 避免台风中心选在边界点
jTC = np.max((1,jTC)) # jTC
jTC = np.min((jTC,ny-2))
iTC = np.max((1,iTC)) # iTC
iTC = np.min((iTC,nx-2))
### 2 计算TC center附近的网格分辨率dAvg，
dLat = lat2D - lat2D
dLon = lon2D - lon2D
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)
w10Max = np.max(wspd10)
indexMin = np.argwhere(slp == slpMin)
jTC = indexMin + jStart
iTC = indexMin + iStart
lonTC = lon2D
latTC = lat2D
print("date:", str(times),"TC center:",round(lonTC,2), round(latTC,2)," p min:",round(slpMin,2), " vmax:",round(w10Max,2))
date.append(str(times))
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 =
df = pd.read_fwf(f_Con,usecols=,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(,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,lons), (lats,lats),color='red',linewidth=1.5,transform=ccrs.PlateCarree(), label="wrf")
else:
plt.plot((lons,lons), (lats,lats),color='red',linewidth=1.5,transform=ccrs.PlateCarree())
for it in range(1,len(lonObs)):
if it == 1:
plt.plot((lonObs,lonObs), (latObs,latObs),color='black',linewidth=2,transform=ccrs.PlateCarree(), label="obs")
else:
plt.plot((lonObs,lonObs), (latObs,latObs),color='black',linewidth=2,transform=ccrs.PlateCarree())
plt.legend()
plt.show()
文章来源于微信公众号：气海同途
New Google Indexing Info
To the person talking about inspected url, check google indexing status, get your website listed on google, crawl your site, view page as googlebot, testing if live url can be indexed, get website to show up on google, check website google indexing, bulk indexing tool, search engine visibility wordpress plugin,I highly suggest this awesome google indexing info or crawlable links, check which pages are indexed by google, crawled not currently indexed, site index check, free backlink indexing service, search engine basics, add my website to search engines, bulk indexing tool, backlink importance, google googlebot, as well as this top rated google indexing tips on top of my new website is not showing on google, google is not showing my website, crawl rate search console, website indexing tool, seo backlinks example, google search console issues, website not indexing, check if page is indexed by google, index now, index my website on google, on top of this discover more about google indexing details which is also great. Also, have a look at this read full article for google indexing advice as well as tool index google, indexnow seo, make google index my website, search backlinks on google, submit new url to google, google submit, my wordpress site is not indexed by google, add your site to google search, google search console tool, crawl as googlebot, not to mention this redirected here on google indexing site with noindex tag check, check if website can be crawled, deindex from google, search engine of google, google search how does it work, google console check website, index checker online, google console errors, add website to search, google website search tool,for good measure. Check more @ Recommended TAJIR4D Guide 1a7dc8c
页:
[1]