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 数组，而不关心元数据。通过以下两种方式可以禁用元数据。
- 完全禁用 xarray
- 设置 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 函数按照以下流程执行：
- 如果没有缺省值或填充值，那么将直接调用 xarray.DataArray.values 属性返回值
- 如果有缺省值或填充值，那么会用 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)
# 获取 x轴的边界数组
xlims = cartopy_xlim(slp)
# 获取y轴的边界数组
ylims = cartopy_ylim(slp)对于 basemap 和 pyngl 系统直接使用 get_basemap 和 get_pyngl 即可。