GrADS格点数据文件解析

目录

本文介绍GrADS格点数据的描述文件的解析和二进制数据文件的读取,关于GrADS格点数据格式,请参照我之前一篇博文《GrADS格点数据格式

ctl文件解析

从ctl文件格式来看,每个独立的ctl语句行都以一个关键字开头。有些关键字只有一行语句,比如dset、title等;有些关键字可能有多行语句,比如zdef、vars等,后续语句行数都可以从第一行语句中得到。
所以解析ctl文件是顺序读取每行,第一个单词为关键字,根据关键字采取相应的操作,比如dset只需读当前行,而vars需要读取后面所有的变量记录加上endvars关键字,返回最后读取的行数。接下来继续从下一行开始解析文件。
框架图后续补充。
代码如下:

lines = f.readlines()
self.ctl_file_lines = [l.strip() for l in lines]
self.cur_no = 0
total_lines = len(lines)
while self.cur_no < total_lines:
    cur_line = lines[self.cur_no]
    first_word = cur_line[0:cur_line.find(' ')]
    if first_word.lower() in self.parser_mapper:
        self.parser_mapper[first_word](self)
    self.cur_no += 1

每个关键字对应一个解析器:

parser_mapper = {
        'ctl_file_name': ctl_file_name_parser,
        'dset': dset_parser,
        'options': options_parser,
        'title': title_parser,
        'undef': undef_parser,
        'xdef': dimension_parser,
        'ydef': dimension_parser,
        'zdef': dimension_parser,
        'tdef': tdef_parser,
        'vars': vars_parser,
    }

下面介绍一些解析器

dset

替换{shell}^{/shell}为绝对路径。

options

查看是否设置big_endian等参数

dimension

适用于xdef和ydef,根据映射类型,采用不同的方式设置层次。
对于linear类型,使用层数、起始值和步进值生成层次列表

levels = [start+step*n for n in range(count)]

对于levels类型,读取后续层次值组成列表。

tdef

对于起始时间,目前只考虑业务系统中没有分钟的形式。

start_date = datetime.datetime.strptime(start_string.lower(), '%Hz%d%b%Y')

对于增量时间,目前不支持月份、年份增量(因为{python}datetime.timedelta{/python}没有显式支持月份和年份)。
最后生成时间列表:

values = [start_date + time_step * i for i in range(count)]

vars

目前只支持v2.0.1及之前版本的格式,不考虑附加码。读取每条记录,分析层次类型(是否使用zdef设置的层次),生成变量列表。
再根据变量列表,加上相应的层次信息,生成数据记录列表,对应数据文件中的每条记录。

二进制文件读取

从格式说明中可以知道,二进制文件分为两种:带fortran长度记录和不带fortran长度记录。通过ctl文件的options中是否设置sequential判断。
fortran长度记录为4字节,首尾各一个。每个数据都是4字节的浮点数。
读取数据时首先要找到记录的开头。因为ctl文件中每个记录的网格相同,所以每条记录占用的空间也相同,所以只需通过记录的索引号即可得到记录的偏移。
带fortran长度记录:

offset = (nx*ny*4+2*4)*a_record_index + 4

不带长度记录

offset = nx*ny*4*a_record_index

之后就是nx*ny个四字节float型数据,先沿x轴存放,再沿y轴存放。
数据的解码与字节顺序有关,通过ctl文件的options中是否设置big_endian或little_endian判断。
使用Python的struct模块解码,大端使用’>f’,小端使用’<f’:

if self.grads_ctl.data_endian == 'big':
    data_format = '>f'
else:
    data_format = '<f'
var_list = [struct.unpack(data_format, data_file.read(4))[0] for i in range(0, y_count*x_count)]

得到数据列表var_list。

参考

详细解析程序可以参见我写的一个小程序:https://github.com/perillaroc/porter