请各位帮忙看看这个正距平概率合成有没有什么问题

qxsw2016

按照国家气候中心http://cmdp.ncc-cma.net/pred/cn_enso.php?product=cn_enso_lanina&season=decaying_JJA&area=asia&elem=prep#Lanina做个合成图。脚本代码如下
f = addfile('f:/precip.mon.mean.2.5x2.5.nc')
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer,edgecolor=(0,0,255))
prepJJA = []
for i in range(1981,2018):
    tidx = (i-1948)*12
    t = f.gettime(tidx)
    data = f['precip'][tidx+5,:,:]+f['precip'][tidx+6,:,:]+f['precip'][tidx+7,:,:]   
    prepJJA.append(data)
prepmJJA = sum(prepJJA[0:30])/30.0

for t in range(0,37):
    prepJJA[t] = (prepJJA[t]-prepmJJA)*100/prepmJJA
    prepJJA[t] = prepJJA[t] /abs(prepJJA[t])
    prepJJA[t][prepJJA[t]<0] = 0
years = [1985,1989,1996,2000,2001,2008,2011,2012]
pos_ana_probs = []
for year in years:
    yid = year - 1981
    #print yid   
    pos_ana_probs.append(prepJJA[yid])
#print pos_ana_probs[2]
pos_ana_prob = sum(pos_ana_probs[0:])*100/len(years)

clevs =[0,10,20,30,40,50,60,70,80,90,100]
cols=[(182,106,40),(205,133,63),(225,165,100),(245,205,132),(245,224,158),\
    (255,245,186),(205,255,205),(153,240,178),(83,189,159),(110,170,200),\
    (5,112,176),(2,56,88)]
#layer = contourfm(com_ana_pre,20)  
layer = imshow(pos_ana_prob,clevs, colors=cols)
colorbar(layer,orientation='horizontal')
xlim(60,150)
ylim(0,60)
title('Positive Anomaly Probability(%)')

MeteoInfo

别的没仔细看，不过imshow需要改为imshowm

qxsw2016

imshowm报错

MeteoInfo

qxsw2016 发表于 2018-4-18 21:18
imshowm报错

你可以加入MeteoInfo QQ群方便交流。

qxsw2016

国家气候中心合成图。降水距平百分率脚本如下：
f = addfile('f:/precip.mon.mean.2.5x2.5.nc')
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer,edgecolor=(0,0,255))
prepJJA = []
for i in range(1981,2018):
    tidx = (i-1948)*12
    t = f.gettime(tidx)
    data = f['precip'][tidx+5,:,:]+f['precip'][tidx+6,:,:]+f['precip'][tidx+7,:,:]   
    prepJJA.append(data)
prepmJJA = sum(prepJJA[0:30])/30.0
for t in range(0,37):
    prepJJA[t] = (prepJJA[t]-prepmJJA)*100/prepmJJA
years = [1985,1989,1996,2000,2001,2008,2011,2012]
com_ana_pres = []
for year in years:
    yid = year - 1981
    #print yid   
    com_ana_pres.append(prepJJA[yid])
com_ana_pre = sum(com_ana_pres[0:])/len(years)
clevs =[-50,-40,-30,-20,-10,0,10,20,30,40,50]
cols=[(182,106,40),(205,133,63),(225,165,100),(245,205,132),(245,224,158),\
    (255,245,186),(205,255,205),(153,240,178),(83,189,159),(110,170,200),\
    (5,112,176),(2,56,88)]
#layer = contourfm(com_ana_pre,20)  
layer = contourfm(com_ana_pre,clevs, colors=cols)
colorbar(layer,orientation='horizontal')
xlim(60,150)
ylim(0,60)
title('Anomaly Percentage(%)')

MeteoInfo

qxsw2016 发表于 2018-4-18 22:39
国家气候中心合成图。降水距平百分率脚本如下：
f = addfile('f:/precip.mon.mean.2.5x2.5.nc')
axesm()
...

嗯，和气候中心的图基本一致。

qxsw2016

在老师的指导下，进行了显著性检验，基本一致，但略有差异。脚本如下：
f = addfile('E:/MeteoInfo/Dataset/precip.mon.mean.2.5x2.5.nc')
print f
startmon = 6
monn= 3

var = 'precip'
yr = '0:60'
xr = '60:150'
vardata = f[var[0,yr,xr

y = vardata.dimvalue(0)
x = vardata.dimvalue(1)
ny = vardata.dimlen(0)
nx = vardata.dimlen(1)


prepJJA = [
for i in range(1981,2018):
    tidx = (i-1948)*12
    starttidx = tidx+startmon-1   
    t = f.gettime(tidx)
   
    data = sum(f[var[starttidx:starttidx+monn,yr,xr,axis=0)   
    #data = f[var][starttidx,yr,xr]*30+f[var][starttidx+1,yr,xr]*31+f[var][starttidx+2,yr,xr]*31
    prepJJA.append(data)
   
prepmJJA = sum(prepJJA[0:30)/30.0  #climate 1981-2010

for t in range(0,37):
    prepJJA[t = (prepJJA[t-prepmJJA)*100/prepmJJA
#print prepJJA[0][0,:]
years = [1985,1989,1996,2000,2001,2008,2011,2012
com_ana_preps = [
for year in years:
    yid = year - 1981
    #print yid   
    com_ana_preps.append(prepJJA[yid)
com_ana_prep = sum(com_ana_preps[0:)/len(years)

tt = zeros([ny,nx)
pval = zeros([ny,nx)
an = len(prepJJA)
bn = len(com_ana_preps)
cn = an-bn
#print cn
for i in range(ny):
    for j in range(nx):
        arra = zeros(an)
        kk =0
        for k in range(an):
            if(not (k+1981) in years):
              arra[kk = prepJJA[k[i,j
              kk=kk+1
              #print k+1981
        #print kk
        arra = zeros(an)
        arrb = zeros(bn)
        for k in range(bn):
            arrb[k = com_ana_preps[k[i,j
        #print arra
        #print arrb
        t,p = stats.ttest_ind(arra,arrb)
        #t,p = stats.ttest_1samp(arrb,0)
        tt[i,j = t
        pval[i,j = p
        #print p

tt[pval>0.05 = nan
#print tt

axesm()
mlayer = shaperead('E:\MeteoInfo\MeteoInfo\Map\country.shp')
geoshow(mlayer,edgecolor=(0,0,255))

clevs =[-50,-40,-30,-20,-10,0,10,20,30,40,50
cols=[(182,106,40),(205,133,63),(225,165,100),(245,205,132),(245,224,158),\
    (255,245,186),(205,255,205),(153,240,178),(83,189,159),(110,170,200),\
    (5,112,176),(2,56,88)
#layer = contourfm(com_ana_pre,20)  
layer = contourfm(com_ana_pre,clevs, colors=cols)
#colorbar(layer,orientation='horizontal')
lon,lat = meshgrid(x,y)
lon[tt==nan = nan
lat[tt==nan = nan
layer1 = scatterm(lon,lat,tt,color='k',size = 3.0,marker = '+')
#colorbar(layer1)
xlim(60,150)
ylim(0,60)
title('Anomaly Percentage(%)')