Meteva笔记：加载GRIB 2要素场

September 18, 2020 (最后修改: August 22, 2021)

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

使用 meb.set_griddata_coords() 函数修改坐标

meb.set_griddata_coords(
    t850_grid,
    name="t",
    member_list = ["GRAPES_GFS"],
    gtime = [start_time],
    dtime_list = [24]
)
t850_grid

使用 Meteva 内置的绘图函数 meb.tool.plot_tools.contourf_2d_grid() 绘制等值线图

meb.tool.plot_tools.contourf_2d_grid(
    t850_grid,
)

读取本地 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

转换

将 time，step 和 pl 都扩展为维度，并将单位转为摄氏度

field = field.expand_dims(["time", "step", "pl"])
field = field - 273.15
field

使用 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

裁剪成与 t850_grid 相同的网格

cropped_grid_data = grid_data.sel(
    lat=slice(-9.875, 79.88),
    lon=slice(0.0, 180.0)
)
cropped_grid_data

验证

对比本地提取的要素场和从 GDS 中获取的要素场是否相同。

从下图中可以看到，member 的名称不同，同时 time 也不同。 GDS 使用的是北京时间，而本地文件使用世界时。

修改坐标值，与 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

求两个场的偏差

diff_t850 = cropped_grid_data - t850_grid
diff_t850

求偏差场中最大偏差

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

预报场

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

计算指标

计算均方根误差 RMSE

使用 squeeze 方法删掉长度为 1 的维度，将数据变为二维矩阵

mem.rmse(
    anal_grid.squeeze(), 
    fcst_grid.squeeze(),
)

计算多个预报数据的指标

加载另一个数据：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

计算均方根误差

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 项目：

《nwpc-data 库简介》

https://github.com/nwpc-oper/nwpc-data