最近比較忙好幾天沒有更新，今天就更新一個簡單且實用的數(shù)據(jù)處理方法吧。
在使用python時大家會發(fā)現(xiàn)有時候程序運(yùn)行的異常緩慢。數(shù)據(jù)讀取環(huán)節(jié)是我們不能控制的，那么只有在數(shù)據(jù)處理環(huán)節(jié)，我們才能盡可能的簡化它。讓python計算加快的關(guān)鍵是向量化，就是對數(shù)據(jù)進(jìn)行矩陣運(yùn)算，而不是循環(huán)運(yùn)算。
python語言本身就決定了它不能像C或者Fortran一樣實行動態(tài)運(yùn)算，而在氣象科研中又常常遇到大量的循環(huán)和邏輯判斷，這就讓程序的運(yùn)行異常緩慢，我提供兩個解決辦法，第一個就是前面提到的向量化，舉一個簡單的例子，這個例子很蠢，但是比較容易理解。

#向量運(yùn)算：
a = np.array([[1,2,3],[4,5,6],[7,8,9]])
b = np.array([[3,2,1],[6,5,4],[9,8,7]])
c = a + b

#非向量運(yùn)算：
a = np.array([[1,2,3],[4,5,6],[7,8,9]])
b = np.array([[3,2,1],[6,5,4],[9,8,7]])
c = np.zeros((3,3))
for i in range(3):
    for j in renge(3):
        c[i,j] = a[i,j] + b[i,j]

相信大家都不會用第二種去運(yùn)算，但是我這里只是一個簡單的例子，其實大家翻翻以前寫過的循環(huán)，應(yīng)該都會發(fā)現(xiàn)自己使用過非向量運(yùn)算，最簡單的比如說計算相關(guān)系數(shù)時：

from scipy.stats import pearsonr
r = np.zeros((a.shape[1],a.shape[2]))
p = np.zeros((a.shape[1],a.shape[2]))

for i in range(a.shape[1]):
    for j in range(a.shape[2]):
            r[i,j], p[i,j] = pearsonr(b,a[:,i,j])

其實這里遠(yuǎn)沒有NCL的相關(guān)系數(shù)計算方便，但是我目前了解的python計算相關(guān)系數(shù)的函數(shù)均沒有指定維度的語句，也就是說它們都是求兩個序列的相關(guān)系數(shù)，當(dāng)我們計算一個場與一個序列的相關(guān)系數(shù)時，只能通過循環(huán)實現(xiàn)。
那么，在這種情況下，就可以用到我要講的第二種方法，使用Numba加速！
關(guān)于Numba是什么的問題，大家感興趣可以百度深入了解，通俗地講，就是一個可以提高python遍歷運(yùn)算性能的庫。安裝方法：

conda install numba

使用numba很簡單，就是自己定義一個python函數(shù)，然后添加一個裝飾器來提高性能。還是同樣的計算相關(guān)系數(shù)的程序：

import numba as nb
#為相關(guān)系數(shù)計算函數(shù)塊提供@jit裝飾器
@nb.jit()
def r_caculate(b,a):
    r = np.zeros((a.shape[1],a.shape[2]))
    p = np.zeros((a.shape[1],a.shape[2]))
    for i in range(a.shape[1]):
        for j in range(a.shape[2]):
            r[i,j], p[i,j] = pearsonr(b, a[:,i,j])
    return r,p

f_hadsic =  xr.open_dataset('HadISST_ice.nc')
sic_had = np.array(f_hadsic['sic'])                     
r,p = r_caculate(pc, sic_had)

運(yùn)算結(jié)果是4.39秒，而不適用numba加速時：

f_hadsic =  xr.open_dataset('HadISST_ice.nc')
sic_had = np.array(f_hadsic['sic'])                     
r,p = r_caculate(pc, sic_had)
r = np.zeros((sic_had.shape[1],sic_had.shape[2]))
p = np.zeros((sic_had.shape[1],sic_had.shape[2]))
for i in range(sic_had.shape[1]):
    for j in range(sic_had.shape[2]):
        r[i,j], p[i,j] = pearsonr(pc, sic_had[:,i,j])

計算花費了12.3秒?？梢姴罹嗍志薮蟆?br> 當(dāng)數(shù)據(jù)的循環(huán)和邏輯判斷很大時，使用numba加速可以說會為我們節(jié)省大量的等待時間。
但是numba并不是萬能的，據(jù)我目前使用的結(jié)果，它只支持numpy，scipy的一部分函數(shù)。具體還需要大家在使用過程中逐漸嘗試。
再舉一個例子(這個函數(shù)實現(xiàn)是在一個(時間，緯度，經(jīng)度)數(shù)組中判斷任意時刻任意一個點是否大于周圍八個點的)：

@nb.jit()
def max_in_8surrounds(a):
    b = np.zeros((a.shape[0],a.shape[1],a.shape[2]))
    for t in range(a.shape[0]):
            for i in range(1,a.shape[1]-1):
                for j in range(1,a.shape[2]-1):
                    tmp = np.array([a[t,i-1,j-1],a[t,i+1,j+1],a[t,i-1,j+1],
                                   a[t,i+1,j],a[t,i-1,j],a[t,i,j+1],a[t,i,j-1],a[t,i+1,j-1]])
                    if(a[t,i,j]>np.max(tmp)):
                        b[t,i,j] = 1
    return b
point = max_in_8surrounds(a)

我計算的數(shù)組十分巨大，加入了修飾器后，比未使用numba加速可以說快了不止幾十倍。其實這個程序還涉及了一個小技巧，就是我在判斷一個點是否小于周圍八個點時，我將周圍八個點先同時放進(jìn)了一個數(shù)組，然后用中心點去與數(shù)組的最大值去判斷大小，這樣遠(yuǎn)比進(jìn)行八次邏輯判斷的速度要快的多。

if ((a>a1) & (a>a2) & (a>a2) & (a>a2) & .................):

如果numba仍不能滿足你的需要，那還是去學(xué)C或者Fortran吧QAQ