请教：如何正确的对sst资料（ERSSTV3)做t检验？

Wziqi

做对几个EL Nino年事件的合成分析，需要用ncl进行t检验，之前用ncep的再分析资料做，结果可以判断显著性，但ncep的再分析资料不是很可靠，所以想换成ERSSTV3的sst资料再画一下，结果出不来通过95%显著水平的的区域了，但是我只是把脚本里的资料换了一下，别的都没改，为什么就会出错了？资料的差异不可能达到让这么大块的显著区域都没有的程度啊。我想了很久，感觉应该是ERSSTV3的资料把陆地和极地都是没有数值的，应该是这块的问题，ncl调试也说warning:ttest:encountered 1176 cases where var1 and/or equal to 0. 0utput set to missing values in these cases,之后又把sx和sy用了dim_num函数，但还是不对，请教大神如何对这种海温数据正确的做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"

f = addfile("sst.mnmean.nc","r")
sst = f->sst
delete(f)

startyear=1854                   ;; data begins 1854 Jan
index1=(1948-startyear)*12
index2=(2010-startyear)*12+11
epyears = (/1951,1958,1966,1973,1977,1998,2009/)
cpyears = (/1969,1978,1987,1995,2003,2010/)
epyears = epyears-1948                        
cpyears = cpyears-1948

sst2 = sst(index1:index2,:,:)    ;; chop the data

djf = month_to_season(sst2,"DJF")

ep = djf(1:dimsizes(epyears),:,:)      ;; copy the array dimensions
do i=0,dimsizes(epyears)-1
  ep(i,:,:) = djf(epyears(i),:,:)
end do
cp = djf(1:dimsizes(cpyears),:,:)      ;; copy the array dimensions
do i=0,dimsizes(cpyears)-1
  cp(i,:,:) = djf(cpyears(i),:,:)
end do

djf_mean = dim_avg_n_Wrap(djf,0)
djf_variance = dim_variance_n(djf,0)
ep_mean = dim_avg_n_Wrap(ep,0)
ep_variance = dim_variance_n(ep,0)
cp_mean = dim_avg_n_Wrap(cp,0)
cp_variance = dim_variance_n(cp,0)

siglvl = 0.05
iflag = False
probep = ep_mean   ;; copy meta
probcp = cp_mean
probep = ttest(ep_mean,ep_variance,dimsizes(epyears),djf_mean,djf_variance,dimsizes(djf(:,0,0)),iflag,False)
probcp = ttest(cp_mean,cp_variance,dimsizes(cpyears),djf_mean,djf_variance,dimsizes(djf(:,0,0)),iflag,False)

ep_mean = ep_mean-djf_mean
cp_mean = cp_mean-djf_mean

wks = gsn_open_wks("png","mean")
gsn_define_colormap(wks,"BlWhRe")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@cnInfoLabelOn = False
res@cnFillOn = True
res@cnLineLabelsOn = False
res@gsnSpreadColors = True
res@lbLabelBarOn = False
res@tiMainString = ""
res@gsnRightString = ""
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF = -2.
res@cnMaxLevelValF = 2.
res@cnLevelSpacingF = 0.4
res@mpCenterLonF = 180

res@gsnLeftString = "EP"
plot1 = gsn_csm_contour_map_ce(wks,ep_mean,res)
res@gsnLeftString = "CP"
plot2 = gsn_csm_contour_map_ce(wks,cp_mean,res)

res2 = res
res2@cnMonoFillPattern = False
res2@cnLevelSelectionMode = "ExplicitLevels"
res2@cnLevels = (/siglvl/)                       ;; set to significance level
res2@cnFillPatterns = (/3,-1/)
res2@gsnLeftString = ""
plot3 = gsn_csm_contour(wks,probep,res2)
overlay(plot1,plot3)
plot4 = gsn_csm_contour(wks,probcp,res2)
overlay(plot2,plot4)

resP = True
resP@gsnPanelLabelBar = True
gsn_panel(wks,(/plot1,plot2/),(/2,1/),resP)

超台十八

资料这里下的 v3b

http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html

超台十八

Wziqi 发表于 2015-2-2 10:02
看不到图，是不是没有传上来？

这样传上来吗

Wziqi

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

;************************************************
; Specify geographical region and time span (year-month start and end
;************************************************

  latS     = -90               
  latN     =  90
  lonL     =   0
  lonR     = 360

  ymStrt   = 194801     
  ymLast   = 201012
  

;************************************************
; Read from netCDF file: variable is type short...unpack
;************************************************
   diri   = "E:/data/"
   fili   = "sst.mnmean.nc"
   f      = addfile(diri+fili,"r")

   YYYYMM = cd_calendar( f->time, -1)

   iStrt  = ind(YYYYMM.eq.ymStrt)
   iLast  = ind(YYYYMM.eq.ymLast)

   x      = short2flt( f->sst(iStrt:iLast,{latS:latN},{lonL:lonR}) )

   printVarSummary(x)                            ; [time| 720]x[lat| 91]x[lon| 180]

   yyyymm = cd_calendar(x&time, -1)             ;转换时间格式
   yyyy   = yyyymm/100                          ;年份

   dimx = dimsizes(x)                           ;给各变量分配维数大小
   ntim = dimx(0)               
   nlat = dimx(1)
   mlon = dimx(2)

   year  = ispan(yyyy(0), yyyy(ntim-1), 1)       ;创建整型数组 （1948，。。。。2010）
   nyrs  = dimsizes(year)                         ;年数 63年

;************************************************
   xann=clmMonTLL(x)
   xano=calcMonAnomTLL(x,xann)
   season="DJF"
   xdjf  = month_to_season(xano ,season)                 ; 月平均资料转换为年平均资料
   printVarSummary(xdjf)
;************************************************
   xep=(/xdjf(3,:,:),xdjf(10,:,:),xdjf(18,:,:),xdjf(25,:,:),xdjf(29,:,:),xdjf(50,:,:),xdjf(61,:,:)/)   
   xepsat=dim_avg_n_Wrap(xep,0)
   xcp=(/xdjf(21,:,:),xdjf(30,:,:),xdjf(39,:,:),xdjf(47,:,:),xdjf(55,:,:),xdjf(62,:,:)/)
   xcpsat=dim_avg_n_Wrap(xcp,0)
copy_VarMeta(xdjf(0,:,:),xepsat)
  copy_VarMeta(xdjf(0,:,:),xcpsat)
;************************************************
;ttest
;************************************************
siglvl=0.05
xepave=xepsat
xcpave=xcpsat
yave=dim_avg_n(xdjf,0)
xepvar=dim_variance_n(xep,0)
xcpvar=dim_variance_n(xcp,0)
yvar=dim_variance_n(xdjf,0)
sxep=dim_num_n((.not.ismissing(xep),0)
sxcp=dim_num_n((.not.ismissing(xcp),0)
sy=dim_num_n((.not.ismissing(xdjf),0)
iflag= False               ; population variance similar
probep = ttest(xepave,xepvar,sxep, yave,yvar,sy,iflag, False)
probcp = ttest(xcpave,xcpvar,sxcp, yave,yvar,sy,iflag, False)
copy_VarMeta(xepsat,probep)  
copy_VarMeta(xcpsat,probcp)  
printVarSummary(probep)
;************************************************
; plotting parameters
;************************************************

   wks  = gsn_open_wks("png","E:/ENSO/ELsat")       ; specifies a ps plot
   colors=read_colormap_file("BlWhRe")
  printVarSummary(colors)
   res                       = True     
   res@gsnDraw=False
   res@gsnFrame=False
   res@gsnMaximize           = True             ; make large
   res@cnFillPalette        =colors(20:75,:)
    res@cnFillOn              = True             ; turn on color
   res@cnLinesOn             = False            ; turn off contour lines
   res@cnLineLabelsOn        = False            ; turn off contour line labels
;res@cnFillMode            = "RasterFill"
  
   res@cnLevelSelectionMode  = "ManualLevels"   ; set manual contour levels
   res@cnMinLevelValF        =  -2.4            ; set min contour level
   res@cnMaxLevelValF        =   2            ; set max contour level
   res@cnLevelSpacingF       =   0.2            ; set contour interval

   res@mpFillOn              = False            ; turn off default background gray
    res@mpCenterLonF          = 180

;  res@gsnCenterString       = year(0)+"-"+year(nyrs-1)
   res@gsnLeftString ="EP"
  ; res@tiMainString          = ""    ; fili
    res@lbLabelBarOn=False
;   res@cnLineColors=True
;   res@cnLinePalette="BlWhRe"
;*********************************************
  res2                       = True     
res2@gsnDraw         = False                      ; do not draw  
res2@gsnFrame        = False                      ; do not advance frame
   res2@gsnMaximize           = True             ; make large
;  res2@cnFillPalette        =colors(22:72,:)
   ; res2@cnFillOn              = True             ; turn on color
;   res2@cnLinesOn             = False            ; turn off contour lines
;   res2@cnLineLabelsOn        = False            ; turn off contour line labels  
  res2@cnLevelSelectionMode  = "ManualLevels"   ; set manual contour levels
   res2@cnMinLevelValF        = -1.4            ; set min contour level
res2@cnMaxLevelValF        = 0.6          ; set max contour level
  res2@cnLevelSpacingF       =   0.1           ; set contour interval
   res2@lbLabelBarOn=False
   res2@cnInfoLabelString=""
     
     res2@gsnRightString =""
   res2@gsnLeftString =""
   plot1 = gsn_csm_contour_map_ce(wks,xepsat,res)
  plot3 = gsn_csm_contour(wks,probep,res2)
   overlay(plot1,plot3)
    res@gsnLeftString ="CP"
      
   plot2 = gsn_csm_contour_map_ce(wks,xcpsat,res)
   plot4 = gsn_csm_contour(wks,probcp,res2)
   overlay(plot2,plot4)
    pres=True
    pres@gsnPanelLabelBar=True
    pres@pmLabelBarWidthF=0.8
    pres@lbBoxLinesOn=False
  gsn_panel(wks,(/plot1,plot2/),(/2,1/),pres)
  end

張顧.煒Zane

纬度调成-60~60，就不把极地的值算进来，-70~70也可以试试

Wziqi

張顧.煒Zane 发表于 2015-1-17 17:01
纬度调成-60~60，就不把极地的值算进来，-70~70也可以试试

还是不可以呢 和之前一样的问题 但是很奇怪的是我如果让函数输出t的值再查表就可以得到显著性了 真的是好奇怪啊

shiyan101

高手啊 ，能写出这样的代码

Wziqi

shiyan101 发表于 2015-1-19 19:00
高手啊 ，能写出这样的代码

没有啦，前面的数据读入和时间转换是用官网的例子改的。

shiyan101

请问是大气所的么

超台十八

解决了吗? 你用 djf 是有问题的，程序只算 jf

超台十八

xep=(/xdjf(3,:,:),xdjf(10,:,:),xdjf(18,:,:),xdjf(25,:,:),xdjf(29,:,:),xdjf(50,:,:),xdjf(61,:,:)/)   

   xcp=(/xdjf(21,:,:),xdjf(30,:,:),xdjf(39,:,:),xdjf(47,:,:),xdjf(55,:,:),xdjf(62,:,:)/)


1983 年是哪个?

Wziqi

shiyan101 发表于 2015-1-31 20:34
请问是大气所的么

不是的，坐标南京~