Meteva笔记:加载GRIB 2要素场
Meteva 是由 nmc 开源的全流程检验程序库,提供了常用的各种气象预报检验评估的算法函数,气象检验分析的图片和表格型产品的制作函数,以及检验评估系统示例。
本文介绍如何通过 nwpc-data 库将本地 GRIB 2 文件接入到 Meteva 工具中。
准备
本文代码均在 CMA-PI 高性能计算机上运行。
载入需要的库
import pandas as pd
import numpy as np
import xarray as xr
import meteva.base as meb
import meteva.method as mem
from nwpc_data.grib.eccodes import load_field_from_file
from nwpc_data.data_finder import find_local_file
meb
提供 IO 和绘图相关函数,mem
提供计算检验指标的函数。
从 GDS 加载数据
在加载本地数据文件前,首先使用 Meteva 内置的函数从 GDS 服务中获取要素场,用于后续对比验证。
GDS 服务的相关配置方法请访问同样由 nmcdev 开源的 nmcdev/nmc_met_io 项目。
Meteva 支持 nmc_met_io 项目的配置文件。 本文使用已配置好的文件:
gds_config_file = "/g1/u/USER/.config/.nmcdev/config.ini"
使用 meb.io.read_gds_ip_port
从配置文件中读取 IP 地址和端口号
ip, port = meb.io.read_gds_ip_port(gds_config_file)
格点数据
加载 GRAPES GFS 的格点数据。
构建数据路径
start_time = pd.Timestamp("2020-09-16 08:00:00")
start_time
Timestamp('2020-09-16 08:00:00')
path = meb.get_path("GRAPES_GFS/TMP/850/YYMMDDHH.TTT", start_time, 24)
path
'GRAPES_GFS/TMP/850/20091608.024'
使用 meb.read_griddata_from_gds()
函数从 GDS 服务中获取数据,返回 xr.DataArray
对象
t850_grid = meb.read_griddata_from_gds(ip, port, path)
t850_grid
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/t850-grid.png)
使用 meb.set_griddata_coords()
函数修改坐标
meb.set_griddata_coords(
t850_grid,
name="t",
member_list = ["GRAPES_GFS"],
gtime = [start_time],
dtime_list = [24]
)
t850_grid
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/t850-grid-coords.png)
使用 Meteva 内置的绘图函数 meb.tool.plot_tools.contourf_2d_grid()
绘制等值线图
meb.tool.plot_tools.contourf_2d_grid(
t850_grid,
)
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/t850-grid-plot.png)
读取本地 GRIB 2 数据
载入
文件路径
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091600",
forecast_time="24h",
)
file_path
PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020091600/ORIG/gmf.gra.2020091600024.grb2')
读取 850hPa 温度场
单个要素场默认只有两个维度:latitude 和 longitude
field = load_field_from_file(
file_path,
parameter="t",
level_type="pl",
level=850
)
field
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/field.png)
转换
将 time,step 和 pl 都扩展为维度,并将单位转为摄氏度
field = field.expand_dims(["time", "step", "pl"])
field = field - 273.15
field
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/field-trans.png)
使用 meb.xarray_to_griddata()
函数将要素场对象转为 meb.grid_data()
函数生成的 xr.DataArray
对象
可以看到,对于单个要素场,该函数自动生成了 memeber
维度,坐标值为 0
grid_data = meb.xarray_to_griddata(
field,
level_dim="pl",
time_dim="time",
dtime_dim="step",
lat_dim="latitude",
lon_dim="longitude",
)
grid_data
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/grid-data.png)
裁剪成与 t850_grid
相同的网格
cropped_grid_data = grid_data.sel(
lat=slice(-9.875, 79.88),
lon=slice(0.0, 180.0)
)
cropped_grid_data
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/cropped-grid-data.png)
验证
对比本地提取的要素场和从 GDS 中获取的要素场是否相同。
从下图中可以看到,member
的名称不同,同时 time
也不同。
GDS 使用的是北京时间,而本地文件使用世界时。
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/data-vs.png)
修改坐标值,与 t850_grid
保持一致:
- 添加 member 名称
- 将时间从世界时改为北京时
meb.set_griddata_coords(
cropped_grid_data,
member_list=["GRAPES_GFS"],
gtime=cropped_grid_data.time + pd.Timedelta(hours=8)
)
cropped_grid_data
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/cropped-grid-data-mod.png)
求两个场的偏差
diff_t850 = cropped_grid_data - t850_grid
diff_t850
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/diff.png)
求偏差场中最大偏差
abs(diff_t850).max()
0.01001473
差值可能是因为压缩精度的问题,在可以接受的范围内。 说明本地读取的 GRIB 2 文件可以代替 GDS 中的数据。
计算
计算 024 时效与该时刻分析场的均方根误差
载入数据
整合函数,实现如下功能:
- 使用 nwpc-data 从 GRIB 2 文件中加载要素场
- 将返回的要素场转换为
xr.DataArray
对象
def load_grid(
file_path,
parameter,
level_type=None,
level=None,
) -> xr.DataArray:
field = load_field_from_file(
file_path,
parameter=parameter,
level_type=level_type,
level=level,
)
field = field.expand_dims(["time", "step", "pl"])
grid = meb.xarray_to_griddata(
field,
level_dim="pl",
time_dim="time",
dtime_dim="step",
lat_dim="latitude",
lon_dim="longitude",
)
return grid
分析场
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091700",
forecast_time="0h",
)
file_path
PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020091700/ORIG/gmf.gra.2020091700000.grb2')
anal_grid = load_grid(
file_path,
parameter="t",
level_type="pl",
level=850
)
anal_grid
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/anal-grid.png)
预报场
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091600",
forecast_time="24h",
)
file_path
PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020091600/ORIG/gmf.gra.2020091600024.grb2')
fcst_grid = load_grid(
file_path,
parameter="t",
level_type="pl",
level=850
)
fcst_grid
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/fcst-grid.png)
计算指标
计算均方根误差 RMSE
使用 squeeze
方法删掉长度为 1 的维度,将数据变为二维矩阵
mem.rmse(
anal_grid.squeeze(),
fcst_grid.squeeze(),
)
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/rmse.png)
计算多个预报数据的指标
加载另一个数据:48 小时预报
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091500",
forecast_time="48h",
)
fcst_grid_48 = load_grid(
file_path,
parameter="t",
level_type="pl",
level=850
)
两个数据仅预报时效不同,沿 dtime
维度合并两个要素场
fcst_grids = xr.concat(
[
grid.squeeze().reset_coords(["member", "level", "time"], drop=True)
for grid in [fcst_grid, fcst_grid_48]
],
dim="dtime"
)
fcst_grids
![](https://cdn.perillaroc.wang/image/blog/2020/verify/meteva/io/grib2/fcst-grids.png)
计算均方根误差
mem.rmse(
anal_grid.squeeze(),
fcst_grids.values,
)
array([1.13974815, 1.85073222])
参考
Meteva 项目:
https://github.com/nmcdev/meteva
nmc_met_io 项目:
https://github.com/nmcdev/nmc_met_io
nwpc-data 项目: