NCAR放弃PyNGL后又入新坑？

bugsuse

发布于 2022-09-23 14:54:03

8150

发布于 2022-09-23 14:54:03

文章被收录于专栏：气象杂货铺

最近在NCAR的GitHub官方源中又出现了一个新的小项目:

A lightweight interface for reading in output from the Weather Research and Forecasting (WRF) model into xarray Dataset

利用Python进行WRF模式后处理的应该都知道，wrf-python用来处理WRF模式结果比较方便，但又太笨重了，经常需要编写很多代码。xarray是目前地球科学领域使用非常多的库，集成度非常高，使用非常方便。

由于一些原因，xarray直接处理WRF模式结果一直是痛点：

WRF输出的nc格式文件不是CF兼容的
wrf-python 需要和 netCDF4-python 和 xarray 接口交互
wrf-python 中缺乏接口以充分利用dask的并行能力

salem 库可以处理WRF模式结果，然后和 xarray 对象整合，但是 salem 库在处理诊断量方面不是很方便。

为了充分利用 xarray 的强大功能，NCAR近期发起了名为 xwrf 的小项目。

可视化实例

目前这个库正处于试验开发阶段，可以直接通过如下命令安装:

python -m pip install git+https://github.com/NCAR/xwrf.git

在官方说明以及示例中可以看到 xwrf 可以更方便的进行批量处理，并且可以和交互式可视化库（比如hvplot）协同工作，进一步简化数据处理和分析的流程。

示例

以下是官方提供的简单示例，可以进行批量处理并且进行交互式可视化：

import glob

import holoviews as hv
import hvplot
import hvplot.xarray
import xarray as xr
import xwrf
from distributed import Client
from ncar_jobqueue import NCARCluster

hv.extension('bokeh')

# 创建集群实例
cluster = NCARCluster()
cluster.scale(10)
client = Client(cluster)
client

files = sorted(glob.glob('/glade/scratch/bruyerec/IAG/METGRID/*.nc'))

# 选择文件
file_subset = files[-80:]

# 变量名
variables = ["PRES", "TT", "RH"]

# 预处理函数
def preprocess(ds):
    return ds[variables]

# 批量处理，选择了 xwrf 后端引擎
ds = xr.open_mfdataset(
    file_subset,
    engine="xwrf",
    parallel=True,
    concat_dim="Time",
    combine="nested",
    preprocess=preprocess,
    chunks={'Time': 1},
)

读取好数据之后，剩下的数据处理和可视化操作与常规使用 xarray 的方式类似，比如查看文件内容：

然后可以选择指定变量指定维度的数据：

ds.PRES.isel(num_metgrid_levels=0, south_north=0, west_east=0).values

也可以选择指定气压层的数据：

ds = ds.isel(num_metgrid_levels=range(1, 21))

可以重命名气压层，并进行计算：

ds_to_plot = ds.rename({'num_metgrid_levels': 'plev'})

plevs = ds_to_plot.PRES.isel(south_north=0, west_east=0).values[0, :]
ds_to_plot['plev'] = plevs / 100
ds_to_plot['plev'].attrs['units'] = 'hPa'

重命名及处理后Dataset

最后可以结合 hvplot 进行交互式绘图：

plot = ds_to_plot.hvplot.quadmesh(
    x='XLONG_M',
    y='XLAT_M',
    z=['TT', 'RH'],
    groupby=['Time', 'plev'],
    widget_location='bottom',
    rasterize=True,
    coastline=True,
    cmap='inferno',
)
plot