GrADS格点数据文件解析
本文介绍 GrADS 格点数据的描述文件的解析和二进制数据文件的读取。
关于 GrADS 格点数据格式的详细介绍，请浏览《GrADS格点数据格式》。
数据描述文件解析
从数据描述文件 (简称 ctl 文件) 格式来看，每个独立的 ctl 语句行都以一个关键字开头。
有些关键字只有一行语句，比如 dset、title 等；有些关键字可能有多行语句，比如 zdef、vars 等，后续语句行数都可以从第一行语句中得到。
所以解析 ctl 文件是顺序读取每行，第一个单词为关键字，根据关键字采取相应的操作，比如 dset 只需读当前行，而 vars 需要读取后面所有的变量记录加上 endvars 关键字，返回最后读取的行数。
接下来继续从下一行开始解析文件。
框架图后续补充。
首先将 ctl 文件内容保存为字符串数组
self.ctl_file_path = Path(ctl_file_path)
with open(ctl_file_path) as f:
lines = f.readlines()
self.ctl_file_lines = [l.strip() for l in lines]
self.cur_no = 0
按顺序解析，读取一行，得到该行对应的解析器，运行解析器，解析器可能会读取更多的行，再从新一行继续解析
total_lines = len(self.ctl_file_lines)
while self.cur_no < total_lines:
cur_line = self.ctl_file_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.cur_no += 1
每个关键字对应一个解析器：
self.parser_mapper = {
'ctl_file_name': self._parse_ctl_file_name,
'dset': self._parse_dset,
'options': self._parse_options,
'title': self._parse_title,
'undef': self._parse_undef,
'xdef': self._parse_dimension,
'ydef': self._parse_dimension,
'zdef': self._parse_dimension,
'tdef': self._parse_tdef,
'vars': self._parse_vars,
}
下面介绍一些解析器
dset
替换 ^ 为绝对路径，例如
dset ^postvar2021080200_024
判断路径是否为模板，例如
dset ^postvar2021092700%f3%n2
options
查看是否设置 big_endian、little_endian、yrev 等参数。
dimension
适用于 xdef 和 ydef，根据映射类型，采用不同的方式设置层次。
对于 linear 类型
xdef 1440 linear 0.0000 0.2500
使用层数、起始值和步进值生成层次列表
levels = [start + step * n for n in range(count)]
对于 levels 类型
zdef 27 levels
1000.000
925.0000
850.0000
700.0000
...
读取后续层次值组成列表。
tdef
对于起始时间，目前只考虑业务系统中没有分钟的形式。
start_date = pd.to_datetime(start_string.lower(), format='%Hz%d%b%Y')
对于时间增量 vvkk，使用 pd.Timedelta 生成时间增量对象 time_step，目前仅支持分钟、小时和天 (即 kk 为 mn，hr 或 dy)。
最后生成时间列表：
values = [start_date + time_step * i for i in range(count)]
vars
目前只支持 v2.0.1 及之前版本的格式，不考虑附加码。
读取每条记录，分析层次类型 (是否使用 zdef 设置的层次)，生成变量列表 self.grads_ctl.vars。
再根据变量列表，加上相应的层次信息，生成数据记录列表 self.grads_ctl.record，对应数据文件中的每条记录。
二进制文件读取
从格式说明中可以知道，二进制文件分为两种：
- 带 fortran 长度记录
- 不带 fortran 长度记录
通过 ctl 文件的 options 中是否设置 sequential 判断。
Fortran 长度记录为 4 字节，首尾各一个。 每个数据都是 4 字节的浮点数。
读取数据时首先要找到记录的开头。 因为 ctl 文件中每个记录的网格相同，所以每条记录占用的空间也相同，所以只需通过记录的索引号即可得到记录的偏移。
带 Fortran 长度记录：
offset = (nx * ny + 2) * 4 * record_index
不带 Fortran 长度记录：
offset = nx * ny * 4 * record_index
之后就是 nx*ny 个四字节 float 型数据，先沿 x 轴存放，再沿 y 轴存放。
数据的解码与字节顺序有关，通过 ctl 文件的 options 中是否设置 big_endian 或 little_endian 判断。
设置文件读取模式，大端使用 >f，小端使用 <f，默认与本地相同 (f)：
if self.grads_ctl.data_endian == 'big':
data_format = '>f'
elif self.grads_ctl.data_endian == 'little':
data_format = '<f'
else:
# print("Data endian is not found. Use local endian to unpack values.")
data_format = 'f'
使用 numpy 的 np.fromrile() 函数从数据文件的 offset 偏移处读取一维数组，并重组为二维数组：
data_file.seek(offset)
values = np.fromfile(
data_file,
dtype=np.dtype(data_format),
count=x_count * y_count
)
values = values.reshape((y_count, x_count))
参考
详细解析程序可以参见我写的一个小程序：
上述解析程序已集成到 nwpc-oper/nwpc-data 项目中。
