Python气象数据处理与绘图(17):加快循环的运算速度

最近比较忙好几天没有更新,今天就更新一个简单且实用的数据处理方法吧。
在使用python时大家会发现有时候程序运行的异常缓慢。数据读取环节是我们不能控制的,那么只有在数据处理环节,我们才能尽可能的简化它。让python计算加快的关键是向量化,就是对数据进行矩阵运算,而不是循环运算。
python语言本身就决定了它不能像C或者Fortran一样实行动态运算,而在气象科研中又常常遇到大量的循环和逻辑判断,这就让程序的运行异常缓慢,我提供两个解决办法,第一个就是前面提到的向量化,举一个简单的例子,这个例子很蠢,但是比较容易理解。

#向量运算:
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

#非向量运算:
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]

相信大家都不会用第二种去运算,但是我这里只是一个简单的例子,其实大家翻翻以前写过的循环,应该都会发现自己使用过非向量运算,最简单的比如说计算相关系数时:

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])

其实这里远没有NCL的相关系数计算方便,但是我目前了解的python计算相关系数的函数均没有指定维度的语句,也就是说它们都是求两个序列的相关系数,当我们计算一个场与一个序列的相关系数时,只能通过循环实现。
那么,在这种情况下,就可以用到我要讲的第二种方法,使用Numba加速!
关于Numba是什么的问题,大家感兴趣可以百度深入了解,通俗地讲,就是一个可以提高python遍历运算性能的库。安装方法:

conda install numba

使用numba很简单,就是自己定义一个python函数,然后添加一个装饰器来提高性能。还是同样的计算相关系数的程序:

import numba as nb
#为相关系数计算函数块提供@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)

运算结果是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秒。可见差距十分巨大。
当数据的循环和逻辑判断很大时,使用numba加速可以说会为我们节省大量的等待时间。
但是numba并不是万能的,据我目前使用的结果,它只支持numpy,scipy的一部分函数。具体还需要大家在使用过程中逐渐尝试。
再举一个例子(这个函数实现是在一个(时间,纬度,经度)数组中判断任意时刻任意一个点是否大于周围八个点的):

@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)

我计算的数组十分巨大,加入了修饰器后,比未使用numba加速可以说快了不止几十倍。其实这个程序还涉及了一个小技巧,就是我在判断一个点是否小于周围八个点时,我将周围八个点先同时放进了一个数组,然后用中心点去与数组的最大值去判断大小,这样远比进行八次逻辑判断的速度要快的多。

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

如果numba仍不能满足你的需要,那还是去学C或者Fortran吧QAQ

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,711评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,932评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,770评论 0 330
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,799评论 1 271
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,697评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,069评论 1 276
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,535评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,200评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,353评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,290评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,331评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,020评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,610评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,694评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,927评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,330评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,904评论 2 341

推荐阅读更多精彩内容

  • python由于它动态解释性语言的特性,跑起代码来相比java、c++要慢很多,尤其在做科学计算的时候,十亿百亿级...
    python大数据分析阅读 1,617评论 2 27
  • “想什么呢?这么出神。”不知道什么时候,边伯贤已经洗完了澡,来到她的身边坐下,他的头发上还挂着水,带着刚刚沐浴完的...
    霸道毅阅读 184,799评论 1 9
  • “江湖有个小浪子, 尝遍天下金波子。 要说这还差啥子, 独缺貌美小娘子~”妖风坐在一家名为味风的酒肆里,一边喝着杯...
    棠玥阅读 1,188评论 9 10
  • 1.你以为放手可以成全我的幸福,可你不知道,我最大的幸福就是能和你手牵手。 2.快乐,不过是给伤口找一个笑着流泪的...
    duodortlol阅读 171评论 0 0
  • 时光荏苒,岁月如梭,转眼到了和孩子说再见的时候,思绪万千,心里不觉怅然,当看到孩子们那一张张稚嫩的脸,听到那一声声...
    白窑小学尹晓丽阅读 95评论 0 0