xarray系列｜WRF模式前处理和后处理
距离上次xarray的更新已经过去两个多星期了...，关于xarray插值方法的介绍官方文档已经给的比较详细了，也有公众号推送过相关文章 xarray指南：插值 基于xarray的气象场站点和格点插值，所以xarray的插值部分就不单独说了。
这一篇主要来说一下WRF模式的前处理和后处理部分，后处理分为：数据提取、投影转换、插值和可视化。
- WRF模式前处理
- WRF模式后处理
- 数据提取
- 投影转换
- 插值
- 可视化
本文除了xarray之外，主要使用了 salem 和 xesmf 这两个库，salem 主要是进行前处理和部分后处理操作，xesmf 主要是进行投影转换。以下是正文：
WRF模式前处理
这里所说的前处理就是通常指的WPS操作，即确定模拟域及可视化。
这部分通常都是使用WPS提供的NCL脚本来完成，但这里我们使用Python来实现。不需要从零编写脚本，只需要通过 pip install salem 安装 salem 即可。
import salem
import shapely.geometry as shpg
g, maps = salem.geogrid_simulator('namelist.wps')
fig, ax = plt.subplots()
p = shpg.point.Point([105.2, 38.3])
maps[0].set_geometry(p, crs=g[0].proj, zorder=1, linewidth=1, marker='*', facecolor='red')
maps[0].set_rgb(natural_earth='lr')
s = maps[0].visualize(title='WPS Domain')
WRF模式嵌套模拟域
WRF模式后处理
后处理的内容通常取决于你需要分析什么。这里就数据提取、投影转换、插值和可视化几个部分说一下。
由于WRF模式的输出并不完全兼容NetCDF格式的CF标准，所以无法直接利用 xarray 的很多函数。这里同样需要用到 salem 来进行转换。
所以这里读取数据的时候需要先用 netCDF4 读取，然后 salem 进行转换。
import netCDF4 as nc
wrfin = nc.Dataset('wrfout_d01')
ds = salem.open_wrf_dataset('wrfout_d01')
转换之后就可以利用 xarray 的函数进行处理和分析了。
数据提取
数据提取和之前说的类似，主要是利用 .sel 和 .isel 等函数。这里还是以提取站点数据为例，强调一个数据提取需要注意的问题。
由于WRF的坐标问题，所以这里不能直接通过经纬度选择，需要将经纬度转换为对应的索引。
sta = pd.read_csv('stations.csv')
xy = wrf.ll_to_xy(wrfin, sta.lat, sta.lon)
var_sta = ds.isel(south_north=xy[0], west_east=xy[1])
注意：这里要注意，这里的
xy是DataArray对象，xy[0]和xy[1]也是DataArray对象，得到的才是对应的站点信息。如果改为xy[0].data和xy[1].data得到的将是站点数x站点数的网格信息。这是因为在提取站点信息时，.sel这些函数接受的参数应该是DataArray对象。
投影转换
一般情况下是不需要进行投影转换的，除非在需要和其它投影的数据进行对比分析。这里我们使用 xesmf 进行网格的转换。
import xesmf as xe
## 创建经纬度网格
ds_out = xe.util.grid_2d(112, 121.01, 0.01, 36, 43.01, 0.01).rename_dims({'y': 'south_north', 'x': 'west_east'})
ds_wrf = xr.open_dataset('wrfout_d01')[['T2']].rename({'XLONG': 'lon', 'XLAT': 'lat'}).isel(Time=0)
## 创建转换器
regridder = xe.Regridder(ds_wrf, ds_out, 'bilinear')
## 转换投影
ds_out = regridder(ds_wrf)
这里需要注意：如果分辨率和范围比较大，创建投影的转换器时会非常耗时，如果需要经常使用这种固定的转换操作，可以存储权重以便重复使用。
proj = ccrs.LambertConformal(central_latitude=39.3, central_longitude=116.48,
standard_parallels = (30, 60))
fig, ax = plt.subplots(figsize=(6, 4), subplot_kw=dict(projection=proj))
ds_wrf['T2'].plot(vmin=250, vmax=275, ax=ax, x='lon', y='lat')
fig, ax = plt.subplots(figsize=(6, 4), subplot_kw=dict(projection=ccrs.PlateCarree()))
ds_out['T2'].plot(vmin=250, vmax=275, x='lon', y='lat', ax=ax)
投影转换前(左)后(右)的2m温度分布
上述仅提供了单个变量的转换示例。如果要应用到其它变量还需要进行调整。这部分算是很久以前立的flag，看这篇 网格气象场插值-NCL版，拖了一年多这次就在这里填上了...
插值
上述投影转换也可以理解成是一种插值。这里也可以使用 xarray 自带的插值方法进行插值，或者使用 salem 提供的函数进行插值，比如 .wrf_zlevel 进行垂直插值：
ds.isel(time=1).salem.wrf_zlevel('WS', levels=100)
salem 也提供了快速可视化的函数，如下：
ds.isel(time=1).salem.wrf_zlevel('WS', levels=100).salem.quick_map()
可视化
可视化的部分其实之前就提到过了，xarray 的 Dataset 和 DataArray 都有 plot 方法可以进行快速绘图，也可以非常方便的绘制多幅子图。比如：
p = ds['T2'].plot.pcolormesh(col_wrap=4, x='lon', y='lat', col='time')
# p.axes.set_xlabel('Longitude(degree)')
# p.axes.set_ylabel('Latitude(degree)')
p.fig.savefig('t2.png', dpi=300, bbox_inches='tight')
所有时刻的2m温度分布图（点击看大图）
除了这种一键可视化之外，也可以进行单个时刻的绘图，或者提取某一个站点的数据绘制时间序列图：
ds['T2'].isel(south_north=120, west_east=50).plot(color='red', linewidth=1.5, label='Temperature 2m')
plt.legend()
plt.savefig('t2.png', dpi=300, bbox_inches='tight')
单站点时间序列图
这次就说这些，salem 的后处理功能没有 wrf-python丰富，尤其是一些诊断变量和绘图的功能，但是目前wrf-python还没有提供 xarray 的兼容接口，很难利用其 xarray 很多便利的函数。每个库都有各自的优势，发挥优势才能更好的提高效率。
目前打算把平时处理WRF模式的脚本合并成命令行工具，以便平时进行快速数据处理和可视化。后续添加完成后会开源，不知道大家平时都有哪些处理操作是经常需要用到的，可以考虑一起加进去，欢迎留言提出，也可以点赞、在看给个鼓励快速开源，也可以赞赏加个鸡腿给更新加个速~
—END—
腾讯云开发者
