气候分析脚本：计算距平

虫儿飞

最近在琢磨numpy，总觉得应该有一个小目标才有动力，所以就打算实现一些常用的气候分析方法。气候分析，距平算是比较重要的概念之一。
先上脚本：
import numpy as np
import numpy.ma as ma

def SeparateIdate(idate):

    if isinstance(idate,int):
        sdate = str(idate)
    else:
        sdate = idate
    if len(sdate) == 4:
        year = int(sdate)
        month = 0
        day = 0
        n = 0
    elif len(sdate) == 6:
        year = int(sdate[0:4])
        month = int(sdate[4:6])
        day = 0
        n = month-1
    elif len(sdate) == 8:
        year = int(sdate[0:4])
        month = int(sdate[4:6])
        day = int(sdate[6:8])
        days=[31,28,31,30,31,30,31,31,30,31,30,31]
        if (year%4 == 0 and year%100 != 0) or year%400 == 0:
            days=[31,29,31,30,31,30,31,31,30,31,30,31]
        n = 0
        for imon in range(1,month):
            n = n+days[imon-1]
        n=n+day-1
    else:
        print 'SeparateIdate-Warning(1): The length of idate should be four(4) or six(6) or eight(8). The current length is',len(sdate),'.'
        exit()
    #print year,month,day,n
    return year,month,day,n

def CalculateAnomaly(data,InitIdate,ClimatePeriod='1981 2010',Type='M'):
    ###data: data[time,lat,lon]
    ###Type: D day      InitIdate: '19810101'
    ###      M month    InitIdate: '198101'
    ###      S season   InitIdate: '198101' YYYY01 spring YYYY02 summer YYYY03 autumn YYYY04 winter
    ###      Y year     InitIdate: '1981'

    ### data check
    if len(data.shape) != 3:
        print 'CalculateAnomaly-Warning(1): The dimension of array <data> should be three(3). The current dimension is',len(data.shape),'.'
        exit()
    else:
        #print type(data),data.shape
        #print data[0,...]
        n = data.shape[0]
        lats = data.shape[1]
        lons = data.shape[2]
    ### date check
    if not isinstance(InitIdate,(str,int)):
        print 'CalculateAnomaly-Warning(2): The format of idate should be string(str) or integer(int).'
        exit()
    if ( not isinstance(ClimatePeriod,str) ) or len(ClimatePeriod) != 9:
        print 'CalculateAnomaly-Warning(3): The format of period should be string(str), and its length should be nine(9).'
        exit()
    else:
        begin_year = int(ClimatePeriod.split()[0])
        end_year = int(ClimatePeriod.split()[1])
    init_year,init_month,init_day,init_n = SeparateIdate(InitIdate)
    if Type == 'D':
        interval = 365   
    elif Type == 'M':
        interval = 12
    elif Type == 'S':
        interval = 4
    elif Type == 'Y':
        interval = 1
    if init_year > begin_year or ((end_year-init_year+1)*interval-init_n) > n:
        print 'CalculateAnomaly-Warning(4): The period of <data> do not contain the climate state period.'
        exit()
    for i in range(1,interval+1):
        t1 = (begin_year-init_year)*interval-init_n+i-1
        t2 = (end_year-init_year)*interval-init_n+i
        x = data[t1:t2:interval,:,:]
        t1 = i-(init_n+1)
        if t1 < 0:
            t1=t1+interval
        t2 = n
        #print t1,t2
        #print range(t1,t2,interval)
        #data[t1:t2:interval,:,:] = ma.subtract(data[t1:t2:interval,:,:],ma.mean(x,axis=0))
        data[t1:t2:interval,:,:] = data[t1:t2:interval,:,:].__sub__(ma.mean(x,axis=0))
    #print data[1824,::5,::5]
    return data


脚本保存的名字叫做ClimateAnalysis.py。
data为numpy的Masked array，需要事先进行mask处理，InitIdate格式形如为‘1854’、‘185401’、‘18540101’，或者已经相应的1854、185401、18540101，ClimatePeriod为所选取的气候态，并不强制为三十年，形如'1981 2010'，必须是长度为9的字符串，中间为空格，依赖此进行split,Type为年、季节、月、天尺度的选取参数，参看源代码说明。
这个脚本可以实现年、季节、月、天（事先进行去闰年）尺度的距平计算，调用格式，以月尺度的海温数据（'sst.mnmean.nc'）为例，如下所示：
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from ClimateAnalysis import *

infile = 'sst.mnmean.nc'
fopen = Dataset(infile)
lons = fopen.variables['lon'][:]
lats = fopen.variables['lat'][:]
sst = fopen.variables['sst'][:]
missing_value = fopen.variables['sst']._FillValue
sst=CalculateAnomaly(sst,'185401')
plt.pcolormesh(lons, lats, sst[1826,:,:])
plt.show()


如果在实践过程中，发现错误，求及时留言，谢谢大家。

王先生

你好，例子中sst=CalculateAnomaly(sst,'185401') 这句是什么意思啊？我读的时候没看懂。

weinihou

感谢分享~~~


看到楼主里面的退出直接用的 exit()  ，建议换成 sys.exit() ，避免可能出现的文件运行方式下的问题

wsjqxjy

{:eb502:}{:eb502:}

虫儿飞

weinihou 发表于 2017-8-6 21:54
感谢分享~~~

好哒，谢谢~

milletppp

大好人啦。

badboy1204

这个很好用，以后计算降水或者是植被变化，都很有参考价值，谢谢楼主分享

lgz

以后计算降水或者是植被变化

王先生

这个帖子就厉害了啊。膜拜一下。顺便问一下脚本print后面没有（）提示错误怎么办呢？加还是不管啊

虫儿飞

王先生 发表于 2019-4-21 23:00
这个帖子就厉害了啊。膜拜一下。顺便问一下脚本print后面没有（）提示错误怎么办呢？加还是不管啊

我的脚本基于Python2写的，你是用的Python3吧？如果是就需要加