WRF后处理：模拟结果插值到站点(python版)_寻找离站点最近的格点_气海同途的博客

link之家

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

相关文章推荐

粗眉毛的柑橘 · .MP4、.ARF 與 .WRF ...· 10 小时前 ·

鼻子大的汤圆 · Colab带进度条的拷贝和解压 - 知乎· 1 年前 ·

冷静的咖啡 · Spring-Security基于源码扩展- ...· 2 年前 ·

鼻子大的松鼠 · 创建并配置事实逻辑表_智能数据建设与治理 ...· 2 年前 ·

之前一篇文章介绍了如何使用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 degreeslon2D       = to_np(getvar(wrfout, "lon"  ))  # units: decimal degreestimes       = 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,:,:,:] - teru_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] 
以上代码的完整版本可以识别下方二维码，关注气海同途公众号，后台回复"python站点"关键字获取。 
                                    WRF模式、空气质量预报、生态、水资源、风资源、碳中和、水文、气候变化、WRF-CHEM、资料同化（WRF-DA）、水文过程（WRF-HYDRO）、城市化（URBAN）、环境
WRF气象预报模型、python、NCL
                                    wrf.getvar(wrfnc, varname, timeidx=0, method=u’cat’, squeeze=True, cache=None, meta=True, **kargs)
Returns basic diagnostics from the WRF ARW model output.
Below is a table of available diagnostics.
Variable Name	Description	                               
                                    WRF模式运行结束后会输出wrfout_d01*结果文件，利用NCL直接进行读取。addfile()函数打开文件看一眼变量和维度信息，包含wrf设定的参数信息等。利用wrf_user_getvar()函数读取文件中的气象变量资料，并简单绘制出图。......
                                    Arcpy简介之离散点插值离散点插值一直在使用Python的读者应该会发现一个问题：大部分的离散点插值为网格数据的Python第三方库都相当不好用(其中包括了Scipy)，经常会由于数据量过大、内存不足等原因造成程序崩溃。那么我们是否能够调用Arcgis工具箱里各种强大的插值函数来完成该任务呢？答案是肯定的。Arcpy为用户提供的主要模块之一名为arcpy.sa，是spatialanalyst的简...
in=addfile("wrfout_d02","r")target_pres=500;目标气压层500hPalon_sta=116;站点经纬度lat_sta=35u...
                                    1.格点数据有时需要插值到另一个范围和分辨率的网格上，可以用 linint2 函数来插值。
http://bbs.06climate.com/forum.php?mod=viewthread&tid=91352
2.格点插值到站点
格点数据插值到站点主要有两种方法：双线性插值和最近距离，算法都很简单，MeteoInfoLab中插值到站点有几种方法：（a）利用DimDataFile的tost...
                                    在我使用的show final solution的答案之后编辑了这个问题我有来自不同来源的非结构化2D数据集,例如： 这些数据集是3 numpy.ndarray(X,Y坐标和Z值).我的最终目标是在网格上插入这些数据以转换为图像/矩阵.所以,我需要找到插入这些数据的“最佳网格”.而且,为此,我需要在该网格的像素之间找到最佳的X和Y步长.根据点之间的欧氏距离确定步骤：使用每个点与其最近邻居之间的欧几...
                                    本文主要介绍python对wrfout结果文件的初步后处理操作，以及基础绘图。
wrfout后处理包括：【读取wrfout文件、读取wrfout文件中变量metadata及数据、对高度场进行500hPa插值、输出nc文件】
基础绘图操作包括：【设置投影和范围、绘制等值线contour和等值线标值、副高区域填色contourf】
仅展示初步评估模拟的效果，若精美绘制需要进一步的设置、细化。
                                    You are given a sequence of numbers a1, a2, ..., an, and a number m.
Check if it is possible to choose a non-empty subsequence aij such that the sum of numbers in this subsequence is divisible by m