ecCodes 学习 利用ecCodes fortran90 api对GRIB文件进行读写

参考 https://www.ecmwf.int/assets/elearning/eccodes/eccodes2/story_html5.html
https://confluence.ecmwf.int/display/OPTR/ecCodes%3A+GRIB+data+decoding+and+encoding+software+2018

基本解码流程

1. 指定打开方式（“读”或“写”），打开一个或多个GRIB文件;

2. 根据不同加载方式，加载一个或多个GRIB messages到内存：
　　有两种函数：codes_grib_new_from_file 和 codes_new_from_index。调用后会返回一个唯一的identifier，用于对已加载的GRIB messages进行操纵。

3. 调用codes_get函数对已加载的GRIB messages进行解码; (可以解码需要的数据)

4. 释放已经加载的GRIB messages：
　　codes_release　　

5. 关闭打开的 GRIB 文件.

此外，eccodes还有以下功能：

　　eccodes的主要目的是提供一个高水平的方法，用以从一个加载的GRIB messages对象去提取和计算出更多额外的信息：

· 返回平均，最小，最大，特定的经纬度等的关键字； 

· 计算经纬度和值：codes_grib_get_data； 

· 提取值的子例程：
　　codes_grib_find_nearest 去提取出距离给定地理点最近位置的点的值
　　codes_get_element 从列表中提取值

· 基于索引访问(indexed access)的子例程：这些是随机访问的方法，通常比顺序访问更快。

 

-----------------------------------------------------------------

顺序访问方式：

大致思路：

-> codes_open_file 

-> codes_grib_new_from_file -> codes_get -> codes_release

     …

-> codes_grib_new_from_file -> codes_get-> codes_release 

-> codes_close_file

 

示例代码：

 

 

索引访问方式（通常比顺序访问快）：

注意，eccodes中的index文件(后缀为.idx)与GrADS中后缀为.idx的文件不能通用！

大致思路：

-> codes_index_create(从grib文件创建index) 或 codes_index_read（读取已有index）

-> codes_index_select 选取键值

-> codes_new_from_index -> codes_get -> codes_release

     …

-> codes_new_from_index -> codes_get -> codes_release

-> codes_index_release

 

 示例代码：

-----------------------------------------------------------------

一些ecCodes命令行工具：

参考 https://confluence.ecmwf.int/display/GRIB/GRIB+tools

grib_filter 筛选

grib_filter [options] rules_file grib_file gribfile ...

options

-f　　强制执行

-o　　输出index文件。如果没指定输出文件，那么输出GRIB文件写在filtered.out

-M　　关闭多场支持。关闭在单GRIB message中多要素场的支持。

-V　　版本

-g　　复制GTS头

-G　　兼容GRIBEX模式

-T T|B　　message类型。T -> GTS, B -> BUFR, 输入文件类型根据messgae解释

-7　　当message长度错误时候不报错

-v　　verbose模式

例子

1. grib_filter程序对于输入文件中的所有GRIB message顺序处理，对其中各项均应用rules。输入消息可以使用“write”声明写入输出文件。write语句可以parameterised以便输出发送到多个文件，根据键值定义输出文件的名字。如果我们写一个rules_file包含唯一的声明：

write "../data/split/[centre]_[dataDate]_[dataType]_[levelType].grib[edition]";

将这个rules_file应用到 "../data/tigge_pf_ecmwf.grib2"文件将会在 ../data/split 目录下得到几个文件，包含根据键值分割的字段

> grib_filter rules_file ../data/tigge_pf_ecmwf.grib2
> ls ../data/split
ecmf_20060619_pf_sfc.grib2
ecmf_20060630_pf_sfc.grib2
ecmf_20070122_pf_pl.grib2
ecmf_20070122_pf_pt.grib2
ecmf_20070122_pf_pv.grib2
ecmf_20070122_pf_sfc.grib2

2. 通过明确指示冒号后所需的类型，也可以以不同的格式获取文件名中的键值。

• :s　　整型
• :d　　双精度
• :s　　字符串

以下语句的工作方式与前一个示例略有不同，包括输出文件名中的center和dataType的整数值。

write "../data/split/[centre:i]_[dataDate]_[dataType:i]_[levelType].grib[edition]";

再次运行相同的命令，我们获得了不同的文件列表。

> grib_filter rules_file ../data/tigge_pf_ecmwf.grib2
> ls ../data/split
98_20060619_4_sfc.grib2
98_20060630_4_sfc.grib2
98_20070122_4_pl.grib2
98_20070122_4_pt.grib2
98_20070122_4_pv.grib2
98_20070122_4_sfc.grib2

3. grib_filter语法中允许使用其他语句：

• 以＃开头的注释
• print "string to print also with key values like in the file name"
  
• transient keyname1 = keyname2;
• set keyname = keyvalue;
• 定义（keyname）以检查message中是否定义了关键字
• 缺少（keyname）来检查关键字的值是否设置为MISSING
• 要将键值设置为MISSING，请使用“set key = MISSING;” （见示例）
• if（condition）{block of rules} else {block of rules} 
  条件可以使用==，！=，并使用||和&& 连接单个块条件 
  语句可以是任何有效语句也是另一个嵌套条件
• 您也可以使用“assert（condition）”进行断言。如果condition为false，它将中止filter。
  例如：assert（edition == 1）;

grib_filter规则的一个复杂示例如下是在GRIB版本1文件中更改温度。

# This filter should only be run on GRIB edition 1