GRIB笔记：为GRIB 2文件创建要素场描述表格
目录
简介
本文介绍如何从 GRIB 2 文件读取元信息，创建文件内容描述表格，表格中每行代表文件中的一个 GRIB 2 要素场。
数据
本文面向 GRAPES 系列模式数据创建要素场信息描述表格。
GRAPES 系列模式的相关发文中已提供详细的预报产品清单，不过清单中给出的序号并不是要素场在 GRIB 2 文件中的序号。 虽然我们从 GRIB 2 文件中载入某个要素场也不应该依赖要素场在文件中的序号，但按要素场序号排序的清单有助于对自定义编码变量的理解。
方法
使用 ecCodes 解码 GRIB 2 格式数据。 ecCodes 定义一系列 GRIB Key 用于从 GRIB 2 消息中读取元信息，非常方便易用。
GRAPES 系列模式 GRIB 2 产品单个文件中仅包含单个时次单个时效单个预报成员的要素场，所以清单中仅需要包含与要素场本身相关的信息。
本文从读取如下三个方面的元数据：
- 要素：要素场表示什么物理量，例如温度、湿度、风等
- 垂直层次：要素场在垂直方向上的位置，例如地面、850hPa 等等
- 时间信息：要素场在时间上的类型和范围，比如要素场可能表示瞬时值 (10m 风)，或区间值统计值 (10m 最大风速)
具体读取的 GRIB Key 列表如下：
| 类型 | key | 含义 |
|---|---|---|
| 序号 | count | 序号 |
| 要素 | shortName | 简称 (可能为 unknown) |
| discipline | 学科 (均为 0) | |
| parameterCategory | 分类 | |
| parameterNumber | 变量代码 | |
| 垂直层次 | typeOfLevel | 层次类型 |
| level | 层次值 | |
| typeOfFirstFixedSurface | 第一层类型 | |
| scaleFactorOfFirstFixedSurface | 第一层因子 | |
| scaledValueOfFirstFixedSurface | 第一层值 | |
| typeOfSecondFixedSurface | 第二层类型 | |
| scaleFactorOfSecondFixedSurface | 第二层因子 | |
| scaledValueOfSecondFixedSurface | 第二层值 | |
| 时间信息 | stepRange | 时效范围 |
| stepType | 时效类型 |
实现
ecCodes 中 GRIB Key 的值可以用多种类型表示
key_type_mapper = {
"str": str,
"int": int,
"float": float,
"": None,
}
将上述 Key 放到列表中，使用 : 定义值类型
keys = [
# 序号
nt",
# 要素
"shortName",
"discipline",
"parameterCategory",
"parameterNumber",
# 垂直层次信息
"typeOfLevel",
"level",
"typeOfFirstFixedSurface:str",
"scaleFactorOfFirstFixedSurface",
"scaledValueOfFirstFixedSurface",
"typeOfSecondFixedSurface:str",
"scaleFactorOfSecondFixedSurface",
"scaledValueOfSecondFixedSurface",
# 时间信息
"stepRange",
"stepType"
]
从文件中载入所有值，其中使用 eccodes.codes_is_missing() 检查值是否未定义。
values = []
f = open(file_path, "rb")
gid = eccodes.codes_grib_new_from_file(f)
while gid is not None:
record = []
for key in keys:
token = key.split(":")
key_name = token[0]
key_type = None
if len(token) == 2:
key_type = key_type_mapper[token[1]]
value = eccodes.codes_get(gid, key_name, key_type)
# check whether value is MISSING
flag = eccodes.codes_is_missing(gid, key_name)
if flag:
value = np.NaN
record.append(value)
eccodes.codes_release(gid)
values.append(record)
gid = eccodes.codes_grib_new_from_file(f)
将结果表示成表格，并计算两个层次的值
df = pd.DataFrame(data=values, columns=keys)
df["first_level"] = np.power(10, (-df["scaleFactorOfFirstFixedSurface"])) * df["scaledValueOfFirstFixedSurface"]
df["second_level"] = np.power(10, (-df["scaleFactorOfSecondFixedSurface"])) * df["scaledValueOfSecondFixedSurface"]
将上述过程封装成一个函数 read_file_content()
GRAPES GFS 的调用示例，将表格保存为 CSV 文件：
from nwpc_data.data_finder import find_local_file
system = "grapes_gfs_gmf"
name = "orig"
forecast_time = "3h"
file_path = find_local_file(
f"{system}/grib2/{name}",
start_time=pd.to_datetime("2021-05-06 00:00:00"),
forecast_time=forecast_time
)
print(file_path)
df = read_file_content(file_path)
print(df)
df.to_csv(f"{system}-{name}-{forecast_time}.csv", index=False)
结果
GRAPES MESO 3KM 原始分辨率 GRIB 2 产品生成的清单如下所示：
注：“2月3日”是 excel 解析错误，实际上是 “2-3”
