关于用ncl做回归分析的显著性检验

lleoiu

之前没看到rc和rcc区别.
我想应该用tval = onedtond(rc@tval , dimsizes(rc))  
用三维的来定tval比较准确.
plot的时候,可以用gsn_csm_pres_hgt函数.

xhb3256789

学习学习

李小毛123

lleoiu 发表于 2017-11-5 18:24
之前没看到rc和rcc区别.
我想应该用tval = onedtond(rc@tval , dimsizes(rc))  
用三维的来定tval比较准 ...

好的，那我先改改看。

伍孚ms

码一下自己改的回归及显著性检验脚本：
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"
load "/home/sma/data/work/Iran_pre/FP/shapefile_utils.ncl"

begin

  f  = addfile("/data1/data/ttxie/CRU/cru_ts3.26.1901.2017.pre.dat.nc", "r")
  f1 = addfile("/data1/data/ttxie/ERA/pressure/2.5_1979-2017/2.5_2.5 1979-2017.nc", "r")
  ;   varnames = getfilevarnames(f1)
  ; if(.not.any(ismissing(varnames))) then
  ;   do i=0,dimsizes(varnames)-1
  ;     printFileVarSummary (f1,varnames(i))
  ;   end do
  ; end if
  ; exit
  ymstart = 197901
  ymlast  = 201712
  YYYYMM  = cd_calendar(f->time,-1)

  iStart  = ind(YYYYMM.eq.ymstart)
  iLast   = ind(YYYYMM.eq.ymlast)

  p = f->pre(iStart:iLast,{35.25:45},{46.25:80})
  p&lat@units = "degrees_north"
  p&lon@units = "degrees_east"
  printVarSummary(p)

  hgt = f1->z(:,{700},:,:)
  hgt&lat@units = "degrees_north"
  hgt&lon@units = "degrees_east"
  printVarSummary(hgt)

  hgt_DJF   = month_to_season(hgt, "DJF")
  printVarSummary(hgt_DJF)
  p_winter = month_to_season(p,"DJF")

  pre_areaavg_DJF = wgt_areaave_Wrap(p_winter,1.0,1.0,0)
  p_winter_st = dim_standardize_Wrap(pre_areaavg_DJF, 1)
  p_winter_st_1 = tofloat(p_winter_st)

  R = regCoef(p_winter_st_1,hgt_DJF({lat|:},{lon|:},{time|:}))
  copy_VarMeta(hgt_DJF(0,:,:),R)
  R!0 = "lat"
  R!1 = "lon"
  lat  = R&lat
  lon  = R&lon
  R&lon@units = "degrees_east"
  R&lat@units = "degrees_north"
  printVarSummary(R)

  prob    = student_t(R@tval, R@nptxy-2)
  printVarSummary(prob)

  ;copy_VarMeta(R(:,:), prob(:,:))
  printVarSummary(prob)
  P = onedtond(prob, (/73,144/))
  copy_VarMeta(R(:,:), P(:,:))
  P@_FillValue  = -999.99
  printVarSummary(P)

  R1 = new((/73,144/), float)
  R2 = new((/73,144/), float)

  R1 = where(R.gt.0, R, 0)
  R2 = where(R.gt.0, 0, R)
  copy_VarMeta(R,R1)
  copy_VarMeta(R,R2)
  print(R1)


;*****************************************************         
wks   = gsn_open_wks ("pdf","pre_hgt_reg")                 
  
  ;显著性的填色图及回归场的等值线分开绘制，先设置两图的公用绘图参数      
  res                   = True   
  res@gsnAddCyclic      = True
  res@gsnDraw           = False         
  res@gsnFrame          = False        
  res@gsnLeftString     = ""
  res@gsnRightString    = ""  

  res@gsnTickMarksOn = False ; 关闭经度标签。虽然默认是绘制经度标签，但由于其经度单位没有“度”符号，因此关闭其经度标签，若绘制标准的
                             ; 经度标签，可利用函数gsn_add_text以及文本符号“~”进行手动添加
  resc = res          ;¸创建resc，用以绘制回归场的等值线      
  resc2 = res
  res@mpDataBaseVersion           = "Ncarg4_0" ; or "Ncarg4_1"
  res@mpOutlineBoundarySets       = "National"
  res@mpLandFillColor             = "white"
  res@mpInlandWaterFillColor      = "white"
  ;vcres@mpUSStateLineThicknessF     = 1  
  res@mpGeophysicalLineThicknessF = 1
  res@mpMaxLatF                   = 90        
  res@mpMinLatF                   = -10
  res@mpMaxLonF                   = 120
  res@mpMinLonF                   = -90  

  res@mpShapeMode                 = "FreeAspect"
  res@vpHeightF                   =  0.55   
  res@vpWidthF                    =  0.8

  res@pmTickMarkDisplayMode      = "Always"
  res@tmXTOn                     = False
  res@tmYROn                     = False
  
  res@cnLevelSelectionMode  = "ExplicitLevels"
  res@cnLevels              = (/0.1/) ; -20与20均对应着0.05置信度，但前者对应负值异常，后者对应正值异常
  res@cnFillColors          = (/"gray","white"/)
  res@cnFillOn              = True
  res@lbLabelBarOn          = False   ; 在表示显著性时，一般不需要绘制等值线线条及其数值，因此关闭以下几项         
  res@cnLinesOn             = False   ;  
  res@cnInfoLabelOn         = False   ;
  res@cnLineLabelsOn        = False   ;
   
  base = gsn_csm_contour_map(wks, P, res)

  ;; 回归场等值线设置
  resc@cnLevelSelectionMode  = "ManualLevels"
  resc@cnMaxLevelValF         = 150
  resc@cnMinLevelValF         = -0  
  resc@cnLevelSpacingF        = 10
  resc@cnLineDashPattern    = 0
  ;resc@cnLevels              =  1.*ispan(0,100,1)   
  resc@cnFillOn              = False   
  resc@cnInfoLabelOn         = False
  resc@gsnContourZeroLineThicknessF = 0.
  resc@cnLineThicknessF             = 1.5
  resc@cnLineLabelsOn               = False
  plot = gsn_csm_contour(wks,R1,resc)

  
  overlay(base, plot) ; 图层叠加

  ; ;; 负相关虚线
  resc2@cnLineDashPattern    = 1
  resc2@cnLevelSelectionMode  = "ManualLevels"
  resc2@cnMaxLevelValF         = 0
  resc2@cnMinLevelValF         = -80  
  resc2@cnLevelSpacingF        = 12
  ;resc2@cnLevels              =  1.*ispan(-100,0,1)   
  resc2@cnFillOn              = False   
  resc2@cnInfoLabelOn         = False
  resc2@gsnContourZeroLineThicknessF = 0.
  resc2@cnLineThicknessF             = 1.5
  resc2@cnLineLabelsOn               = False
  plot2 = gsn_csm_contour(wks,R2,resc2)
   
  overlay(base,plot2)
  
  draw(base)
  frame(wks)

end

阿羊出没请注意

伍孚ms 发表于 2020-5-10 11:25
码一下自己改的回归及显著性检验脚本：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load ...

赞！！！

vicky-zz

李小毛123 发表于 2017-11-5 16:31
tval = onedtond(rcc@tval , dimsizes(rcc))  
printVarSummary(rcc)
tval@FillValue=  -9.96921e+36 ...

请问一下这个问题解决了，我画了一下打点区域也有问题

everglow

请问一下我要是用循环画很多个rc，但属性tval只有最后一个，我要怎么对这么多个rc检验呢？