WRF模拟时,在namelist.input文件中,有一个参数叫做e_vert,表示模拟时垂直的层数。但因为模拟时下垫面的海拔高度不一样,有的为平原,海拔高度较低。有的为高原,海拔高度较高。于是,WRF模式在垂直上主要采用的eta坐标。
eta坐标是通过以下公式定义的:
e
t
a
=
(
p
−
p
top
)
/
(
p
surf
−
p
top
)
eta = (p-p_{\text{top}})/(p_{\text{surf}} - p_{\text{top}})
eta=(p−ptop)/(psurf−ptop)
其中:
- p p p是在ETA坐标上的压强(Pa)。
- p top p_{\text{top}} ptop是大气柱的顶部压强(Pa)。
- p surf p_{\text{surf}} psurf是地面处的压强(Pa)。
- e t a eta eta是非均匀网格上的坐标,它的值在大气柱内从1到0变化。
但是,我们在画图时可能更想画某一高度的要素值的空间变化,这是就需要进行一下数据处理,得到其相对应的每一层的位势高度,主要有两种方法:
第一种:高度的公式转换
g
m
p
=
(
P
H
+
P
H
B
)
/
9.81
−
H
G
T
gmp=(PH+PHB)/9.81-HGT
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))