WRF后处理 |任意高度绘图

WRF模拟时，在namelist.input文件中，有一个参数叫做e_vert，表示模拟时垂直的层数。但因为模拟时下垫面的海拔高度不一样，有的为平原，海拔高度较低。有的为高原，海拔高度较高。于是，WRF模式在垂直上主要采用的eta坐标。
eta坐标是通过以下公式定义的：
$eta=(p-p_{\text{top}})/(p_{\text{surf}}-p_{\text{top}})$
其中：

• $p$是在ETA坐标上的压强（Pa）。
• $p_{\text{top}}$是大气柱的顶部压强（Pa）。
• $p_{\text{surf}}$是地面处的压强（Pa）。
• $eta$是非均匀网格上的坐标，它的值在大气柱内从1到0变化。

但是，我们在画图时可能更想画某一高度的要素值的空间变化，这是就需要进行一下数据处理，得到其相对应的每一层的位势高度，主要有两种方法：
第一种：高度的公式转换
$gmp=(PH+PHB)/9.81-HGT$
相关代码：

P=ph+phb
P = destagger(P,0,meta=True)
gmp=P/9.81-hgt

第二种：利用wrf-python包
利用wrf-python的诊断变量，可以轻松获得位势高度。

z=getvar(ncfile, "height_agl")

获得位势高度后，进行插值，即可获得目标高度的要素值，完整参考代码如下：

from netCDF4 import Dataset
import numpy as np
from wrf import getvar,interpz3d,destagger


# 打开NetCDF文件
file_path = r'H:\wrfout_d01_2023-04-08_00_00_00'
ncfile = Dataset(file_path, 'r')

ph=getvar(ncfile, "PH")
phb=getvar(ncfile, "PHB")
hgt=getvar(ncfile, "HGT")
p = getvar(ncfile, "p")

wspd = getvar(ncfile, "wspd_wdir10", units="m/s")[0,:]


# 按高度层插值
P=ph+phb
P = destagger(P,0,meta=True)
gmp=P/9.81-hgt

z_list=[100,200,300] #目标高度unit:m
      
p_z = interpz3d(p,gmp,np.array(z_list))