- 积分
- 323
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-12-15
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 Haryu 于 2017-4-8 17:38 编辑
本人做了EP nino几年降水合成，但t检验通过的少，不知为什么，请大神我看看，十分感谢，程序如下：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
; ===========================================
begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
latS = 10.
latN = 20.
lonL = 50.
lonR = 180.
ym82s = 19820416
ym82l = 19820730
ym87s = 19870416
ym87l = 19870730
ym97s = 19970416
ym97l = 19970730
ym15s = 20150416
ym15l = 20150730
;Normal year
ym81s = 19810416
ym81l = 19810730
ym84s = 19840416
ym84l = 19840730
ym85s = 19850416
ym85l = 19850730
ym89s = 19890416
ym89l = 19890730
ym90s = 19900416
ym90l = 19900730
ym93s = 19930416
ym93l = 19930730
ym96s = 19960416
ym96l = 19960730
ym99s = 19990416
ym99l = 19990730
ym00s = 20000416
ym00l = 20000730
ym01s = 20010416
ym01l = 20010730
ym08s = 20080416
ym08l = 20080730
ym11s = 20110416
ym11l = 20110730
ym12s = 20120416
ym12l = 20120730
ym13s = 20130416
ym13l = 20130730
leftname = "(a)"
title = "The precipitation of EP El Nino"
outname = "precipitationet2"
; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
f = addfile ("precip.pentad.mean.nc","r")
time = f->time
ym = cd_calendar(time,-2)
;print(ym)
precip = f->precip(:,:,:)
;precip = dim_rmvmean_n_Wrap(precip1,0)
;printVarSummary(precip)
;print(precip(0,{latS:latN},{lonL:lonR}))
precipe = new((/88,72,144/),float)
y82 = ind(ym.ge.ym82s.and.ym.le.ym82l)
y87 = ind(ym.ge.ym87s.and.ym.le.ym87l)
y97 = ind(ym.ge.ym97s.and.ym.le.ym97l)
y15 = ind(ym.ge.ym15s.and.ym.le.ym15l)
precipe(0:21,:,:) = precip(y82,:,:)
precipe(22:43,:,:) = precip(y87,:,:)
precipe(44:65,:,:) = precip(y97,:,:)
precipe(66:87,:,:) = precip(y15,:,:)
precipeave = new((/22,72,144/),float)
do i = 0,21,1
precipeave(i,:,:) = dim_avg_n_Wrap(precipe(i:87:22,:,:),0)
end do
; trick to copy coord information
tme = ispan(22,43,1)
;print(tme)
precipeave!0 = "time"
precipeave&time = tme
;calculate varprecipavg
varprecipeave = new((/22,72,144/),float)
do j = 0,21,1
varprecipeave(j,:,:) = dim_variance_n_Wrap(precipe(j:87:22,:,:),0)
end do
copy_VarCoords(precipeave, varprecipeave) ; trick to copy coord information
precipeavedim = dim_avg_n_Wrap(precipeave(:,{latS:latN},:),1)
; read in Normal year
y81 = ind(ym.ge.ym81s.and.ym.le.ym81l)
y84 = ind(ym.ge.ym84s.and.ym.le.ym84l)
y85 = ind(ym.ge.ym85s.and.ym.le.ym85l)
y89 = ind(ym.ge.ym89s.and.ym.le.ym89l)
y90 = ind(ym.ge.ym90s.and.ym.le.ym90l)
y93 = ind(ym.ge.ym93s.and.ym.le.ym93l)
y96 = ind(ym.ge.ym96s.and.ym.le.ym96l)
y99 = ind(ym.ge.ym99s.and.ym.le.ym99l)
y00 = ind(ym.ge.ym00s.and.ym.le.ym00l)
y01 = ind(ym.ge.ym01s.and.ym.le.ym01l)
y08 = ind(ym.ge.ym08s.and.ym.le.ym08l)
y11 = ind(ym.ge.ym11s.and.ym.le.ym11l)
y12 = ind(ym.ge.ym12s.and.ym.le.ym12l)
y13 = ind(ym.ge.ym13s.and.ym.le.ym13l)
precipn = new((/308,72,144/),float)
precipn(0:21,:,:) = precip(y81,:,:)
precipn(22:43,:,:) = precip(y84,:,:)
precipn(44:65,:,:) = precip(y85,:,:)
precipn(66:87,:,:) = precip(y89,:,:)
precipn(88:109,:,:) = precip(y90,:,:)
precipn(110:131,:,:) = precip(y93,:,:)
precipn(132:153,:,:) = precip(y96,:,:)
precipn(154:175,:,:) = precip(y99,:,:)
precipn(176:197,:,:) = precip(y00,:,:)
precipn(198:219,:,:) = precip(y01,:,:)
precipn(220:241,:,:) = precip(y08,:,:)
precipn(242:263,:,:) = precip(y11,:,:)
precipn(264:285,:,:) = precip(y12,:,:)
precipn(286:307,:,:) = precip(y13,:,:)
precipnave = new((/22,72,144/),float)
do x = 0,21,1
precipnave(x,:,:) = dim_avg_n_Wrap(precipn(x:307:22,:,:),0)
end do
;calculate varprecipavg
varprecipnave = new((/22,72,144/),float)
do y = 0,21,1
varprecipnave(y,:,:) = dim_variance_n_Wrap(precipn(y:307:22,:,:),0)
end do
copy_VarCoords(precipeave, precipnave) ; trick to copy coord information
copy_VarCoords(precipeave, varprecipeave) ; trick to copy coord information
;==============================================
;Returns an estimate of the statistical significance and, optionally, the t-values.
;=============================================
iflag = True
prob = ttest(precipeave,varprecipeave,4, precipnave,varprecipeave,14, iflag, False )
;print(prob(0,:,:,:))
copy_VarCoords(precipeave, prob) ; trick to copy coord information
;printVarSummary(prob)
probdim = dim_avg_n_Wrap(prob(:,{latS:latN},:),1)
;print(probdim)
;==============================================
; create color plot
;=============================================
wks = gsn_open_wks ("png", outname ) ; open ps file
gsn_define_colormap(wks,"CBR_wet") ; choose colormap
res = True ; plot mods desired
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@gsnMaximize = True ; large format
res@tiMainString = title ; title
res@tmXBMode = "Explicit"
res@tmXBValues = (/55,70, 90, 110, 130, 150,170/)
res@tmXBLabels = (/"55E","70E","90E","110E","130E","150E","170E"/)
res@tmYLMode = "Explicit"
res@tmYLValues = (/22.,24., 26., 28.,30.,32.,34.,36.,38.,40.,42./)
res@tmYLLabels = (/"22","24","26","28","30","32","34","36","38","40","42"/)
res@tmXTOn = False
res@tmYROn = False
res@tmXBMinorOn = False
res@tmYLMinorOn = False
res@txFontHeightF = 0.005
res@pmLabelBarHeightF = 0.035 ; height of the LabelBar
res@cnInfoLabelOn = False
res@vpHeightF = 0.5
pres =res
res@cnFillOn = False ; turn on color fill
res@cnLineLabelsOn = True ; True is default
res@cnMinLevelValF = 5.
;res@cnLevelSelectionMode = "ManualLevels"
res@gsnLeftString = leftname
res@gsnLeftStringFontHeightF = 0.02
res@gsnContourZeroLineThicknessF = 2. ; doubles thickness of zero contour
res@gsnContourNegLineDashPattern = 1 ; sets negative contours to dash pattern 1
; set symmetric plot min/max
pres@cnLinesOn = False ; True is default
pres@cnLineLabelsOn = False ; True is default
opt = True
opt@gsnShadeLow = 8
plot1 = gsn_csm_hov(wks, precipeavedim(:,{lonL:lonR}),res)
plot2 = gsn_csm_hov(wks, probdim(:,{lonL:lonR}),pres)
plot2 = gsn_contour_shade(plot2,0.1,-999,opt)
overlay(plot2,plot1)
draw(plot2)
frame(wks)
end
|
-
合成的降水图，阴影为通过0.1检验
-
挑选出的四年EP elnino年的降水异常图
-
挑选的非el nino年的降水异常图
|
窥视卡
雷达卡
发表于 2017-4-8 17:29:28
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
