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

可以看到,常用的变量,如 uvt 等已经自动识别,但某些变量无法自动识别。

默认配置的问题

实际上,在运行上面的 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 文件中部分变量名,如 q2q3 等,不是内置的 shortName,无法被自动识别。

当然,GRIB API 支持添加自定义的 shortName。详细方法请参考之前的文章《GRIB学习笔记:添加自定义表格》。

不过使用自定义的表格同样存在无法推广应用的问题。

cdo 提供额外的机制,可以自定义输入与输出变量的对应关系。

自定义参数表

cdo 提供大量修改数据信息的方法。 其中 SETPARTAB 系列方法 setpartabpsetpartabn 用于设置参数表。

语法

< operator >,table[,convert]  infile outfile

operatorsetpartabpsetpartabn

  • setpartabp 使用变量标识符查找变量
  • setpartabn 使用变量名查找变量

table 是参数表文件。

convert 用于指示转换单位。

参数表

参数表是一个文本文件,由多个如下所示的参数描述单元构成。语法类似 Fortran 语言的 namelist 文件。

&parameter 
 name         = t 
 out_name     = ta 
 standard_name = air_temperature 
 units        = "K" 
 missing_value = 1e+20 
 valid_min    = 157.1 
 valid_max    = 336.3 
/

字段如下

EntryTypeDescription
nameWORName of the variable
out_nameWORDNew name of the variable
paramWORDParameter identifier (GRIB1: code[.tabnum]; GRIB2: num[.cat[.dis]])
out_paramWORDNew parameter identifier
typeWORDData type (real or double)
standard_nameWORDAs defined in the CF standard name table
long_nameSTRINGDescribing the variable
unitsSTRINGSpecifying the units for the variable
commentSTRINGInformation concerning the variable
cell_methodsSTRINGInformation concerning calculation of means or climatologies
cell_measuresSTRINGIndicates the names of the variables containing cell areas and volumes
missing_valueFLOATSpecifying how missing data will be identified
valid_minFLOATMinimum valid value
valid_maxFLOATMaximum valid value
ok_min_mean_absFLOATMinimum absolute mean
ok_max_mean_absFLOATMaximum absolute mean
factorFLOATScale factor
deleteINTEGERSet to 1 to delete variable
convertINTEGERSet to 1 to convert the unit if necessary

设置 shortName

ctl 中的某些变量在 GRIB API 的设置中已有 shortName,只不过与 ctl 中的名字不同。 例如 10 米风在 ctl 中为 u10mv10m,但在 GRIB API 中为 10u10v

使用下面的参数表

&parameter 
 name         = u10m 
 out_name     = 10u
/
&parameter 
 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

可以看到最后两个变量已被识别为 uv

但层次类型依然没有识别正确,可能是我用的方法不对。

设置变量标识符

如果 postvar 的某些变量在 GRIB API 中没有 shortName,需要直接定义变量标识符。

下面以云水混合比(q2)为例说明。

注:实际上 GRIB API 中 clwmr 表示云水混合比。

&parameter 
 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