GRIB API学习笔记01——GRIB简介
ECWMF提供完整的计算机用户培训课程，参见
<//www.ecmwf.int/services/computing/>
其中关于GRIB API的培训资料参见：
<//www.ecmwf.int/services/computing/training/material/com_grib.html>
感谢ECWMF提供如此详细的培训资料！
本学习笔记内容主要由以上讲义内容的翻译，再加上我在单位大型机和工作电脑上的虚拟机的一些练习结果构成，省略一些我认为不太关键的部分。鉴于我糟糕的文学素养，未经润色的翻译文字不太连贯。
1. GRIB简介
1.1. GRIB概述
1.1.1. GRIB版本1与版本2
GRIB1格式
SECTION 0 Indicator
SECTION 1 Product Definition
SECTION 2 [ Grid Description ]
SECTION 3 [ Bitmap ]
SECTION 4 Binary Data
SECTION 5 End (7777)
GRIB2格式
SECTION 0 Indicator
SECTION 1 Identification
SECTION 2 [ Local Use ]
SECTION 3 Grid Definition
SECTION 4 Product Definition
SECTION 5 Data Representation
SECTION 6 Bitmap
SECTION 7 Binary Data
SECTION 8 End (7777)
1.1.2. GRIB1与GRIB2的区别
编码原则相似，但实现不同
消息结构不同
GRIB2的某些变量有更高的精度。
参数的编码不同
GRIB2对数据的描述基于模板和表格。
1.1.3. 从GRIB1迁移到GRIB2情况
略
1.2. 使用GRIB工具
1.2.1. 介绍GRIB API
提供高层访问方式，隐藏二进制细节。
使用相同的函数编解码两个GRIB版本的数据。
提供命令行工具。
提供Fortran 90、C和Python接口。
1.2.1.1. GRIB API key
使用key/value方法来访问GRIB消息。
key的设置与GRIB消息本身相关（版本、内容）
改变某些键值，可能会使其它键失效，使另一些键可用。
coded and computed keys
某些键（computed keys）根据其他键计算得出。
为某些键提供别名（aliases）
l Coded keys：
直接设置GRIB消息的字节，数值由解码该字节得到。如indicatorOfParameter。
l Computed keys：
通过组合其它键得到，提供消息的综合信息，提供一个访问复杂属性的简便方式
设置计算键会同时设定级联的所有相关键。
MARS关健词也作为Computed keys
GRIB API Keys 参考连接：
GRIB 1 keys
GRIB 2 keys
Edition independent keys
1.2.1.2. namespace
一组键的集合，包括以下几种类型：
l parameter
l time
l geography
l vertical
l statistics
l mars
GRIB API keys – parameter
两个版本的定义大不一样
<td valign="top">
<b>GRIB 2 keys</b>
</td>
<td valign="top">
discipline
</td>
<td valign="top">
parameterCategory
</td>
<td valign="top">
parameterNumber
</td>
<td valign="top">
typeOfFirstFixedSurface
</td>
<td valign="top">
scaleFactorOfFirstFixedSurface
</td>
<td valign="top">
scaledValueOfFirstFixedSurface
</td>
<td valign="top">
typeOfSecondFixedSurface
</td>
<td valign="top">
scaleFactorOfSecondFixedSurface
</td>
<td valign="top">
scaledValueOfSecondFixedSurface
</td>
<td valign="top">
productDefinitionTemplateNumber
</td>
<td valign="top">
…
</td>
也有两个版本通用的key
<td valign="top">
<b>Example value</b>
</td>
<td valign="top">
151
</td>
<td valign="top">
msl
</td>
<td valign="top">
ecmf (or 98)
</td>
<td valign="top">
Mean sea level pressure
</td>
<td valign="top">
Pa
</td>
GRIB API keys – time
预报起始时间
<td valign="top">
<b>Example values </b>
</td>
<td valign="top">
20130306 (YYYYMMDD)
</td>
<td valign="top">
0, 600, 1200, 1800
</td>
预报步长
<td valign="top">
<b>Example values </b>
</td>
<td valign="top">
instant, accum, avg, max, min, …
</td>
<td valign="top">
s, m, h, 3h, 6h, 12h, D, M, Y, 10Y, 30Y, C
</td>
<td valign="top">
0, 3, …
</td>
<td valign="top">
0, 3, ..
</td>
<td valign="top">
3-6, 6 (“startStep-endStep” , “endStep” )
</td>
预报时效
<td valign="top">
<b>Example values </b>
</td>
<td valign="top">
20130306 (YYYYMMDD)
</td>
<td valign="top">
0, 300, 1200, 1800
</td>
GRIB API Keys – MARS
包含所有MARS关键字
<td valign="top">
<b>Example values </b>
</td>
<td valign="top">
20120306 (YYYYMMDD)
</td>
<td valign="top">
0000, 0600, 1200, 1800
</td>
<td valign="top">
3, 6, 9, 12, …
</td>
<td valign="top">
od, …
</td>
<td valign="top">
oper, enfo,…
</td>
<td valign="top">
0001
</td>
<td valign="top">
an, fc, cf, pf, …
</td>
<td valign="top">
sfc, pl, ml
</td>
<td valign="top">
500, 850, …
</td>
<td valign="top">
151.128
</td>
查看GRIB文件并寻找关键字是否有效的最简单方式就是使用GRIB Tools：
grib_ls：摘要信息
grib_dump：更详细的信息
1.2.2. 检查GRIB消息的内容
1.2.2.1. 基本概念
提供完成常规GRIB数据处理的简便解决方案，推荐尽可能使用命令行工具。
基本语法：
命令概述：
提供GRIB API本身消息：
girb_info, grib_keys
查看、比较GRIB消息
grib_dump, grib_ls, grib_get, grib_get_data, grib_compare
统计数目、复制消息
grib_count, grib_copy
修改消息内容
grib_set, grib_filter
1.2.2.2. grib_dump
详细显示文件中包含的一个或多个GRIB消息内容。支持多种显示方式：
l Octet mode：WMO文档格式，支持十六进制
l Debug mode：打印所有GRIB文件中可用的key
以上两个模式无法同时使用。
可以打印Key的别名和类型信息。
用法：
参数
<td valign="top">
<b>Octet mode (WMO Documentation style) </b>
</td>
<td valign="top">
<b>Debug mode </b>
</td>
<td valign="top">
<b>Print key alias information </b>
</td>
<td valign="top">
<b>Print key type information </b>
</td>
<td valign="top">
<b>Print octet content in hexadecimal </b>
</td>
<td valign="top">
<b>Where option </b>
</td>
<td valign="top">
<b>Print all data values </b>
</td>
操作示例：
显示：
以#开头的是只读key key的名字大小写敏感。
例子2：
[bash]grib_dump -O gmf.639.2013112800003.grib2[/bash]
显示
例子3
[bash]grib_dump -OtaH gmf.639.2013112800003.grib2[/bash]
显示：
例4：debug模式，显示更多的computed keys
例5：C代码示例
[bash]>grib_dump -w count=1 -C -d rmf.gra.2013103112000.grb2 > 1.c[/bash]
得到如下的代码
[c]
include <grib_api.h>
/* This code was generated automatically */
int main(int argc,const char** argv)
{
grib_handle h = NULL;
size_t size = 0;
double vdouble = NULL;
long* vlong = NULL;
FILE* f = NULL;
const char* p = NULL;
const void* buffer = NULL;
if(argc != 2) {
fprintf(stderr,”usage: %s out\n”,argv[0]);
exit(1);
}
h = grib_handle_new_from_samples(NULL,”GRIB2″);
[/c]
1.2.2.3. grib_ls
显示GRIB文件的摘要信息。 可以设定显示的key 可以排序 可以用来寻找最接近某个经纬度的网格点，并打印该点的值。例如寻找最近的一个或四个网格点。 还可以输出一个C程序示例，生成带数据或不带数据的GRIB文件。该代码可以用作生成类似GRIB文件的模板。 用法：
基本选项：
<td valign="top">
<b>Keys to print </b>
</td>
<td valign="top">
<b>Additional keys to print </b>
</td>
<td valign="top">
<b>Where option </b>
</td>
<td valign="top">
<b>Order by: “step asc, centre desc” </b>
</td>
<td valign="top">
<b>Print keys for namespace </b>
</td>
<td valign="top">
<b>Print MARS keys </b>
</td>
<td valign="top">
<b>Print data value at given index </b>
</td>
<td valign="top">
<b>Value(s) nearest to lat-lon point </b>
</td>
<td valign="top">
<b>Format for floating point values </b>
</td>
<td valign="top">
<b>Minimum column width (default 10) </b>
</td>
例子1：基本例子
[bash]grib_ls gmf.639.2013112800003.grib2[/bash]
gmf.639.2013112800003.grib2
edition centre date dataType gridType step levelType lev short_name packingType valuesCount
2 babj 20131128 fc regular_ll 3 pl 10 gh grid_jpeg 820480
2 babj 20131128 fc regular_ll 3 pl 20 gh grid_jpeg 820480
2 babj 20131128 fc regular_ll 3 pl 50 gh grid_jpeg 820480
例子2：指定需要打印的key
例子3：不存在的key，返回not_found，程序并不异常退出。
例子4：使用where选项，可以用于所有GRIB工具中。
条件格式：
IS key=value
NOT key!=value
AND key1=value1,key2=value2
OR key=value1/value2/value3
例子5：使用grib_ls寻找最近网格点
MODE：
4 打印最近邻的4个点，默认
1 最近的一个点
file 指定一个GRIB文件作为mask，打印最近的land point(with mask>=0.5)。
1.2.2.4. GRIB Examiner
Metview组件
1.2.2.5. grib_get
从一个或多个GRIB文件中获取一个或多个键值，与grib_ls类似。
默认在错误发生（比如key没找到）的情况下，grib_get会运行失败，返回非零的退出码，适合在脚本中使用，也可以强制关闭出错运行失败。
浮点数支持C风格的格式控制符。
用法：
选项：
<td valign="top">
<b>Keys to get </b>
</td>
<td valign="top">
<b>Additional keys to get with –m, -n </b>
</td>
<td valign="top">
<b>Where option </b>
</td>
<td valign="top">
<b>Keys to set </b>
</td>
<td valign="top">
<b>Get all keys for namespace </b>
</td>
<td valign="top">
<b>Get all MARS keys </b>
</td>
<td valign="top">
<b>Value(s) nearest to lat-lon point </b>
</td>
<td valign="top">
<b>Format for floating point values </b>
</td>
<td valign="top">
<b>Do <i>not </i>fail on error </b>
</td>
例1：打印第一条消息，string和long
例2：错误退出
例3：Mars
例4：打印整个namespace的key
例5：控制输出格式
使用-F参数，后接C风格输出字符串。
默认格式为 –F “%.10e”
例6：stepRange和stepUnits
以integer格式输出，默认单位为小时。
用其他单位表示，需要用-s选项设置stepUnits。stepUnits 可以为 s, m, h, 3h, 6h, 12h, D, M,Y, 10Y, 30Y, C。
例7：使用grib_get寻找最近网格点
与grib_ls相同
grib_get -w count=1 -l 39.92,116.46 rmf.gra.2013103112000.grb2
23.5779 23.7622 23.6317 23.9018
例8：获取某个网格点的值
使用-i加序号。例如，使用grib_ls找到最近的点，并用grib_get创建该点的值列表。
1.2.2.6. grib_get_data
打印数值
使用-F控制value显示格式，不控制经纬度坐标。
默认不现实missing values
默认出错时返回错误码，可以用-f强制不出错。
用法：
选项
<td valign="top">
<b>Keys to print </b>
</td>
<td valign="top">
<b>Where clause </b>
</td>
<td valign="top">
<b>Specify missing value string </b>
</td>
<td valign="top">
<b>C-style format for output values </b>
</td>
<td valign="top">
<b>Do <i>not </i>fail on error </b>
</td>
<td valign="top">
<b>Print GRIB API Version </b>
</td>
例1：
例2：用-m设定missing value显示
1.2.2.7. grib_copy
复制grib文件内容，并打印指定key值，默认不打印key值 使用-v打印指定key。 输出可以排序 键值可以用于指定输出文件名。 发生错误时会出错，用-f取消。 用法：
参数
<td valign="top">
<b>Keys to print (only with –v) </b>
</td>
<td valign="top">
<b>Where option </b>
</td>
<td valign="top">
<b>Order by: “step asc, centre desc” </b>
</td>
<td valign="top">
<b>Verbose </b>
</td>
<td valign="top">
<b>Do <i>not </i>fail on error </b>
</td>
例1：
例2：使用-v和-p显示输出
例3：在输出文件名中使用键值
在””中使用[]标出键值。
1.2.2.8. grib_set
设置键值对，进行简单修改。 每个grib消息都写入输出文件，默认包括那些没有修改的消息。使用参数-S（strict）将只拷贝满足where子句所有条件的消息。 可以设置重新打包的方式。 在发生错误时会执行失败。 用法：
常用参数：
<td valign="top">
<b>List of key / values to set </b>
</td>
<td valign="top">
<b>Keys to print (only with –v) </b>
</td>
<td valign="top">
<b>Where option </b>
</td>
<td valign="top">
<b>Set all data values to value </b>
</td>
<td valign="top">
<b>Do <i>not </i>fail on error </b>
</td>
<td valign="top">
<b>Verbose </b>
</td>
<td valign="top">
<b>Strict </b>
</td>
<td valign="top">
<b>Repack data </b>
</td>
例1：设置参数
该句设置的效果举例：
shortName 设为10si
paramId设为207
name / parameterName设为‘10 metre wind speed’
units / parameterUnits设为‘m s ** -1’
indicatorOfParameter设为207
marsParam设为207.128
例2：某些key只读，无法更改，只能通过设置其它相关key来修改
缩放，使用scaleValuesBy乘以一个因子
在输出文件中使用键值，与grib_copy一样。
grib_set无法实现的功能：
数据转换，from spectral to grid-point representation等
更改网格
选择子区域数据
GRIB工具无法用于插值数据。
1.2.3. GRIB消息解码
1.2.4. 操控GRIB消息
1.3. 使用Fortran 90解码GRIB消息
1.3.1. GRIB API 用户接口
C：
Fortran 90：
Python：
1.3.2. 通用框架
调用GRIB_API库的代码通常包含以下几步，以Fortran为例
- 打开一个或多个GRIB文件
在Fortran中使用grib_open_file/grib_close_file打开和关闭文件
- 将一个或多个GRIB消息加载到内存
返回一个唯一的grib identifier，可以用于操控已加载的grib消息
- 调用编解码GRIB消息的函数
只有加载到内存中的消息才可以被编解码
应该只编解码需要的那部分，而不是整个消息。
- 将一个或多个GRIB消息写入文件
- 释放加载的GRIB消息
- 关闭打开的GRIB文件
1.3.3. GRIB API F90接口细节
所有subroutine均以grib_开头
所有routines有一个额外参数用于错误处理。
[fortran]
subroutine grib_new_from_file(ifile, igrib, status)
integer, intent (in) :: ifile
integer, intent (out) :: igrib
integer, optional, intent (out) :: status
[/fortran]
如果status不存在且发生错误，程序会停止并返回错误码给shell
使用status去处理错误
[fortran]
call grib_new_from_file(ifile, igrib, status)
if (status /= 0) then
call grib_get_error_string(status,err_msg)
print*,’GRIB_API Error: ’,trim(err_msg),’ (err=‘,status,’)’
call mpi_finalize(ierr)
stop
end if
[/fortran]
可以在下面的网站查看错误码的含义：
https://software.ecmwf.int/wiki/display/GRIB/Error+codes
1.3.3.1. 加载或释放GRIB消息
GRIB API只能处理加载过的GRIB消息。
三种主要方式
- grib_new_from_file (ifile, igrib)
从用grib_open_file打开的文件中加载GRIB消息，使用grib_close_file关闭文件。
- grib_new_from_samples(igrib, “GRIB1” )
从样例文件中加载GRIB文件，用于编码
- grib_new_from_index (indexid, igrib)
从index中加载GRIB消息，index需要首先被创建。
返回一个唯一的grib identifier(igrib)，使用该描述符操控消息。
无法直接访问内存中的缓冲区，隐藏在GRIB API内部。
任意GRIB消息占用的缓冲区均保留在内存中，需要使用grib_release(igrib)释放已经加载grib消息的缓冲区。
例：
[fortran]
PROGRAM load_message
USE grib_api
IMPLICIT NONE
INTEGER :: rfile, igrib
CHARACTER(LEN=256), PARAMETER :: input_file=’input.grb‘
CHARACTER(LEN=10), PARAMETER :: open_mode=’r‘
!
! Open GRIB data file for reading.
!
CALL grib_open_file(rfile, input_file, open_mode)
CALL grib_new_from_file(rfile, igrib)
CALL grib_release (igrib)
CALL grib_close_file (rfile)
END PROGRAM load_message
[/fortran]
1.3.3.2. 解码已加载的GRIB消息
无需解码整个消息，除非使用grib_dump。
解码具有统一的格式：
[fortran]
grib_get (igrib, keyname, values, status)
integer, intent (in) :: igrib
character(len=), intent (in) :: keyname
,[dimension(:),] intent (out) :: values
Integer, optional, intent (out) :: status
[/fortran]
type是integer或single/double precision real或string。
具体函数：
获取数组的大小：
[fortran]
grib_get_size (igrib, keyname, size, status)
integer, intent (in) :: igrib
character(len=), intent (in) :: keyname
integer, intent (out) :: size
integer, optional, intent (out) :: status
[/fortran]
输出GRIB消息的内容
[fortran]
grib_dump (igrib, status)
integer, intent (in) :: igrib
integer, optional, intent (out) :: status
[/fortran]
检查某个key是否missing
[fortran]
grib_is_missing (igrib, keyname, missing, status)
integer, intent(in) :: igrib
character(len=), intent (in) :: keyname
integer, intent (out) :: missing
[/fortran]
获取文件中的消息个数：
[fortran]
grib_count_in_file(ifile, count, status)
[/fortran]
例：grib_get
[fortran]
! Load all the GRIB messages contained in file.grib1
call grib_open_file(ifile, ‘file.grib1′,’r’)
call grib_new_from_file(ifile,igrib, iret)
LOOP: do while (iret /= GRIB_END_OF_FILE)
! Decode/encode data from the loaded message
call grib_get( igrib , “dataDate”, date)
call grib_get(igrib, “typeOfLevel”, levtype)
call grib_get(igrib, “level”, level)
call grib_get_size(igrib, “values”, nb_values)
allocate(values(nb_values))
call grib_get(igrib, “values”, values)
print, date, levtype, level, values(1), values(nb_values)
! Release
deallocate(values)
call grib_release(igrib)
! Next message
call grib_new_from_file(ifile,igrib, iret)
end do LOOP
call grib_close_file(ifile)
[/fortran]
1.3.3.3. GRIB API的更多功能
为检索或计算GRIB消息的额外信息提供一套高层keys和子程序。例如：
只读的平均值、最值等
计算经纬度和数值：grib_get_data
提取数值：
grib_find_nearest: 提取指定点最近的数值
grib_get_element: 从一组index中提取数值
1.3.3.4. GRIB API indexed access
子程序：
grib_index_create(indexid, filename, keys, status)
根据文件创建索引
grib_index_get_size(indexid, key, size, status)
获取索引中某个key的维度
grib_index_get(indexid, key, values, status)
获取索引中某key的不同值
grib_index_select(indexid, key, value, status)
选择索引中某key的值为value
grib_new_from_index(indexid, igrib, status)
to load the GRIB message corresponding to the selection made
对于随机存取，索引访问比顺序访问更快。
1.4. GRIB解码总结
尽可能使用GRIB Tool
使用版本无关的key
如果必须写程序，需要注意如何访问元素。
1.4.1. 相关链接
The WMO FM 92 GRIB Manuals can be obtained from
The ECMWF GRIB API manual is available at
The GRIB Tools are documented at
GRIB API Fortran 90 interface:
https://software.ecmwf.int/wiki/display/GRIB/Fortran+package+grib_api
GRIB API examples
https://software.ecmwf.int/wiki/display/GRIB/GRIB+API+examples
