WRF后处理：模拟结果插值到站点(python版)
之前一篇文章介绍了如何使用NCL将WRF模拟结果插值到站点，包括特定的高度层和气压层。尽管NCL仍为WRF模式后处理最佳语言之一，但是随着python的使用逐渐广泛，我们需要逐渐将代码转向python版本。本文介绍如何使用python，实现WRF模拟结果插值到站点，包括不同的气压层和高度层。
实现WRF模拟结果插值到站点主要需要两个功能：一是寻找距离站点最近的网格点，通过编写一个函数实现。二是垂直插值功能，用于插值到特定高度或气压层。垂直层的插值，这里用到了wrf-python库，此库也是NCAR团队开发，用于WRF模式输出的众多诊断变量和部分插值例程，提供了 30 多个诊断计算、多个插值例程和实用程序，以帮助通过 cartopy、底图或 PyNGL 进行绘图，实现许多类似于NCL提供的功能。更多内容可访问：https://wrf-python.readthedocs.io/en/latest/index.html
一 .将u10, v10风插值到站点
给定一个站点经纬度，寻找WRF网格中最近点的索引，NCL中使用了
wrf_user_ll_to_xy函数。为了实现此功能，这里编写了如下函数：
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)上述函数，给定wrf网格的二维经度和纬度，以及指定站点的经纬度，最后返回(y， x)索引。
wrf-python也有内置的函数实现上述功能：
import wrf
wrfout_dir= "wrfout_d01"
lonSta = 96.0# 站点经纬度
latSta = 45.0
#返回的iSta, jSta为X，Y方向的索引
iSta, jSta = wrf.ll_to_xy(wrfout, latSta, lonSta, timeidx=0, meta=False)
具体使用如下：
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)
#寻找站点索引，法1
indexSta = nearest_position(latSta, lonSta, lat2D, lon2D)
jSta = indexSta
iSta = indexSta
#寻找站点索引，法2
iSta, jSta = wrf.ll_to_xy(wrfout, latSta, lonSta, timeidx=0, meta=False)
u10 = wrfout.variables['U10'][:]
v10 = wrfout.variables['V10'][:]
uSta = u10
vSta = v10
二 .提取站点距地面100m高度的U, V风
思路简述：(1)读入数据(三维uv风，位势高度和地形高度)，(2)将位势高度减去地形高度，得到模式各层数据距离地面的高度(AGL，Above Ground Level )，(3)以AGL高度作为参考系，将三维uv风插值到目标高度层(100m)，(4)最后根据前述找到的离站点最近的索引，读取站点特定高度的风速。
#读取3维u，v， 位势高度和地形高度
ua = to_np(getvar(wrfout, "ua" , timeidx=ALL_TIMES, meta=False, units="m s-1"))
va = to_np(getvar(wrfout, "va" , timeidx=ALL_TIMES, meta=False, units="m s-1"))
height = to_np(getvar(wrfout, "z", timeidx=ALL_TIMES, meta=False, units="m" ))
ter = to_np(getvar(wrfout, "ter", timeidx=0 , meta=False, units="m" ))
height_temp = height - ter
u_plane = interplevel(ua, height_temp, target_height, meta=False)
v_plane = interplevel(va, height_temp, target_height, meta=False)
uSta = u_plane
vSta = v_plane
三 .提取站点500hPa高度的U, V风
思路简述：(1)读入数据(三维uv风，气压)，(2)以气压作为参考系，将三维uv风插值到目标气压层(500hPa)，(3)最后根据前述找到的离站点最近的索引，读取站点特定气压层的风速。
ua = to_np(getvar(wrfout, "ua" , timeidx=ALL_TIMES, meta=False, units="m s-1"))
va = to_np(getvar(wrfout, "va" , timeidx=ALL_TIMES, meta=False, units="m s-1"))
pres = to_np(getvar(wrfout, "p", timeidx=ALL_TIMES, meta=False, units="hPa"))
u_plane = interplevel(ua, pres, target_pres, meta=False)
v_plane = interplevel(va, pres, target_pres, meta=False)
uSta = u_plane
vSta = v_plane
文章来源于微信公众号：气海同途
页:
[1]