- 积分
- 7427
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-9-20
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
最近在琢磨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()
如果在实践过程中,发现错误,求及时留言,谢谢大家。
|
|