python精彩編程200例，python 3.6.3自帶的編程調試環境包括了_序章：資料預處理(python3.6 可用fortran unformatted

首先我只是一個接觸Python約半年的菜鳥，開這一個專欄的目的主要是記錄自己所學，以及實踐的一些有用的東西，順便分享一些自己寫的公用代碼段以方便具有同樣想法的朋友。

既然是序章我就多寫一些吧，我本人對人工智能在氣象方面的應用特別感興趣，有一樣想法的你歡迎email到(wenqs@grmc.gov.cn)

為什么是Python?

這個問題其實被問道了很多次了，相對其他行業氣象尤其是天氣業務類里能真充分發揮Python強大處理能力的地方其實比較有限，并且本來氣象學就是一個偏重理解的學科，所以許多大佬們對這種新型的面向對象語言和其他編程語言的區別也不太在意。

python精彩編程200例？但這里既然是我的場地，我要說：python的好處我這菜雞可能體會得不深，不過有三點需要明確：

想要使用最火熱的人工智能算法你繞不過python，應為很多很多的例子就是直接用python書寫的。

我并不是專業編程科班出身，接觸的編程語言并不多，但是在我接觸到的里面，python的語法結構是最清晰，并且本人認為是最接近人類思維邏輯及自然語言的。這導致很多代碼段其實就跟作文有一點點相似了。

Python有非常強大的社區支持，無論你程序編寫時，還是查找錯誤時，有用的幫助信息無處不在，同時網上許多的可用模塊都統樣高效，也保持著較為一致的語法特點。

為什么需要資料預處理呢？

python的開發環境有哪些、——一句話，巧婦難為無米之炊啊！

這一點其實比較尷尬，請注意我這里說的“預處理”還沒有到許多機器學習教程提到的“Data clean proccess”，僅僅只是將數據讀入python。。。。

由于一些歷史原因，國內天氣預報業務用的數值模式預報產品一般采用兩種格式：

標準的NetCDF格式格點資料。(這種資料網上到處是讀取模塊，這里就不贅述了)

Fortran的“無格式二進制順序存取”文件(fortran unformatted sequencial data)。這種文件在不同的操作系統中還細分為big-endian與little-endian版本。而且在存放高位數組集合時，將他們統一的看成很多個二維數組的疊加，然后存放每一個二維數組時會在數組的一頭一尾添加特定的占位符，然后再在更高維度重讀這種操作，所以直接用python二進制文件讀取模塊會因為錯位問題根本讀取不到想要的信息。endian問題和占位符問題也是網上很多文件讀取教程根本無法正常讀取氣象模式預報數據的原因。(本章只針對這種格式)

如何安裝python環境，很不幸，我工作的生產環境采用了第二種，這種格式由于太過時，網上python對這種格式的支持并不好，一般的教程頂多叫你用numpy.fromfile()等等的方法定制特定數據類型再嘗試讀取，但是講得都不夠深入。另一種做飯是通過一些強大的數據格式轉換軟件如：CDO等等，將數據轉化為NetCDF再進行讀取，可是這樣做即不效率又需要雙倍的存儲空間(我也曾經嘗試過這種做法，實在是不好用)。

于是就誕生了自己書寫可用的讀取模塊的沖動

這里首先說明，這個模塊的設計思路來自一篇網頁，但是作者停止了更新，于是我按照這個思路成功的重寫了適合于grapes模式輸出結果的讀取模塊CTLReader,完整的測試數據及和代碼在github中，歡迎大家一起開發完善。

為了能夠讓大家看得懂代碼，我在代碼中進行了詳細的中文注釋，不需要的可以刪除。

下面通過截圖來說明幾個有意思的代碼段

pycharm和python區別、圖片.png

這一塊為大家都會用的import

圖片.png

這這一塊我們進行了big/little endian的轉換，一次性搞定以后就不需要類似>f4等等的類型說明符了。

圖片.png

為了提高python代碼運行速度和進行？代碼中有很多類似這樣的正則表達式定義，網上有很多詳細的說明，使用它們可以很方便的處理提取文本中的有用信息，具體到這個表達式的意思是：

匹配這樣一段字符串“它以任意字母或數字開頭重復一次或多次，后面接著一個或多個空格，再后面接著一個或多個數字，再后面接著一個或多個空格，再后面接著一個或多個數字，再后面接著一個或多個空格，最后面是任意字符串的組合”

具體到我們的test.ctl文件它能匹配到：

圖片.png

紅色圈圈表示正則式里的()中的內容

我們編寫的python代碼在運行過程中。利用正則表達式和python中類型的定義我們愉快的完成了變量的分類

接下來，這一段里：

圖片.png

看起來比較復雜，其實這里體現了python強大的數組管理功能，近乎完美的將這種疊成放置且每一個二維數組頭尾各有一段多余占位符的數據處理了。首先[i:i+size]指定了類似一個“記錄長度的范圍”形成一個具有reshape方法的“子集”。然后，該方法的-1參數表示將這個子集的所有數據按照原本的大小進行遍歷，然后在利用計算出的二維平面大小去迭代這段數據，相當于不用指定層數python自動把一個高維數組(這里是三維或思維)疊成了一疊由二維數組構造的“千層餅”。[:,int(place_hold/2):-int(place_hold/2)]剔除掉了“每一平面層”不需要的一頭一尾，這樣得到的子集再按照應該有的變量維度進行reshape。(真是非常方便呢！)

圖片.png

python編譯器？上面圖片的這一塊按照本人的工作環境進行了特意的布置，一般的思路是既然ctl文件里有時間信息就應該直接獲取之，但是在數值模式的產品中往往存在多個時間節點的數據，這些數據又是單獨以不同文件的形式存放的，這樣如果需要讀取特定文件要不就是遍歷所有的預報時次，要不就是生成新的ctl文件，這樣都不效率，所有這里我直接將從文件名獲取的時間信息寫入變量屬性中，這一段可以根據自己的需要相應修改。

以下放上模塊的主要代碼，不想去github下載的同學直接復制就可以用了：

import pandas as pd #

import numpy as np #

#，擅長處理各種各樣的數據類型，還能以object形式組建數組。

Python開發環境、import dateparser #

# python中的datetime類型對象，實現進一步處理。

import datetime #

import re #

import os #

以下哪項python能正常啟動。#以上是本腳本主體部分需要的功能模塊。

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~~