- 积分
- 1633
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-1-5
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
这是原值相关,判断缺省值如何写:
#打开nino3.4数据文件
table = readtable('f:/nino34.csv', delimiter=',', format='%s%f')
yy = table['yyJJA']
nino34 = table['nino34']
#print nino34
#打开降水数据文件
f = addfile('f:/precip.mon.mean.2.5x2.5.nc')
#prep = f['precip'][0::12,'0:60','60:150']
var = f['precip'][0,'0:60','60:150']
prepJJA = []
for i in range(1981,2016):
tidx = (i-1948)*12
t = f.gettime(tidx)
data = f['precip'][tidx+5,'0:60','60:150']+\
f['precip'][tidx+6,'0:60','60:150']+\
f['precip'][tidx+7,'0:60','60:150']
prepJJA.append(data)
years = arange(1981,2016)
nx = var.dimlen(1)
ny = var.dimlen(0)
cor = zeros([ny,nx])
cor = dim_array(cor,var.dims)
for i in arange(0,nx):
for j in arange(0,ny):
x = []
y = []
for t in syears:
tidx = t -1981
if(prepJJA[tidx][j,i]>-9.96921E36):
x.append(nino34[t-1981])
y.append(prepJJA[tidx][j,i])
#scatter(xx, yy, s=2, color='k', label='PM')
if(len(x)>1):
slope, intercept, r, p, std,number = linregress(x, y)
cor[j,i] = r
else:
cor[j,i] = nan
#print cor
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
clevs =[-0.53,-0.43,-0.33,-0.28,0,0.28,0.33,0.43,0.53]
cols=[(86,49,5),(138,80,0),(186,145,42),(216,179,107),(245,236,212),\
(217,238,235),(117,197,186),(46,145,137),(0,99,91),(0,59,47)]
#layer = contourfm(cor,20)
layer = contourfm(cor,clevs,colors=cols)
colorbar(layer,orientation='horizontal')
xlim(60,150)
ylim(0,60)
title('correlation:nino34(JJA)-prep(JJA) 1981-2015')
|
-
|