使用cdo将GrADS数据转换成Grib2数据
虽然 GRAPES 模式可以直接输出 GRIB 2 格式数据，但 NWPC 的模式试验通常选择输出 GrADS 二进制格式。 当需要与其他模式数据进行比较时，就无法避免需要将 GrADS 格式转成 GRIB2 格式。
虽然 NWPC 已开发一些转码工具（例如 grads2grib.exe等），但大多数工具仅提供二进制可执行程序，没有源码的工具很难推广应用，至少我个人不希望将自己的工作构建在非开源软件之上。
所以，我也一直在寻找一款合适的转码工具。
cdo 提供强大的数据处理功能，支持常见的数据格式，也支持将 GrADS 数据转换成 Grib2 格式。
本文介绍使用 cdo 将 GrADS 数据转换成 Grib2 数据的简单方法。
开始
执行一行命令就可以实现转码。
下面的语句将 GRAPES GFS 输出的地面场和等压面场二进制数据 postvar 转码为 Grib 2 文件。
cdo -v -f grb2 import_binary post.ctl_2020070900_000 a.grb2使用 grib_ls 查看转码后的文件内容：
grib_ls a.grb2以下省略部分输出内容。
a.grb2
edition centre date dataType gridType stepRange typeOfLevel level shortName packingType
2 ecmf 20200709 af regular_ll 0 isobaricInhPa 10 u grid_ieee
2 ecmf 20200709 af regular_ll 0 isobaricInhPa 9 u grid_ieee
...
2 ecmf 20200709 af regular_ll 0 isobaricInPa 20 u grid_ieee
2 ecmf 20200709 af regular_ll 0 isobaricInPa 10 u grid_ieee
2 ecmf 20200709 af regular_ll 0 isobaricInhPa 10 v grid_ieee
...
2 ecmf 20200709 af regular_ll 0 isobaricInPa 10 v grid_ieee
2 ecmf 20200709 af regular_ll 0 isobaricInhPa 10 t grid_ieee
...
2 ecmf 20200709 af regular_ll 0 surface 0 unknown grid_ieee
2 ecmf 20200709 af regular_ll 0 surface 0 unknown grid_ieee
2 ecmf 20200709 af regular_ll 0 surface 0 unknown grid_ieee
257 of 257 grib messages in a.grb2
257 of 257 total grib messages in 1 files
可以看到，常用的变量，如 u，v，t 等已经自动识别，但某些变量无法自动识别。
默认配置的问题
实际上，在运行上面的 cdo 命令时，会输出类似下面的信息：
OpenMP: num_procs=32 max_threads=1
status 0
cdo import_binary: Opening file: postvar2020070900_000
cdo import_binary: Reading timestep: 1 2020-07-09 00:00:00
gribapiDefLevel : Changed zaxis type from generic to pressure
...
GRIB_API ERROR : concept: no match for shortName=u10m
!!! grib_set_string( grib_handle* h, "shortName", "u10m") !!!
Warning (gribapiDefParam): grib_api: No match for shortName=u10m
gribapiEncode : *** GRIB2 shortName does not correspond to chosen variable name: "unknown" ("u10m").
GRIB_API ERROR : concept: no match for shortName=v10m
!!! grib_set_string( grib_handle* h, "shortName", "v10m") !!!
Warning (gribapiDefParam): grib_api: No match for shortName=v10m
gribapiEncode : *** GRIB2 shortName does not correspond to chosen variable name: "unknown" ("v10m").
GRIB_API ERROR : concept: no match for shortName=pblh
!!! grib_set_string( grib_handle* h, "shortName", "pblh") !!!
Warning (gribapiDefParam): grib_api: No match for shortName=pblh
gribapiEncode : *** GRIB2 shortName does not correspond to chosen variable name: "unknown" ("pblh").
cdo import_binary: Processed 59 variables [3.72s 296MB]可以看到，cdo 内部使用 GRIB API 进行 GRIB 2 编码，并将 GrADS 数据 ctl 文件中的变量名当成 GRIB API 中的 shortName。
GRAPES 模式的 postvar 文件中部分变量名，如 q2，q3 等，不是内置的 shortName，无法被自动识别。
当然，GRIB API 支持添加自定义的 shortName。详细方法请参考之前的文章《GRIB学习笔记：添加自定义表格》。
不过使用自定义的表格同样存在无法推广应用的问题。
cdo 提供额外的机制，可以自定义输入与输出变量的对应关系。
自定义参数表
cdo 提供大量修改数据信息的方法。
其中 SETPARTAB 系列方法 setpartabp 和 setpartabn 用于设置参数表。
语法
< operator >,table[,convert] infile outfileoperator 是 setpartabp 或 setpartabn。
setpartabp使用变量标识符查找变量setpartabn使用变量名查找变量
table 是参数表文件。
convert 用于指示转换单位。
参数表
参数表是一个文本文件，由多个如下所示的参数描述单元构成。语法类似 Fortran 语言的 namelist 文件。
¶meter
name = t
out_name = ta
standard_name = air_temperature
units = "K"
missing_value = 1e+20
valid_min = 157.1
valid_max = 336.3
/字段如下
| Entry | Type | Description |
|---|---|---|
| name | WOR | Name of the variable |
| out_name | WORD | New name of the variable |
| param | WORD | Parameter identifier (GRIB1: code[.tabnum]; GRIB2: num[.cat[.dis]]) |
| out_param | WORD | New parameter identifier |
| type | WORD | Data type (real or double) |
| standard_name | WORD | As defined in the CF standard name table |
| long_name | STRING | Describing the variable |
| units | STRING | Specifying the units for the variable |
| comment | STRING | Information concerning the variable |
| cell_methods | STRING | Information concerning calculation of means or climatologies |
| cell_measures | STRING | Indicates the names of the variables containing cell areas and volumes |
| missing_value | FLOAT | Specifying how missing data will be identified |
| valid_min | FLOAT | Minimum valid value |
| valid_max | FLOAT | Maximum valid value |
| ok_min_mean_abs | FLOAT | Minimum absolute mean |
| ok_max_mean_abs | FLOAT | Maximum absolute mean |
| factor | FLOAT | Scale factor |
| delete | INTEGER | Set to 1 to delete variable |
| convert | INTEGER | Set to 1 to convert the unit if necessary |
设置 shortName
ctl 中的某些变量在 GRIB API 的设置中已有 shortName，只不过与 ctl 中的名字不同。
例如 10 米风在 ctl 中为 u10m 和 v10m，但在 GRIB API 中为 10u 和 10v。
使用下面的参数表
¶meter
name = u10m
out_name = 10u
/
¶meter
name = v10m
out_name = 10v
/使用下面的命令转码
cdo -v -f grb2 -setpartabn,wind.param,convert \
-import_binary post.ctl_2020070900_000 b.grb2使用 grib_ls 查看新生成的文件内容
....
2 ecmf 20200709 af regular_ll 0 surface 0 unknown grid_ieee
2 ecmf 20200709 af regular_ll 0 surface 0 u grid_ieee
2 ecmf 20200709 af regular_ll 0 surface 0 v grid_ieee
2 ecmf 20200709 af regular_ll 0 surface 0 unknown grid_ieee
257 of 257 grib messages in b.grb2可以看到最后两个变量已被识别为 u 和 v。
但层次类型依然没有识别正确，可能是我用的方法不对。
设置变量标识符
如果 postvar 的某些变量在 GRIB API 中没有 shortName，需要直接定义变量标识符。
下面以云水混合比（q2）为例说明。
注：实际上 GRIB API 中 clwmr 表示云水混合比。
¶meter
name = q2
out_param = 22.1.0
/使用下面的命令转码
cdo -v -f grb2 -setpartabn,q2.param,convert \
-import_binary post.ctl_2020070900_000 c.grb2使用 grib_ls 查看新生成的文件内容
....
2 ecmf 20200709 af regular_ll 0 isobaricInPa 10 h grid_ieee
2 ecmf 20200709 af regular_ll 0 isobaricInhPa 10 clwmr grid_ieee
2 ecmf 20200709 af regular_ll 0 isobaricInhPa 9 clwmr grid_ieee
...可以看到位势高度（h）后面的 q2 被正确地识别为 clwmr。
讨论
cdo 是一个功能强大的命令行工具，基本能覆盖常用的数据处理功能。 当我们想编程实现某种数据操作时，需要考虑运行效率是否优于 cdo 之类的现有命令行工具。 如果使用 cdo 能满足应用需求，往往就不用继续开发新的工具。
cdo 更多用于处理 NetCDF 格式的数据，感觉对 GrADS 格式的数据支持不太好。
不过，NWPC 正在致力于将 GRAPES 模式输出改为 GRIB 2 格式，也许明年会正式进入业务系统。 所以，对于非标准的 GrADS 数据格式，也就没有继续研究的必要。
昨天 NMIC 的国际交流报告会中提到 CMA 目前在数据格式方面缺乏与国际标准的对接，缺乏标准制定上的中国声音和中国方案。
确实如此，至少 NWPC 就缺乏公开透明的 GRIB 2 编码表，更不用说提供像 NCEP WMO GRIB2 Documentation 之类的文档，或者像 ECMWF Parameter database 之类的在线查询工具。
也许这是因为数据管理不是 NWPC 的职责，只能寄希望于其他单位能拿出类似的工具，给 CMA 带来更多开源共享的精神。
参考
https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-2210002.6.2
