wrf-python 详解之如何使用

近几年，python在气象领域的发展也越来越快，同时出现了很多用于处理气象数据的python包。比如和NCL中的 WRF_ARWUser库类似的 wrf-python模块。

wrf-python是用于WRF模式后处理的python模块，其中提供了很多有用的函数，下面就来详细说一下其用法：

基本用法

• 计算诊断变量

wrf.getvar 函数的主要作用是返回需要计算的诊断变量，因为WRF本身不会返回这些变量。比如：CAPE（对流有效位能），SRH（风暴螺旋度）等等。

下例计算海平面压力：

from __future__ import print_function
from netCDF4 import Dataset
from wrf import getvar
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# 获取海平面压力
slp = getvar(ncfile, "slp")

• 提取WRF netCDF 变量

除了输出诊断变量外，wrf.getvar函数也可以用来提取常规的WRF输出的netCDF 变量。

p = getvar(ncfile, "P")

• 关闭 xarray 和 metadata

有时候你只需要返回常规的 numpy 数组，而不关心元数据。通过以下两种方式可以禁用元数据。

1. 完全禁用 xarray
2. 设置 meta 参数为 False

# 方法a
disable_xarray()
p_no_meta = getvar(ncfile, "P")
print (type(p_no_meta))
enable_xarray()
# 方法b
p_no_meta = getvar(ncfile, "P", meta=False)
print (type(p_no_meta))

• 从DataArray中提取 numpy 数组

如果你需要将 xarray.DataArray 转换为 numpy.ndarray， wrf-python中的 wrf.to_np 函数可以帮助你完成这一操作。尽管 xarray.DataArray 对象已经包含了 xarray.DataArray.values 属性用以提取 numpy 数组，但是用于编译扩展时仍会存在问题。因为 xarray 会将缺失值填充为 NaN，当用于编译扩展时会出错。还有就是一些程序可能可以用于 numpy.ma.MaskedArray，但含有 NaN 的numpy数组可能并不能工作。

wrf.to_np 函数按照以下流程执行：

1. 如果没有缺省值或填充值，那么将直接调用 xarray.DataArray.values 属性返回值
2. 如果有缺省值或填充值，那么会用 xarray.DataArray.attrs 属性 _FillValue 值替代 NaN 并返回 numpy.ma.MaskedArry

# 获取3D对流有效位能(包含缺省值)
cape_3d = getvar(ncfile, "cape_3d")
# 因为包含缺省值，因此返回mask数组
cape_3d_ndarray = to_np(cape_3d)
print(type(cape_3d_ndarray))

返回：

<class 'numpy.ma.core.MaskedArray'>

文件序列

• 使用 cat 方法合并多个文件

cat 方法会将序列中所有文件沿着 'Time' 维进行合并，时间维度将作为返回数组的最左侧维度。为了在输出数组中包含所有文件中的所有时间，设置 timeidx 参数为 wrf.ALL_TIMES（或设置为 None）。如果 timeidx 是单个值，那么将假设时间索引取自所有文件所有时间的连接。

注意 ：执行 wrf.getvar 时并不会进行排序，也就是说在执行函数之前应在序列中按时间对文件进行排序。

from __future__ import print_function
from netCDF4 import Dataset
from wrf import getvar, ALL_TIMES
# 创建文件序列
wrflist = [Dataset("wrfout_d01_2016-10-07_00_00_00"),
           Dataset("wrfout_d01_2016-10-07_01_00_00"),
           Dataset("wrfout_d01_2016-10-07_02_00_00")]
# 提取所有时刻的P变量
p_cat = getvar(wrflist, "P", timeidx=ALL_TIMES, method="cat")
print(p_cat)

返回：

<xarray.DataArray u'P' (Time: 3, bottom_top: 50, south_north: 1059, west_east: 1799)>
Coordinates:
    XLONG        (south_north, west_east) float32 -122.72 -122.693 -122.666 ...
    XLAT         (south_north, west_east) float32 21.1381 21.1451 21.1521 ...
  * Time         (Time) datetime64[ns] 2016-10-07 2016-10-07 2016-10-07
    datetime     (Time) datetime64[ns] 2016-10-07T00:00:00 ...

• 使用 join 方法组合多个文件

使用join方法合并一系列文件时，会将文件/序列索引作为新数组的最左侧维度。当有多个文件并且每个文件具有多个时间时，如果最后一个文件的时间数少于之前文件的时间数，那么剩余的数组将用缺省值填充。wrf-python中有算法会对缺省值数组进行检查，但是当你编译模块时，如果模块代码中使用了wrf-python，那么就要小心了，应尽量避免出现上述情况。

大多数情况下，timeidx 参数应设置为 wrf.ALL_TIMES。如果指定值的话，那么从每个文件中提取变量时，指定值将应用于每个文件。在具有多个时刻的多个文件中，这样做可能是没有意义的，因为每个文件的第 n 个索引可能表示不同的时刻。

通常，join 方法很少能用到。大部分情况下，只需要使用 cat 方法即可。

from __future__ import print_function
from netCDF4 import Dataset
from wrf import getvar, ALL_TIMES
wrflist = [Dataset("wrfout_d01_2016-10-07_00_00_00"),
           Dataset("wrfout_d01_2016-10-07_01_00_00"),
           Dataset("wrfout_d01_2016-10-07_02_00_00")]
# 提取所有时刻的 P 变量
p_join = getvar(wrflist, "P", timeidx=ALL_TIMES, method="join")
print(p_join)

返回(大部分信息省略)：

<xarray.DataArray u'P' (file: 3, bottom_top: 50, south_north: 1059, west_east: 1799)>

由于 Numpy 会自动对单 'Time' 维度进行自动压缩，因此，'Time' 维度被 'file' 维度替代了。如果想要维持 ‘TIme’ 维度，只需设置 squeeze 参数为 False即可。

p_join = getvar(wrflist, "P", timeidx=ALL_TIMES, method="join", squeeze=False)

<xarray.DataArray u'P' (file: 3, Time: 1, bottom_top: 50, south_north: 1059, west_east: 1799)>

• WRF 文件序列字典

wrf.getvar 函数同样可以接受字典作为参数，当使用集合时是非常有用的。然而，在字典中所有的WRF文件都应包含相同的维度。结果是一个数组，最左侧的维度是字典中的键。同样允许使用嵌套字典。

wrf_dict = {"ens1" : [Dataset("ens1/wrfout_d01_2016-10-07_00_00_00"),
                      Dataset("ens1/wrfout_d01_2016-10-07_01_00_00"),
                      Dataset("ens1/wrfout_d01_2016-10-07_02_00_00")],
            "ens2" : [Dataset("ens2/wrfout_d01_2016-10-07_00_00_00"),
                      Dataset("ens2/wrfout_d01_2016-10-07_01_00_00"),
                      Dataset("ens2/wrfout_d01_2016-10-07_02_00_00")]
p = getvar(wrf_dict, "P", timeidx=ALL_TIMES)
print(p)

返回：

<xarray.DataArray 'P' (key_0: 2, Time: 2, bottom_top: 50, south_north: 1059, west_east: 1799)>
Coordinates:
    XLONG     (south_north, west_east) float32 -122.72 -122.693 -122.666 ...
    XLAT      (south_north, west_east) float32 21.1381 21.1451 21.1521 ...
  * Time      (Time) datetime64[ns] 2016-10-07 ...
    datetime  (Time) datetime64[ns] 2016-10-07 ...
  * key_0     (key_0) <U6 u'label1' u'label2'

插值

wrf.interplevel 函数可以插值3D场到水平层上，通常是压力层或是高度层。

from wrf import getvar, interplevel
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# 提取位势高度和压力场
z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
# 计算 500 mb 位势高度
ht_500mb = interplevel(z, p, 500.)

• 垂直剖面插值

wrf.vertcross 函数可以用来创建垂直剖面图。为了定义垂直剖面，需要指定剖面的起始和终止点。当然，也可以提供中心点和角度来进行剖面。可以使用 wrf.CoordPair 对象指定起始，终止或中心点。坐标点也可以是 (x, y) 网格点或是经纬度坐标点。当使用经纬度坐标时，需要提供 netCDF文件对象或是wrf.WrfProj 对象。

垂直层也可以通过 levels 参数指定，如果未指定，将以 1% 的增量选择大约100层。

• 使用起始和终止点

from __future__ import print_function, division
from netCDF4 import Dataset
from wrf import getvar, vertcross, CoordPair
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# 获取位势高度和压力
z = getvar(ncfile, "z")
p = getvar(ncfile, "pressure")
# 定义起始和重点坐标为网格坐标
start_point = CoordPair(x=0, y=(z.shape[-2]-1)//2)
end_point = CoordPair(x=-1, y=(z.shape[-2]-1)//2)
# 计算垂直剖面。通过设置 latlon = True，将沿着剖面线计算经纬度坐标 
# 并且添加经纬度坐标到 xy_loc 元数据，从而帮助绘图
p_vert = vertcross(p, z, start_point=start_point, end_point=end_point, latlon=True)

• 使用中心点和角度

# 在网格坐标中定义中心点和角度, 中心点在网格的中心
pivot_point = CoordPair(x=(z.shape[-1]-1)//2, y=(z.shape[-2]-1)//2)
angle = 90.0
p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle, latlon=True)

• 使用经纬度坐标

lats = getvar(ncfile, "lat")
lons = getvar(ncfile, "lon")
# 使用相同的水平线，但是分别为纬度和经度
start_lat = lats[(lats.shape[-2]-1)//2, 0]
end_lat = lats[(lats.shape[-2]-1)//2, -1]
start_lon = lons[(lats.shape[-2]-1)//2, 0]
end_lon = lons[(lats.shape[-2]-1)//2, -1]
# 创建剖面线
start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)
# 使用经纬度坐标时，必须提供netcdf文件对象或投影对象
p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)

返回：

<xarray.DataArray u'pressure_cross' (vertical: 100, idx: 1798)>

• 指定垂直层

# 指定垂直层
levels = [1000., 2000., 3000.]
p_vert = vertcross(p, z, wrfin=ncfile, levels=levels, start_point=start_point, end_point=end_point, latlon=True)

返回：

<xarray.DataArray u'pressure_cross' (vertical: 3, idx: 1798)>

对比上述两个插值后返回的结果可以发现，此例中只返回3各垂直层，而使用经纬度坐标的返回了100个垂直层。

• 插值2D场到一条线

使用 wrf.interpline 函数可以沿着一条线对2D场进行插值，这类似3D场的垂直剖面插值。为了定义插值的线，可以是线的起始和终止点。当然，也可以提供中心点和角度来进行剖面。可以使用 wrf.CoordPair 对象指定起始，终止或中心点。坐标点也可以是 (x, y) 网格点或是经纬度坐标点。当使用经纬度坐标时，需要提供 netCDF文件对象或是wrf.WrfProj 对象。

• 给定起始和终止点

from wrf import getvar, interpline, CoordPair
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
# 获取 2m 温度
t2 = getvar(ncfile, "T2")
# 使用起始和终止点在区域中心创建南北线
start_point = CoordPair(x=(t2.shape[-1]-1)//2, y=0)
end_point = CoordPair(x=(t2.shape[-1]-1)//2, y=-1)
# 通过设置 latlon = True，将沿着线计算经纬度坐标 
# 并且添加经纬度坐标到 xy_loc 元数据，从而帮助绘图
t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True)

• 给定中心点和角度

# 使用中心点和角度创建南北线
pivot_point = CoordPair((t2.shape[-1]-1)//2, (t2.shape[-2]-1)//2)
angle = 0.0
t2_line = interpline(t2, pivot_point=pivot_point, angle=angle, latlon=True)

• 使用经纬度坐标

lats = getvar(ncfile, "lat")
lons = getvar(ncfile, "lon")
# 选择经纬度使其通过区域中心作为插值线
start_lat = lats[0, (lats.shape[-1]-1)//2]
end_lat = lats[-1, (lats.shape[-1]-1)//2]
start_lon = lons[0, (lons.shape[-1]-1)//2]
end_lon = lons[-1, (lons.shape[-1]-1)//2]
# 创建 CoordPairs
start_point = CoordPair(lat=start_lat, lon=start_lon)
end_point = CoordPair(lat=end_lat, lon=end_lon)
t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=Tru

• 插值3D场为面类型

wrf.vinterp 函数用于插值一个场为面类型。可用的面是 压力，位势高度，theta，theta-e。要插值的表面层同样需要指定。

from wrf import getvar, vinterp
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
tk = getvar(ncfile, "tk")
# 插值 tk 到 theta-e 层
interp_levels = [200, 300, 500, 1000]
interp_field = vinterp(ncfile,
               field=tk,
               vert_coord="eth",
               interp_levels=interp_levels,
               extrapolate=True,
               field_type="tk",
               log_p=True)

Lat/Lon 和 ij 坐标互转

wrf-python 提供了一些函数用于经纬度坐标和 xy 坐标的转换。 wrf.xy_to_ll() , wrf.xy_to_ll_proj() , wrf.ll_to_xy() , wrf.ll_to_xy_proj() 等方法。如果要转换多个点时，传递序列即可。

比如：

• 单坐标转换

from wrf import getvar, interpline, CoordPair, xy_to_ll, ll_to_xy
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
lat_lon = xy_to_ll(ncfile, 400, 200)

• 多坐标转换

lat_lon = xy_to_ll(ncfile, [400,105], [200,205])
x_y = ll_to_xy(ncfile, lat_lon[0,:], lat_lon[1,:])

添加地图

wrf-python 包含了一些辅助绘图函数，主要用于获取 cartopy，basemap，PyNgl使用的地图对象。对这三种绘图系统，当使用 xarray 时通过变量可直接确定地图对象，如果没有使用 xarray，可从 WRF 输出文件获取。

还包括直接从 xarray 切片中获取地理边界的函数。这在当你想要使用一个大区域的子集，而不想在此子集区域定义地图对象时非常有用。

• Cartopy用于变量

from wrf import getvar, get_cartopy, latlon_coords, geo_bounds
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")
slp = getvar(ncfile, "slp")
# 获取 cartopy 地图对象
cart_proj = get_cartopy(slp)
print (cart_proj)
# 获取经纬度坐标。绘图需要
lats, lons = latlon_coords(slp)
# 获取 SLP 的地理边界
bounds = geo_bounds(slp)
# 获取区域子集的地理边界
slp_subset = slp[150:250, 150:250]
slp_subset_bounds = geo_bounds(slp_subset)

• Cartopy 用于 WRF 输出文件

# 从 netcdf 文件中获取地图对象
cart_proj = get_cartopy(wrfin=ncfile)
# 从文件中获取地理边界，默认使用 XLAT, XLONG
# 提供变量名，可以获取其栅格边界
bounds = geo_bounds(wrfin=ncfile)

对于 basemap 和 pyngl 系统来说，只需将 get_cartopy 换为 get_basemap 或 get_pyngl 即可。

• 移动嵌套

当嵌套区域是移动的时候，使用 cat 方法合并多个文件后，区域边界将是时间的函数；当使用 join 方法合并多个文件后，区域边界将是文件和时间的函数。因此，当检测到多个时间或是文件时，依赖于地理边界的方法将返回对象数组而不是单个对象。

wrf.get_cartopy 获取的地图对象中并不包含地理边界信息。但可以使用 wrf.cartopy_xlim 和 wrf.cartopy_ylim 获取 matplotlib 轴边界数组（轴投影坐标）。

• 获取地理边界

from glob import glob
from netCDF4 import Dataset as nc
from wrf import getvar, ALL_TIMES, geo_bounds
# 获取所有区域2文件
wrf_filenames = glob("wrf_files/wrf_vortex_multi/wrfout_d02_*")
ncfiles = [nc(x) for x in wrf_filenames]
slp = getvar(ncfiles, "slp", timeidx=ALL_TIMES)
# 获取地理边界
bounds = geo_bounds(slp)

• 使用 cartopy

从变量中获取 cartopy 地图对象。因为cartopy 地图对象并不包含地理边界信息，因此仅返回一个 cartopy 对象。然而，如果需要轴边界，可以使用wrf.cartopy_xlim 和 wrf.cartopy_ylim 获取轴投影坐标中的移动边界数组。

from wrf import getvar, ALL_TIMES, get_cartopy, cartopy_xlim, cartopy_ylim
wrf_filenames = glob("wrf_files/wrf_vortex_multi/wrfout_d02_*")
ncfiles = [nc(x) for x in wrf_filenames]
slp = getvar(ncfiles, "slp", timeidx=ALL_TIMES)
# 获取 cartopy 地图对象
cart_proj = get_cartopy(slp)
print (cart_proj)