IDL 读取grib数据的一些经验之谈

uranuscc

    第一次用IDL尝试读取grib数据，难免会感到各种困惑，比如不知道grib数据内有哪些物理量，该如何选择需要的物理量。以下是我遇到的各种问题及解决经验，有不足之处还请多多指教：

1 IDL是如何读取GRIB数据信息的？  IDL与NCL和grads不同，NCL(ncl_filedump)和Grads读取grib很方便，主要在于他们能生成一个数据描述性文件，通过解读描述性文件能一眼看到数据中包包含的数据量，维度等基本信息。IDL似乎没有这个功能，IDL8.3以后有两个函数grib_list和grib_getdata能获取部分数据的描述信息，但是当针对某些grib数据中的头文件中一些关键字（key）不识别的时候，这两个函数就不是那么好用了，会报错，而且错误信息会让人摸不着头脑，此时只有具体问题具体分析。
  IDL 读取grib数据是根据数据的存储顺利按照数据（record）一条一条读取的，每一条记录(record)是有头文件信息(header)和数据(value)组成，读取一条计入，首先要获得该条记录的句柄（handle，grib_new_from_file（）），然后可通过grib_get_values获得该条记录数据（只有数值，不包括数据描述），该条记录的数据描述存放在header里。

2 获取物理量信息(即每条记录的header信息)

   a.IDL8.2是没有grib_list和grib_getdata函数的，所以只能按照官方例子的方式获取header来查看，数据记录的写入是有规律的，所以只要查看前几条记录就能大致了解数据是如何写入的，然后根据需要取数据就好了。
   b.IDL8.3以后有grib_list和grib_getdata函数，运气好，数据不报错的话，能方便获取整个文件的信息，但是针对一些可能会报错的数据，也只能按照az中方式获取。

3 按照官方例子读header时要注意的一个小问题
   头文件的读取是通过获取每一个关键字(key)及其值来获取该条记录的信息，官方例子中grib_keys_iteartor_new(h,/computed)，这个computed参数设置后可能会滤掉一些我们需要的信息（比如我读的文件中物理量(shortName)关键字就被过滤掉了，找了很久才发现是这个computed参数在搞事.......），选取适当的参数(比如，all)就能获得你想要的信息，然后根据关键字来判断是否要取出该条记录的数据信息。

知道数据是怎么写入的，以及如何选取自己所需要的信息，读出来就不是个问题了。

总结： IDL 读取grib数据不是那么方便，主要还在于读取数据的时候没有描述文件，函数的运用抽象不直观，有点像读二进制数据。但归根结底，读取数据还是数据怎么存入的就怎么读。

初次摸索，抛砖引玉，不足之处还请多多指教。


IDL8.2中的官方例子如下：

[size=13.3333px]FUNCTION GRIB_READ_EX, filename, HEADER=header

[size=13.3333px]  ON_ERROR, 2

[size=13.3333px]  IF(filename eq !null) THEN MESSAGE, 'File is undefined.'

[size=13.3333px]  f = GRIB_OPEN(filename)

[size=13.3333px]  data = PTRARR(GRIB_COUNT(filename))

[size=13.3333px]  if(ARG_PRESENT(header)) THEN header = MAKE_ARRAY(GRIB_COUNT(filename), /OBJ)

[size=13.3333px]  h = GRIB_NEW_FROM_FILE(f)   获取该条记录句柄

[size=13.3333px]  i=0

[size=13.3333px]    WHILE(h NE !NULL) DO BEGIN

[size=13.3333px]    ; Get values array

[size=13.3333px]    values = GRIB_GET_VALUES(h)  获取该条记录数据

[size=13.3333px]    data = PTR_NEW(values)

[size=13.3333px]    ; Get header information if requested

[size=13.3333px]    IF (ARG_PRESENT(header)) THEN BEGIN

[size=13.3333px]      kiter = GRIB_KEYS_ITERATOR_NEW(h, /COMPUTED)  获取关键字的迭代序号（/computed参数要注意）

[size=13.3333px]      header = LIST()

[size=13.3333px]      res = GRIB_KEYS_ITERATOR_NEXT(kiter)    如果下一个关键字存在则...

[size=13.3333px]      WHILE (res EQ 1) DO BEGIN

[size=13.3333px]        key = GRIB_KEYS_ITERATOR_GET_NAME(kiter)

[size=13.3333px]        IF (STRCMP(key, 'values', /FOLD_CASE) EQ 0) THEN BEGIN

[size=13.3333px]          IF (GRIB_GET_SIZE(h, key) GT 1) THEN $

[size=13.3333px]            val = GRIB_GET_ARRAY(h, key)ELSE val = GRIB_GET(h, key)  关键字变量是一个数组则grib_get_array，否者grib_get获取关键字的值

[size=13.3333px]          IF (STRCMP(key, '7777', /FOLD_CASE) EQ 1) THEN key = 'end_section' 如果是最后一个关键字，则end_section

[size=13.3333px]          key_value = CREATE_STRUCT(key, val)

[size=13.3333px]          header.add, key_value

[size=13.3333px]        ENDIF

[size=13.3333px]        res = GRIB_KEYS_ITERATOR_NEXT(kiter)

[size=13.3333px]      ENDWHILE

[size=13.3333px]      GRIB_KEYS_ITERATOR_DELETE, kiter

[size=13.3333px]    ENDIF

[size=13.3333px]    GRIB_RELEASE, h  释放句柄

[size=13.3333px]    h = GRIB_NEW_FROM_FILE(f) 获取下一个句柄

[size=13.3333px]    i++

[size=13.3333px]  ENDWHILE

[size=13.3333px]  GRIB_CLOSE, f

[size=13.3333px]  RETURN, data

[size=13.3333px]END

mrxiaoan

这段代码如果用来读fnl资料，可能是会有问题的，问题主要出在：Parameter information、grib 2 Section 5 DATA REPRESENTATION SECTION、grib 2 Section 6 BIT-MAP SECTION、grib 2 Section 7 data这4个字段上，这4个字段在GRIB_GET(h, key)，会出现错误，提示“Invalid type”，如果忽略这几个字段，数据也就没了...

uranuscc

mrxiaoan 发表于 2018-11-27 17:49
这段代码如果用来读fnl资料，可能是会有问题的，问题主要出在：Parameter information、grib 2 Section 5 D ...

所以说还是不要用idl来读grib了，不是很方便

gege

请问楼主解决invalid type 这个问题吗？