WRF byp.com.cn

WRF
Back

WRF

WRF模式的结果一般为三个输出文件（如下图），对应着三个网格的结果。

d01为最外一层网格（45公里），d02为第二层网格（9公里），d03为第三层网格（3公里）。预报时效为72小时逐时

如果使用python做后处理，wrf-python还是很方便的。例如使用wrf.getvar就可以方便读出诊断变量和原始变量（源码），但是要注意使用参数timeidx（默认为0，全部时次为timeidx=wrf.ALL_TIMES）

wrf-python的安装有点麻烦，在windows环境下一直提示“出错，要安装“numpy+mkl”。因此，需要pythonlib手动下载与你python版本匹配的安装包（[文件：numpy-1.22.4+mkl-cp310-cp310-win_amd64.whl]）安装（245MB，下载速度很慢）。不过在linux环境下，最简单的方法就是安装和使用conda了！（见conda）

geogrid生成的geo_em.d01.nc文件，可以使用程序可以读取变量生成网页文件，使用程序画出以下地形图：

WRF的三个结果文件，可以用netCDF4读取。每个结果文件都包含有191个变量（variables）。该程序可以读取这些变量生成网页文件。

变量的读取有几种方法：

1、直接读取

例如读取风的变量，其实有5个变量：U、V、W、U10、V10，单位都是m s-1。

使用print(data_nc.variables['U'])就可以读取变量的信息。

使用print(data_nc.variables['U'].shape)可以看到该变量（d01网格）的维数是(73, 36, 69, 70)，73个时次、36层、69×70个网格。

使用print(data_nc.variables['U'][1][1][:].shape)可以查看第1个时次、第1层所有变量的维数是(69, 70)。

2、使用wrf.getvar

方法1：print(data_nc.variables['U'][1][1][1][1])
方法2：
u = wrf.getvar(data_nc, 'U', timeidx=1)
print(u.values[1][1][1])
以上两个方法的读数是一致的。

wrf.getvar读取的数据格式是xarray.DataArray，一种使用标签的多维数组。数组主要包括四种关键属性：

（1）values：数组值的numpy.ndarray；（2）dims：每个坐标轴的维度名称；（3）coords：包含数组坐标的类似字典的容器；（4）attrs：存放元数据属性的字典；

各种变量的读取

1、风

WRF的输出变量中有关风的变量就包括：U、V、W、U10、V10等5种变量。而wrf-python本身也提供的诊断变量（就是直接用wrf.getvar读取出来的数值）就包括：ua、va、wa、uvmet10、uvmet、wspd_wdir、wspd_wdir10、uvmet_wspd_wdir、uvmet10_wspd_wdir等9种诊断变量，并且直接提供了m s-1、km h-1、mi h-1、kt、ft s-1等5种单位（无需另外转换）。

以诊断变量（wrf-python计算处理的）uvmet10为例。

uvmet10按照官方的说明是：10m U and V Components of Wind Rotated to Earth Coordinates，即经过地球坐标投影转换的10米高度风的U和V分量值。

我们可以直接使用一下命令读取和打印结果（结果见截图）：

u = wrf.getvar(data_nc, 'uvmet10', timeidx=idx)
print(u)

使用print(u.attrs['projection'])就可以打印出“地球投影”：

需要注意的是一个数组是包含了所有预报时次的，时次需要使用timeidx参数来设置。例如下面的例子，第一个是选择了所有时间（timeidx=wrf.ALL_TIMES）；第二个是选择了第一个时次（timeidx=0）

u = wrf.getvar(data_nc, 'uvmet10', timeidx=wrf.ALL_TIMES)
Time: 73, south_north: 69, west_east: 69)>

u = wrf.getvar(data_nc, 'uvmet10', timeidx=0)

wrf-python诊断变量里面有关10米风的就有三个变量：uvmet10、wspd-wdir10、uvmet10_wsps_wdir，官方说明如下：

uvmet10:
10 m U and V Components of Wind Rotated to Earth Coordinates
wspd_wdir10:
10m Wind Speed and Direction (wind_from_direction) in Grid Coordinates
uvmet10_wspd_wdir:
10m Wind Speed and Direction (wind_from_direction) Rotated to Earth Coordinates

以上区别主要是坐标系统（Coordinates）。但是如果我们输出这三个变量的结果来看，其实坐标系是一样的！？（见下图|源码）。

其实读取WRF的默认变量U10和V10的投影也是和上述诊断变量一样的！

uvmet10的u和v可以经过计算得到风速值，与wspd_wdir10和uvmet10_wspd_wdir的风速值是一致的（见程序运行结果）。

总结10米风的变量：

U10和V10（注意是大写）是风的两个方向的分量，是WRF输出的2个独立的“自有变量”；而uvmet10（注意是小写）也是风的两个方向的分量，但是u和v的值都在一个读取结果里面，并且是wrf-python的“诊断变量”。

wspd_wdir10和uvmet10_wspd_wdir都是风的风速和风向的“诊断变量”。虽然在官网介绍中wspd_wdir10是Grid Coordinates，而uvmet10_wspd_wdir是Rotated to Earth Coordinates，但是打印出它们的内容时（如下图），是没有实际的数据差异的！

读取风对应的经纬度范围

经纬度的获取可以用wrf_getvar直接读取诊断变量XLONG和XLAT，也可以从读取风的诊断变量的坐标信息中获取，见程序。

风场作图

cartopy可以很方便的做风场风羽图（barbs，见程序），命令里面的参数可以调节风羽的size（间隔spacing、height为风羽上面的“羽毛”的长度），包括可以直接对空值的处理（见下图）。cartopy的风羽其实也是使用了matplotlib的命令quiver.Barbs。

下图是10米风场的风羽和风速等值线图（源码|叠加Logo）

读取wrfout的variable：源码01|源码02
读取nc文件信息：源码
读取slp的DataArray关键属性：源码01(values、dims、coords、attrs)|
读取风uv10：源码01|生成js
模型网格：提取最临近网格(wrf.ll_to_xy)源码01|生成海量点

Grib数据的读取

WRF的运行需要下载NCEP的FNL数据，而这些数据的格式是GRIB格式。

GRIB格式的数据主要是欧洲中长期预报中心使用，以前一般用pygrib读取。不过一般安装会出错（见错误信息)

比较简单的做法是使用xarray来读取，但是需要有cfgrib引擎（安装时提醒|示例|参考）。

不过如果直接使用cfgrib引擎来读，也会报错。

加上：ds = xr.open_dataset(name_file, engine='cfgrib', backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface'}})（程序），依然报错。

根据报错提示：

filter_by_keys={'stepType': 'instant', 'typeOfLevel': 'surface'}
filter_by_keys={'stepType': 'avg', 'typeOfLevel': 'surface'}
filter_by_keys={'stepType': 'accum', 'typeOfLevel': 'surface'}
更改为：

ds = xr.open_dataset(name_file, engine='cfgrib', backend_kwargs={'filter_by_keys': {'stepType': 'instant', 'typeOfLevel': 'surface'}})即可。

如果选择气压高度（层），typeOfLevel更改为isobaricInPa，并且赋值level，例如500hPa：

ds = xr.open_dataset(name_file, engine='cfgrib', backend_kwargs={'filter_by_keys': {'stepType': 'instant', 'typeOfLevel': 'isobaricInhPa', 'level': 500}})

WRF的模拟结果中不同来源的降雨是分开的，wrfout文件中，降雨的变量如下：

RAINC: 积云深对流过程产生的累积降水量，也就是模式中的积云对流参数化方案导致的降雨( cu_physics)。对于高分辨率的模拟，比如dx<5km，通常会将积云对流参数化方案关闭，此时RAINC为0。
RANNC: 此类降雨来源于云微物理参数化方案(mp_physics)，如大尺度抬升过程产生的凝结等微物理过程降水，也就是非对流产生的降水。
RAINSH: 积云对流参数化方案主要是反映深对流的降水过程，但是一些积云对流参数化方案，能够支持浅对流导致的降水，此时总降水还需要加上RAINSH。WRF中支持浅对流的参数化方案(cu_physics)有以下几种：KF，SAS，G3，BMJ，Tiedtke。WRF中也有独立于深对流过程的浅对流方案，通过namelist中设置shcu_physics。一般情况下，浅对流产生的降水量较小。

此外，固态降水例如雪(SNOWC/SNOWNC)、霰(GRAUPELC/GRAUPELNC)等降水，它们是降水的不同相态，已经都包含在RAINC/RAINNC中，不需要额外添加。

总结为一句话：总降水=RAINNC + RAINC + RAINSH，其中RAINC和RAINSH根据物理参数化方案的设置可能为0。

需要注意的是，降雨是个累积量，而不是瞬时值，在wrfout中如果需要得到模拟起止日期内的降雨，需要将最后一个时刻的降雨量减去第一个时刻的降雨量。