近几年,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.getvar函数也可以用来提取常规的WRF输出的netCDF 变量。
p = getvar(ncfile, "P")
有时候你只需要返回常规的 numpy 数组,而不关心元数据。通过以下两种方式可以禁用元数据。
# 方法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))
如果你需要将 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 函数按照以下流程执行:
# 获取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 方法会将序列中所有文件沿着 '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方法合并一系列文件时,会将文件/序列索引作为新数组的最左侧维度。当有多个文件并且每个文件具有多个时间时,如果最后一个文件的时间数少于之前文件的时间数,那么剩余的数组将用缺省值填充。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.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个垂直层。
使用 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
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() 等方法。如果要转换多个点时,传递序列即可。
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 切片中获取地理边界的函数。这在当你想要使用一个大区域的子集,而不想在此子集区域定义地图对象时非常有用。
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)
# 从 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 对象。然而,如果需要轴边界,可以使用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)