之前一篇文章介绍了如何使用NCL将WRF模拟结果插值到站点，包括特定的高度层和气压层。尽管NCL仍为WRF模式后处理最佳语言之一，但是随着python的使用逐渐广泛，我们需要逐渐将代码转向python版本。本文介绍如何使用python，实现WRF模拟结果插值到站点，包括不同的气压层和高度层。

实现WRF模拟结果插值到站点主要需要两个功能：一是寻找距离站点最近的网格点，通过编写一个函数实现。二是垂直插值功能，用于插值到特定高度或气压层。垂直层的插值，这里用到了wrf-python库，此库也是NCAR团队开发，用于WRF模式输出的众多诊断变量和部分插值例程，提供了 30 多个诊断计算、多个插值例程和实用程序，以帮助通过 cartopy、底图或 PyNGL 进行绘图，实现许多类似于NCL提供的功能。更多内容可访问：https://wrfpython.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)索引。

具体使用如下：

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)
indexSta = nearest_position(latSta, lonSta, lat2D, lon2D)
jSta = indexSta[0]
iSta = indexSta[1]
u10 = wrfout.variables['U10'][:]
v10 = wrfout.variables['V10'][:]
uSta = u10[it, jSta, iSta]
vSta = v10[it, jSta, iSta]

二 提取站点距地面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[it,:,:,:] - ter
u_plane = interplevel(ua[it,:,:], height_temp, target_height, meta=False)
v_plane = interplevel(va[it,:,:], height_temp, target_height, meta=False)
uSta = u_plane[jSta, iSta]
vSta = v_plane[jSta, iSta]

三 提取站点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[it,:,:], pres[it,:,:,:], target_pres, meta=False)
v_plane = interplevel(va[it,:,:], pres[it,:,:,:], target_pres, meta=False)
uSta = u_plane[jSta, iSta]
vSta = v_plane[jSta, iSta]