GRIB笔记：从GRIB 2文件中加载单个要素场

meteorology/data
目录

之前的文章《GRIB笔记：使用cfgrib加载GRIB文件》介绍如何使用 cfgrib 库读取 GRIB 2 文件。 即使只有一个要素场符合过滤条件，返回的对象仍然是xarray.Dataset，需要经过一层索引才能得到单个要素场的xarray.DataArray。

本文介绍如何使用过滤条件从 GRIB 2 文件中获取单个要素场。 示例代码来自nwpc-oper/nwpc-data项目。

过滤条件

nwpc-data 项目仅面向 GRAPES 系列模式生成的 GRIB 2 数据，单个文件中只保存同一起报时次同一预报时效的要素场。 所以想要从该文件中检索某个要素场最多需要三个限制条件：

  • 要素名
  • 层次类型
  • 层次值

例如检索 850 hPa 温度场需要设置以下三个条件：

  • shortName = "t"
  • typeOfLevel = "isobaricInhPa"
  • level = 850

如果检索2m温度只需要设置shortName = "2t"即可。

nwpc-oper/nwpc-data项目设计了一个从本地文件中检索单个要素场的API，返回xarray.DataArray对象。

def load_field_from_file(
        file_path: str or Path,
        parameter: str or typing.Dict,
        level_type: str or typing.Dict = None,
        level: int = None,
        engine: str = "cfgrib",
        **kwargs
) -> xr.DataArray or None:
    # ..skip...

该函数通过 engine 参数支持使用 cfgrib 或 eccodes-python 库加载要素场。

下面分别介绍如何使用这两个库过滤 GRIB 2 文件。

使用cfgrib筛选要素场

在 backend_kwargs 中的 filter_by_keys 字典对象中设置过滤条件。

例如检索 850 hPa 温度场需要设置：

backend_kwargs["filter_by_keys"] = {
    "shortName": "t",
    "typeOfLevel": "isobaricInhPa",
    "level": 850,
}

使用 xarray.open_dataset() 函数加载 GRIB 2 文件：

data_set = xr.open_dataset(
    file_path,
    engine="cfgrib",
    backend_kwargs=backend_kwargs
)

如本文最开始提到的文章中介绍的，该函数返回的是 xarray.Dataset 对象。 对于温度场来说，我们可以直接用 data_set.t 或者 data_set["t"] 获得该要素场对应的 xarray.DataArray 对象。 但如果检索的要素没有shortName，也就没有固定的名称，可以使用下面的函数从 Dataset 中返回第一个变量。

def _load_first_variable(data_set: xr.Dataset) -> xr.DataArray:
    first_variable_name = list(data_set.data_vars)[0]
    return data_set[first_variable_name]

使用ecCodes筛选要素场

ecCodes 没有直接内置用于过滤的函数（index可以实现类似的功能，但我更关注单个的要素场），需要依次读取文件中的要素场，并根据筛选条件从 GRIB 消息中读取属性值进行对比。

对于 850 hPa 温度场，需要从 GRIB 消息 message_id 中获取三个属性值：

short_name = eccodes.codes_get(message_id, "shortName")
message_type = eccodes.codes_get(message_id, "typeOfLevel", ktype=str)
message_level = eccodes.codes_get(message_id, "level", ktype=int)

load_field_from_file 会找到第一个符合条件的 GRIB 消息，并使用文章《从ecCodes的GRIB2消息构建xarray对象》中介绍的方法，将其封装为 xarray.DataArray 对象。

高级筛选条件

上述三个 GRIB key 仅能检索文件中的部分要素场。 想要检索文件中的所有要素场，需要添加额外的信息。

要素名

之前的文章《GRIB学习笔记：添加自定义表格》中介绍过，GRAPES 模式的 GRIB 2 产品有一些要素场使用数值预报中心自定义的值，使用默认的 ecCodes 解码时没有 shortName，无法直接获取要素名称。

GRIB 2 消息中的要素通过 disicpline, parameterCategory 和 parameterNumber 三个属性唯一确定，所以对于自定义的要素，可以使用这三个属性检索。

load_field_from_file 中的 parameter 参数支持类似下面的字典对象：

parameter = {
    "discipline": 0,
    "parameterCategory": 2,
    "parameterNumber": 225,
}

因为 cfgrib 不会读取上面三个 GRIB Key 的值，所以需要在 backend_kwargs 中的 read_keys 列表对象中加上这三个值。

backend_wargs["read_keys"].extend(["discipline", "parameterCategory", "parameterNumber"])

层次类型

ecCodes中的 typeOfLevel 有比较复杂的定义，而我们常用的层次类型包括等压面层（pl），地面场（sfc）和模式面层（ml）。

load_field_from_file 支持使用上述简化的层次类型，使用下面的函数将其转换为 ecCodes 可以识别的 GRIB Key。

def fix_level_type(level_type: str or typing.Dict or None) -> str or typing.Dict:
    if level_type is None:
        return level_type
    if isinstance(level_type, dict):
        return level_type
    if level_type == "pl":
        return {
            "typeOfLevel": "isobaricInhPa"
        }
    elif level_type == "sfc":
        return {
            "typeOfLevel": "sfc"
        }
    elif level_type == "ml":
        return {
            "typeOfFirstFixedSurface": 131,
        }

GRAPES 系列模式的模式面 GRIB 文件中，使用 typeOfFirstFixedSurface=131 作为层次类型。 经过测试发现，设置typeOfFirstFixedSurface后无法使用 cfgrib 加载。 不过，因为 GRAPES 模式的模式面 GRIB 文件中只有一种层次类型，所以从文件中检索时可以不用设置 typeOfLevel，只需要提供要素名和层次值即可。

参考

nwpc-oper/nwpc-data