NCLtips! (*'O')v

例：SSTトレンド †

1981年〜2015年の年平均SSTトレンド。
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

  in      = addfile("./HadISST_sst.nc","r")  ; HadISSTを使用
  sst     = in->sst
  time    = in->time
  lat     = in->latitude
  lon     = in->longitude
  YYYYMM  = cd_calendar(time,-1)             ; 時間の配列からYYYYMMを整数で取得

  ;;; 1981年〜2015年だけ取り出す
  t1      = ind(YYYYMM.eq.1981*100+1)
  t2      = ind(YYYYMM.eq.2015*100+12)
  sst    := sst(t1:t2,:,:)
  YYYYMM := YYYYMM(t1:t2)
 
  ;;; 年平均とトレンドの計算 
  sst_ann = month_to_annual_weighted(YYYYMM,sst,1)  ; 年平均(YYYYMMを使って各月の日数で重み付け)
  YYYY    = YYYYMM(::12)/100                ; 1981から2015までの配列を作る
                                            ;    ispan(1981,2015,1)でも同じこと
  trend   = regCoef_n(YYYY,sst_ann,0,0)     ; 回帰係数の計算

  ;;; t検定
  rstd    = reshape(trend@rstd,dimsizes(trend))   ; trend@rstdは1次元配列なのでtrendと同じ形に変形
  dof     = new(dimsizes(trend),integer)          ; 自由度の配列
  dof     = dimsizes(YYYY)-2                      ; ここでは簡単のためどこでも年数とした。もちろん実際の解析では検討が必要。
  cdl     = cdft_p(trend/rstd,dof)                ; t値(trend/rstd)からt分布の片側確率を計算
  cdl     = mask(cdl,ismissing(trend),False)      ; 陸地をマスキング
     ;;; cdl には，t分布を-∞から各地点のt値まで積分して得られる確率(0から1まで)が入る
  
  trend       = trend*10.                         ; 単位を変換

  ;;; 格子情報等の付加
  trend!0     = "lat" 
  trend!1     = "lon"
  trend&lat   = lat
  trend&lon   = lon
  trend@unit  = "~S~O~N~C/decade"
  copy_VarCoords(trend,cdl)
 

  ;;; 以降，お絵かきの設定
  wks  = gsn_open_wks("eps","sst_trend")
 
  res                 = True         ; トレンドのためのres
  res@gsnDraw         = False        ; plotを描かない
  res@gsnFrame        = False        ; WorkStationを更新しない
  res@cnLinesOn       = False
  res@cnFillOn        = True
  res@cnInfoLabelOn   = False
  res@tiMainString    = "SST trend (1981-2015)"
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF  = -0.3
  res@cnMaxLevelValF  =  0.3
  res@cnLevelSpacingF =  0.05
  res@cnFillPalette   = "BlWhRe"
  res@mpMaxLatF       =  60
  res@mpMinLatF       = -60
  res@mpCenterLonF    = 210
  res@gsnRightString  = ""
  res@pmLabelBarOrthogonalPosF = 0.2        ; カラーバーの位置を微調整
  ;;; タイトルや軸の文字の大きさを設定
  res@tmYLLabelFontHeightF = 0.016
  res@tmXBLabelFontHeightF = 0.016
  res@tiMainFontHeightF    = 0.024
  res@lbLabelFontHeightF   = 0.016
  ;;; カラーバーにタイトルをつける
  res@lbTitleOn            = True
  res@lbTitleString        = trend@unit
  res@lbTitlePosition      = "Right"
  res@lbTitleDirection     = "Across" 
  res@lbTitleFontHeightF   = 0.016
  
  plot = gsn_csm_contour_map_ce(wks,trend,res)   ; trendを描いたものを一旦plotに収める

  res2                 = True          ; 有意性のためのres
  res2@gsnDraw         = False         ; plotを描かない
  res2@gsnFrame        = False         ; WorkStationを更新しない
  ;;; あとでShadeLtGtContourを用いるため，あらゆるものをFalseにしておく
  res2@cnLinesOn       = False
  res2@cnLineLabelsOn  = False
  res2@cnFillOn        = False
  res2@cnInfoLabelOn   = False
  ;;; ShadeLtGtContourのために，等値線を念のため指定しておく
  res2@cnLevelSelectionMode = "ExplicitLevels"
  res2@cnLevels        = (/0.024,0.025,0.975,0.976/)

  dum  = gsn_csm_contour(wks,cdl,res2)           ; とりあえずcdlを描く
  dum  = ShadeLtGtContour(dum,0.025,6,0.975,17)
      ;;; 0.025のコンターより下を6番のハッチ，0.975のコンターより上を点々(17番)
      ;;; これは有意水準5％の両側検定に対応

  overlay(plot,dum)  ; 有意性を示したdumをplotに重ねる
  
  draw(plot)         ; ここでplotを描く
  frame(wks)         ; WorkStationの更新
   
end