用SOI和PDO指数回归热带和北太平洋海盆SST并作显著性检验

HDdd

闲来无事 分享一下回归分析以及显著性检验的脚本

   ;load"$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"

;***********************************************************

    begin

    ;duqu1987-2016desoi&pdozhishu

    diri = "/cygdrive/G/ncl/qhtj/"

    soizs =asciiread(diri+"SOI-1951-CURRENT.txt", -1, "float")

    printVarSummary(soizs)

    ;duqu1987-2016sst

    fili = "sst.mon.mean.nc"

    f = addfile(diri+fili,"r")

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

    print(yyyymm(0))

    ymstrt  = 198701     

    ymlast  = 201612

    istrt = ind(yyyymm.eq.ymstrt)

    ilast = ind(yyyymm.eq.ymlast)

    temp  = f->sst;

    printVarSummary(temp);

    temper = temp(istrt:ilast,{:},{:})

    printVarSummary(temper)

    sst =temper(lat|:,lon|:,time|:);

    printVarSummary(sst)

    ; jisuanhuiguixishubingzuojianyan

    rc = regCoef(soizs,sst);返回值为tval（t值统计量tval=rc@tval）和nptxy（自由度df=rc@nptxy-2）

    printVarSummary(rc);

    rc@long_name = "Trend"

    rc!0 = "lat"

    rc!1 = "lon"

    rc&lat = sst&lat

    rc&lon = sst&lon

    printVarSummary(rc)

    tval = onedtond(rc@tval,dimsizes(rc));

    df = onedtond(rc@nptxy, dimsizes(rc)) - 2

    b = tval;betainc(x, a, b)的x，a，b的维度必须一致

    b = 0.5

    b!0 = "lat"

    b!1 = "lon"

    b&lat = sst&lat

    b&lon = sst&lon

    printVarSummary(b)

    prob = betainc(df/(df+tval^2),df/2.0,b);

    prob!0 = "lat"

    prob!1 = "lon"

    prob&lat = sst&lat

    prob&lon = sst&lon

  

    wks = gsn_open_wks("png",diri+"soiIII")

   gsn_define_colormap(wks,"BlueRedGray")  

;************************************************

; plotting parameters（画相关系数rc）

;************************************************

   

   res                       = True

   res@gsnDraw               = False

   res@gsnFrame              = False   

   res@gsnMaximize           = True             ; make large

   res@gsnAddCyclic          = False

   

   res@cnFillOn              = True             ; turn on color

   res@cnLinesOn             = False            ; turn off contour lines

   res@cnLineLabelsOn        = False            ; turn off contour linelabels

   ;;res@cnFillMode            = "RasterFill"

   res@cnLevelSelectionMode  = "ManualLevels"   ; set manual contour levels

   res@cnMinLevelValF        = -1              ; set min contourlevel

   res@cnMaxLevelValF        =  1              ; set max contourlevel

   res@cnLevelSpacingF       =  0.1            ; set contourinterval

  

   res@mpOutlineBoundarySets       = "Geophysical" ;

   res@mpMinLatF                   = -30.

   res@mpMaxLatF                   = 60.

   res@mpMinLonF                   = 105.

   res@mpMaxLonF                   = 265.

   res@mpCenterLonF                = 180.

   res@mpGeophysicalLineThicknessF = 2.         ; double the thickness of geophysicalboundaries

   res@mpNationalLineThicknessF    = 2.       ; double the thickness of national boundaries

   res@mpOutlineOn                 = True       ; turn off outline

   res@mpFillOn                    = False    ; turn off map fill   

   res@tiMainString          = "Regression Test of SOI andSST(1987-2016)"

  

   base = gsn_csm_contour_map(wks,rc,res)  

;============================================<

;add 画通过显著性检验的区域

;============================================<

   res2 = True

   res2@gsnDraw = False;

   res2@gsnFrame =False;

   

   res2@cnFillOn = True

   res2@cnLinesOn = False;

   res2@cnLineLabelsOn = False

   res2@cnInfoLabelOn = False;

   res2@lbLabelBarOn = False;

   res2@cnMonoFillPattern = False

   res2@cnLevelSelectionMode ="ExplicitLevels"

   res2@cnLevels = (/0.05/);

   res2@cnFillPatterns = (/7,-1/);

   res2@cnFillColors = (/1,-1/);

   res2@gsnLeftString = ""

   

   plot = gsn_csm_contour(wks,prob,res2)

   overlay(base, plot)

   draw(base)

   frame(wks)

   end

人生别来无恙

楼主，你好，我想请教一下脚本当中的这几行：
“   res2@cnLevels = (/0.05/);set to significancelevel设置显著性水平
   res2@cnFillPatterns = (/7,-1/);-1表示透明,res2@cnLevels =(/0.05/)把区域分成了>0.05（用7形状填充）和<0.05（用-1形状填充）两个区域
   res2@cnFillColors = (/1,-1/);-1表示透明,>0.05的区域用1填色，<0.05的区域用-1填色 ”
按照这样的设置，应该是>0.05的用-1填充，<0.05的用形状7和颜色1填充吧，你是不是注释反了啊

HDdd

人生别来无恙 发表于 2020-4-4 17:22
楼主，你好，我想请教一下脚本当中的这几行：
“   res2@cnLevels = (/0.05/);set to significancelevel设 ...

好像是的哈==

我头发乱了吗

请问楼主，这里的色标数值有意义吗？

HDdd

我头发乱了吗 发表于 2020-5-30 22:34
请问楼主，这里的色标数值有意义吗？

你说的哪一句？

我头发乱了吗

HDdd 发表于 2020-5-31 09:33
你说的哪一句？

图中的色标，这里的颜色深浅有含义吗？是有单位的吗？谢谢！

李小毛123

你好，我用的代码一毛一样，只是用来回归的时间系数不一致，为什么这个检验打出来的点感觉有问题

psfhy

你好，请问您SOI指数是怎么计算的