首先我只是一个接触Python约半年的菜鸟，开这一个专栏的目的主要是记录自己所学，以及实践的一些有用的东西，顺便分享一些自己写的公用代码段以方便具有同样想法的朋友。

既然是序章我就多写一些吧，我本人对人工智能在气象方面的应用特别感兴趣，有一样想法的你欢迎email到(wenqs@grmc.gov.cn)

为什么是Python?

这个问题其实被问道了很多次了，相对其他行业气象尤其是天气业务类里能真充分发挥Python强大处理能力的地方其实比较有限，并且本来气象学就是一个偏重理解的学科，所以许多大佬们对这种新型的面向对象语言和其他编程语言的区别也不太在意。

但这里既然是我的场地，我要说：python的好处我这菜鸡可能体会得不深，不过有三点需要明确：

想要使用最火热的人工智能算法你绕不过python，应为很多很多的例子就是直接用python书写的。

我并不是专业编程科班出身，接触的编程语言并不多，但是在我接触到的里面，python的语法结构是最清晰，并且本人认为是最接近人类思维逻辑及自然语言的。这导致很多代码段其实就跟作文有一点点相似了。

Python有非常强大的社区支持，无论你程序编写时，还是查找错误时，有用的帮助信息无处不在，同时网上许多的可用模块都统样高效，也保持着较为一致的语法特点。

为什么需要资料预处理呢？

——一句话，巧妇难为无米之炊啊！

这一点其实比较尴尬，请注意我这里说的“预处理”还没有到许多机器学习教程提到的“Data clean proccess”，仅仅只是将数据读入python。。。。

由于一些历史原因，国内天气预报业务用的数值模式预报产品一般采用两种格式：

标准的NetCDF格式格点资料。(这种资料网上到处是读取模块，这里就不赘述了)

Fortran的“无格式二进制顺序存取”文件(fortran unformatted sequencial data)。这种文件在不同的操作系统中还细分为big-endian与little-endian版本。而且在存放高位数组集合时，将他们统一的看成很多个二维数组的叠加，然后存放每一个二维数组时会在数组的一头一尾添加特定的占位符，然后再在更高维度重读这种操作，所以直接用python二进制文件读取模块会因为错位问题根本读取不到想要的信息。endian问题和占位符问题也是网上很多文件读取教程根本无法正常读取气象模式预报数据的原因。(本章只针对这种格式)

很不幸，我工作的生产环境采用了第二种，这种格式由于太过时，网上python对这种格式的支持并不好，一般的教程顶多叫你用numpy.fromfile()等等的方法定制特定数据类型再尝试读取，但是讲得都不够深入。另一种做饭是通过一些强大的数据格式转换软件如：CDO等等，将数据转化为NetCDF再进行读取，可是这样做即不效率又需要双倍的存储空间(我也曾经尝试过这种做法，实在是不好用)。

于是就诞生了自己书写可用的读取模块的冲动

这里首先说明，这个模块的设计思路来自一篇网页，但是作者停止了更新，于是我按照这个思路成功的重写了适合于grapes模式输出结果的读取模块CTLReader,完整的测试数据及和代码在github中，欢迎大家一起开发完善。

为了能够让大家看得懂代码，我在代码中进行了详细的中文注释，不需要的可以删除。

下面通过截图来说明几个有意思的代码段

图片.png

这一块为大家都会用的import

图片.png

这这一块我们进行了big/little endian的转换，一次性搞定以后就不需要类似>f4等等的类型说明符了。

图片.png

代码中有很多类似这样的正则表达式定义，网上有很多详细的说明，使用它们可以很方便的处理提取文本中的有用信息，具体到这个表达式的意思是：

匹配这样一段字符串“它以任意字母或数字开头重复一次或多次，后面接着一个或多个空格，再后面接着一个或多个数字，再后面接着一个或多个空格，再后面接着一个或多个数字，再后面接着一个或多个空格，最后面是任意字符串的组合”

具体到我们的test.ctl文件它能匹配到：

图片.png

红色圈圈表示正则式里的()中的内容

利用正则表达式和python中类型的定义我们愉快的完成了变量的分类

接下来，这一段里：

图片.png

看起来比较复杂，其实这里体现了python强大的数组管理功能，近乎完美的将这种叠成放置且每一个二维数组头尾各有一段多余占位符的数据处理了。首先[i:i+size]指定了类似一个“记录长度的范围”形成一个具有reshape方法的“子集”。然后，该方法的-1参数表示将这个子集的所有数据按照原本的大小进行遍历，然后在利用计算出的二维平面大小去迭代这段数据，相当于不用指定层数python自动把一个高维数组(这里是三维或思维)叠成了一叠由二维数组构造的“千层饼”。[:,int(place_hold/2):-int(place_hold/2)]剔除掉了“每一平面层”不需要的一头一尾，这样得到的子集再按照应该有的变量维度进行reshape。(真是非常方便呢！)

图片.png

上面图片的这一块按照本人的工作环境进行了特意的布置，一般的思路是既然ctl文件里有时间信息就应该直接获取之，但是在数值模式的产品中往往存在多个时间节点的数据，这些数据又是单独以不同文件的形式存放的，这样如果需要读取特定文件要不就是遍历所有的预报时次，要不就是生成新的ctl文件，这样都不效率，所有这里我直接将从文件名获取的时间信息写入变量属性中，这一段可以根据自己的需要相应修改。

以下放上模块的主要代码，不想去github下载的同学直接复制就可以用了：

import pandas as pd #

import numpy as np #

#，擅长处理各种各样的数据类型，还能以object形式组建数组。

import dateparser #

# python中的datetime类型对象，实现进一步处理。

import datetime #

import re #

import os #

#以上是本脚本主体部分需要的功能模块。

from mpl_toolkits.basemap import Basemap

import matplotlib.pyplot as plt

#以上是出图以测试数据需要用到的模块。

NUMBER = '[-+]?[0-9]*\.?[0-9]+(?:[eE][-+]?[0-9]+)?' #识别一定长度或科学计数法的范例，因经常用到就单独写了

class CTLReader(object):

def __init__(self,filepath,filename,place_hold=2):

self.dimensions = {}

self.variables = {}

self.ctlpath = filepath

self.filename = filename

#将ctl文件信息读入一个巨大的字符串中便于之后处理

with open(self.ctlpath,'r') as f:

self.ctl = f.read()

self._read_data() #读取二进制文件数据

self._read_dimensions() #获取ctl中的维度信息

self._read_vars(place_hold) #将二进制文件数据规整为变量明命名的数组

def _read_data(self):

self.undef = eval(re.search('undef (%s)' % NUMBER , self.ctl).group(1)) #获取CTL文件中的缺省值信息

big_endian = bool(re.search('options.*big_endian',self.ctl,flags=re.I)) #探测数据是否是big_endian

data = np.fromfile(self.filename,'f4') #以4bytes的浮点形式(单精度)读取二进制文件

if big_endian:

data = data.byteswap() #统一将big_endian数据进行位调换

self.data = np.ma.masked_values(data,self.undef) #建立带有默认缺省值的numpy数组并添加到类的自身属性中

def _read_dimensions(self):

if 'xdef' in self.ctl: #探测是否存在xdef关键字

p = re.compile("%s\s+(\d+)\s+linear\s+(%s)\s+(%s)" % ('xdef',NUMBER,NUMBER)) #创建正则维度信息范式

m = p.search(self.ctl)

self.variables['longitude'] = np.linspace(float(m.group(2)),

float(m.group(2))+float(m.group(3))*(int(m.group(1))-1),

int(m.group(1)))

self.dimensions['longitude'] = int(m.group(1))

if 'ydef' in self.ctl: #探测是否存在ydef关键字

p = re.compile("%s\s+(\d+)\s+linear\s+(%s)\s+(%s)" % ('ydef',NUMBER,NUMBER)) #创建正则维度信息范式

m = p.search(self.ctl)

self.variables['latitude'] = np.linspace(float(m.group(2)),

float(m.group(2))+float(m.group(3))*(int(m.group(1))-1),

int(m.group(1)))

self.dimensions['latitude'] = int(m.group(1))

if 'zdef' in self.ctl: #探测是否存在zdef关键字

self.variables['levels'] = Variable('levels',self._parse_dimension('zdef')) #创建“层数”信息变量

self.dimensions['levels'] = len(self.variables['levels'])

if 'grapes' in self.ctl: #探测是否存在grapes关键字

self.variables['time'] = Variable('time',self._parse_dimension('time')) #创建“时间”信息变量

#目前只需要处理“单片”时次的数据

self.dimensions['time'] = 1

def _read_vars(self,place_hold):

read = False #是否识别为目标变量的开关

for line in self.ctl.split('\n'):

if line.startswith('endvars'): #探测目标变量组结束符号

read = False

if read:

p = re.compile('(\w+)\s+(\d+)\s+(\d+)\s+(.*)') #目标变量行的正则范式

m = p.match(line)

name = m.group(1)

var = self.variables[name] = Variable(name) #生成特定的变量类并在本段方法中以"var"的别名进行描述

levels = int(m.group(2))

SPACE = self.dimensions['latitude']*self.dimensions['longitude']

if levels > 0:

var.dimensions = ('time','levels','latitude','longitude') #当变量为四维数组时变量的维度信息

size = self.dimensions['time']*self.dimensions['levels']*(SPACE+place_hold)

else:

var.dimensions = ('time','latitude','longitude') #当变量为三维数组时变量的维度信息

size = self.dimensions['time']*(SPACE+place_hold)

var.shape = tuple(self.dimensions[dim] for dim in var.dimensions) #根据不同的维度信息创建维度宽度提示元组

var.data = self.data[i:i+size].reshape(-1,SPACE+place_hold)[:,

int(place_hold/2):

-int(place_hold/2)].reshape(var.shape)

#以上操作较复杂，主要就是重构数据，去掉头尾的占位符，再次按照维度重构数据

i += size #相当与跳过一定长度的二进制数据字段

units = int(m.group(3)) #单位信息，由于目前阶段处理数据不复杂，暂时不需要添加

if units != 0: #变量的量级转化开关(这种功能交给pandas等模拟自动做吧^_^)

raise NotImplementedError('for now only 0 units will be implemented!')

var.attributes = {

'long_name': m.group(4).strip(),

'units': 'not needed right now'

}

#以上是变量的描述信息，及单位的存放属性

if line.startswith('var'): #探测目标变量组开始符号

i = 0

read = True

def _parse_dimension(self,dim): #用于检索CTL信息中维度相关信息的方法

p = re.compile("%s\s+(\d+)\s+levels([\s\S]+)tdef" % (dim)) #获取层数的具体信息的正则范式

m = p.search(self.ctl)

if m:

return np.fromstring(m.group(2),sep='\n') #以换行符分离目标信息，并生成numpy数组

#time info read from file name

if dim == 'time': #对时间信息的定制处理

filetime = os.path.basename(self.filename)

p = re.compile('mars3km(\d{8})(\d+)')

m = p.search(filetime)

date = m.group(1)

initime = dateparser.parse("20%s %s %s-%s:00:00" % (date[:2],date[2:4],date[4:6],date[6:8]))

endtime = initime + datetime.timedelta(hours=int(m.group(2)))

p = re.compile('\s+\d+\s+linear\s+[:\w]+\s+(\d+)(\w{2})')

m = p.search(self.ctl)

if m:

if m.group(2) == 'mn':

increment = datetime.timedelta(minutes=int(m.group(1)))

else:

increment = datetime.timedelta(hours=int(m.group(1)))

return np.array([initime,endtime,increment])

class Variable(object): #变量类定义

def __init__(self,name,data=None): #创世纪

self.name = name #python说：“要有名字“！于是有了变量

self.data = data #python说：”要有数据“！于是有了变量

def __getitem__(self,index): #python说：”要有方法“！于是有了变量

return self.data[index]

def __getattr__(self,key):

return self.attributes[key]

def __len__(self):

return len(self.data)

最后，我这么做只是希望能方便的将模式数据读取到Python 中方便接下来的人工智能应用，如果下面还有合适分享的公用代码我还是会分享到简书和github上的。

最后的最后，祝福大家狗年吉祥如意！工作这半年确实学习到了不少好东西，希望狗年能尽快将方法应用到实际生产生活中。

顺便帮同学打个广告，我码字这么轻松就能写这么多主要是多亏了有“航天枸杞”保驾护航~~~_

图片.png

图片.png

忘记了最重要的数据读取测试结果了>.<

补充如下：

图片.png

可以看到读取数据画出来的反射率结果完全一致，说明读取数据是成功的~~oh,year~~