NCLtips! (*'O')v

例：降水量の気候値と経年変動標準偏差 †

降水量の気候値と経年変動の標準偏差を年平均および各季節平均について描いた。
標準偏差をシェード，気候値を等値線で示した。
気候値や標準偏差はCMAPの1979年～2015年に基づく。

スクリプト †

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"

;;; ここでは地図の領域等を指定する部分をprocedureとしてあらかじめ定義してみます
undef("set_map_area")
procedure set_map_area(res:logical,slat:numeric,nlat:numeric, \
                       wlon:numeric,elon:numeric,clon:numeric,cyc:logical)
begin
  res@mpMinLatF    = slat  ; 南限の緯度
  res@mpMaxLatF    = nlat  ; 北限の緯度
  res@mpMinLonF    = wlon  ; 西限の経度
  res@mpMaxLonF    = elon  ; 東限の経度
  res@mpCenterLonF = clon  ; 中心の経度
  res@gsnAddCyclic = cyc   ; 東西に周期的にするか？(True or False)
end
;;; ここまでがprocedureの定義。これを書いたスクリプトを別に用意してloadしても良い。

;;; ここからが本体
begin
  out     = "eps"   ; 絵の形式を変数として最初に与えられるようにしている

  in      = addfile("precip.mon.mean.nc","r")
  time    = in->time
  YYYYMM  = cd_calendar(time,-1)        ; 
  t1      = ind(YYYYMM.eq.197901)       ; 期間を決めて切り取るための変数
  t2      = ind(YYYYMM.eq.201512)       ;
  lat     = in->lat
  lon     = in->lon
  NY      = dimsizes(lat)
  NX      = dimsizes(lon)
  precip  = in->precip(t1:t2,:,:)       ; 時間は解析期間に対応する範囲のみ読む
  ann     = month_to_annual(precip,1)   ; 年平均の値
  precip  = runave_n(precip,3,0,0)      ; 3か月の移動平均を計算
  djf     = precip(0::12,:,:)           ; DJF平均の値
  mam     = precip(3::12,:,:)           ; MAM平均の値
  jja     = precip(6::12,:,:)           ; JJA平均の値
  son     = precip(9::12,:,:)           ; SON平均の値

  clm     = new((/5,NY,NX/),float)      ; 気候値が入る配列
  std     = new((/5,NY,NX/),float)      ; 標準偏差が入る配列
      ;;; 一番左の次元は，年間，DJF，MAM，JJA，SON平均を一つの配列に入れるため
 
  ;;; 平均や標準偏差を計算する
  clm(0,:,:) = dim_avg_n(ann,0)
  std(0,:,:) = dim_stddev_n(ann,0)
  clm(1,:,:) = dim_avg_n(djf,0)
  std(1,:,:) = dim_stddev_n(djf,0)
  clm(2,:,:) = dim_avg_n(mam,0)
  std(2,:,:) = dim_stddev_n(mam,0)
  clm(3,:,:) = dim_avg_n(jja,0)
  std(3,:,:) = dim_stddev_n(jja,0)
  clm(4,:,:) = dim_avg_n(son,0)
  std(4,:,:) = dim_stddev_n(son,0)

  ;;; 配列に座標情報を与える
  clm!0 = "season"
  clm!1 = "lat"
  clm!2 = "lon"
  clm&season = (/"Annual","DJF","MAM","JJA","SON"/)
  clm&lat = lat
  clm&lon = lon
  copy_VarCoords(clm,std)  ; clmの座標情報をstdにコピー
 
  ;;; 以降，お絵描き
  wks = gsn_open_wks(out,"precip")

  res = True  ; 標準偏差のプロット用(ベースとなるので地図に関するの設定もresに対して行う)
  res@gsnDraw   = False
  res@gsnFrame  = False
  res@cnFillOn  = True
  res@cnLinesOn = False
  res@cnInfoLabelOn  = False
  res@cnLineLabelsOn = False
  res@cnFillMode     = "AreaFill"
  res@cnFillPalette  = "CBR_wet"
  
  set_map_area(res,-30,30,30,290,180,False)  ; 上で定義したprocedureを使用
 
  res@tmXBLabelFontHeightF  = 0.025
  res@tmYLLabelFontHeightF  = 0.025
  res@pmTickMarkDisplayMode = "Always"       ; 縦軸と横軸の表示を「いい感じ」にする
  res@mpDataBaseVersion     = "MediumRes"    ; 地図の解像度
  res@mpGeophysicalLineThicknessF = 1
  res@mpGeophysicalLineColor = "green"

  res@cnLevelSelectionMode   = "ManualLevels"
  res@cnMinLevelValF         = 0.5
  res@cnMaxLevelValF         = 4.0
  res@cnLevelSpacingF        = 0.5
 
  res@lbLabelBarOn       = False
  res@gsnLeftStringFontHeightF = 0.035
  res@gsnRightString = ""

  resx = True  ; 気候値のプロット用(あとで標準偏差の図に重ねる)
  resx@gsnDraw   = False
  resx@gsnFrame  = False
  resx@cnFillOn  = False
  resx@cnLinesOn = True
  resx@lbLabelBarOn   = False
  resx@cnInfoLabelOn  = False
  resx@cnLineLabelsOn = False
  resx@cnLineColor    = "red"
  resx@cnLineThicknessF = 0.5
  resx@cnLevelSelectionMode   = "ManualLevels"
  resx@cnMinLevelValF         = 3.
  resx@cnMaxLevelValF         = 15.
  resx@cnLevelSpacingF        = 3.
  resx@cnLineDrawOrder = "PostDraw"  ; 等値線が地図よりも上にくるように
  resx@gsnLeftString  = ""
  resx@gsnRightString = ""
 
  ;;; 5枚分の図が入る配列×2を用意
  plot  = new(5,graphic)  ; 標準偏差用
  plotx = new(5,graphic)  ; 気候値用

  ;;; 年間および各季節平均の絵を描いていく
  do i = 0, 4
    res@gsnLeftString = clm&season(i)
    plot(i)  = gsn_csm_contour_map_ce(wks,std(i,:,:),res)
    plotx(i) = gsn_csm_contour(wks,clm(i,:,:),resx)
    overlay(plot(i),plotx(i))
  end do

  resP = True                           ; パネルプロット
  resP@gsnPanelRowSpec  = True          ; 各行の枚数で並び方を指定
  resP@gsnPanelLabelBar = True          ; カラーバーはまとめて一つに 
  resP@lbLabelFontHeightF = 0.015
  resP@lbTitleOn          = True
  resP@lbTitleString      = "  mm/day"
  resP@lbTitleFontHeightF = 0.015
  resP@lbTitlePosition    = "Right"
  resP@lbTitleDirection   = "Across"
  resP@pmLabelBarOrthogonalPosF = -0.03
  resP@txString           = "Precipitation / climatology & Standard deviation"
  resP@gsnPanelXWhiteSpacePercent = 5   ; パネル同士の間隔(ヨコ)
  resP@gsnPanelYWhiteSpacePercent = 10  ; パネル同士の間隔(タテ)
  gsn_panel(wks,plot,(/1,2,2/),resP)    ; 1行目は年平均のみ，2行目と3行目は2枚ずつ
  
end