AMS 2024 Tutorial

plot_cubed.ncl

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

procedure add_box(wks,plot,ll[2],width,cname,oval)
local xbox, ybox, gsres, dumstr
begin
  xbox = (/ll(0),ll(0)+width,ll(0)+width,ll(0),ll(0)/)
  ybox = (/ll(1),ll(1),ll(1)+width,ll(1)+width,ll(1)/)

  gsres                = True
  gsres@gsFillColor    = cname
  gsres@gsFillOpacityF = oval
  dumstr               = unique_string("gon")
;
; Adding it as an attribute is a sneaky way to
; make sure it "lives" outside this procedure.
;
  plot@$dumstr$ = gsn_add_polygon(wks,plot,xbox,ybox,gsres)
end

begin
  ;--- parameters ---
  ntiles = 6

  ;--- select variable ---
  ;vname = "evap"
  ;vname = "tskin"
  vname = "hflx" 
  
  ;--- select time index ---
  tid = 11 

  ;--- turn off missing value check ---
  setfileoption("nc", "MissingToFillValue", False)

  ;--- loop over tiles ---
  do i = 0, ntiles-1
    ;--- open file ---
    lsA = systemfunc ("ls ufs.cpld.lnd.out.*.tile"+sprinti("%d", (i+1))+".nc")
    ncA = addfiles(lsA, "r")
    ListSetType(ncA, "cat")
    ncB = addfile("INPUT/oro_data.tile"+sprinti("%d", (i+1))+".nc", "r")

    ;--- read fields ---
    dum = ncA[:]->$vname$(tid,:,:)
    land_frac = ncB->land_frac

    ;--- read coordinate data ---
    lon2d = (/ ncA[0]->grid_xt /)
    lat2d = (/ ncA[0]->grid_yt /)

    ;--- string for label bar ---
    lstr = vname+" ["+dum@units+"]"

    time = ncA[:]->time(tid)
    time@calendar = "standard"
    time_str = ut_string(time, "%D-%N-%Y_%H:%M")

    ;--- create data array ---
    if (i .eq. 0) then
      dims = dimsizes(lon2d)
      nlat = dims(0)
      nlon = dims(1)
      data = new((/ntiles, nlat, nlon /), "double", dum@missing_value)
      lmsk = new((/ ntiles, nlat, nlon /), "integer")
      lon  = new((/ ntiles, nlat, nlon /), typeof(lon2d))
      lat  = new((/ ntiles, nlat, nlon /), typeof(lat2d))
    end if

    ;--- store coordinate arrays ---
    lon(i,:,:) = (/ lon2d /)
    lat(i,:,:) = (/ lat2d /)
    if (i .eq. 0 .or. i .eq. 5) then
      lon(i,:,:) = where(lon(i,:,:) .ge. 180. .or. lon(i,:,:) .eq.360., lon(i,:,:)-360, lon(i,:,:) )
    end if

    ;--- fill common variable ---
    data(i,:,:) = (/ dum /)
    lmsk(i,:,:) = (/ where(land_frac .gt. 0, 1, 0) /)
    lmsk(i,:,:) = (/ mask(lmsk(i,:,:), lmsk(i,:,:) .gt. 0, True) /)

    ;--- special care for following tiles ---
    if any((i .eq. (/ 3, 4 /))) then
      data(i,:,:) = (/ transpose(data(i,:,:)) /)
      data(i,:,:) = data(i,::-1,:)
      lmsk(i,:,:) = (/ transpose(lmsk(i,:,:)) /)
      lmsk(i,:,:) = lmsk(i,::-1,:)
      lon(i,:,:) = (/ transpose(lon(i,:,:)) /)
      lon(i,:,:) = lon(i,::-1,:)
      lat(i,:,:) = (/ transpose(lat(i,:,:)) /)
      lat(i,:,:) = lat(i,::-1,:)
    end if
    if any((i .eq. (/ 2 /))) then
      data(i,:,:) = (/ transpose(data(i,:,:)) /)
      data(i,:,:) = data(i,:,::-1)
      lmsk(i,:,:) = (/ transpose(lmsk(i,:,:)) /)
      lmsk(i,:,:) = lmsk(i,:,::-1)
      lon(i,:,:) = (/ transpose(lon(i,:,:)) /)
      lon(i,:,:) = lon(i,:,::-1)
      lat(i,:,:) = (/ transpose(lat(i,:,:)) /)
      lat(i,:,:) = lat(i,:,::-1)
    end if

    ;--- maskout data ---
    data(i,:,:) = mask(data(i,:,:), lmsk(i,:,:) .eq. 1, True)
    data(i,:,:) = mask(data(i,:,:), .not. ismissing(lmsk(i,:,:)), True)

    ;--- delete temporary variables ---
    delete([/ dum, lat2d, lon2d, land_frac /])
  end do

  ;--- create plot ---
  wks_type = "pdf"
  wks_type@wkOrientation = "landscape"
  wks = gsn_open_wks(wks_type, "plot_cubed_"+time_str+"_"+vname)

  res1 = True
  res1@gsnDraw = False
  res1@gsnFrame = False
  res1@gsnAddCyclic = False
  res1@cnFillOn = True
  res1@cnLinesOn = False
  res1@cnLineLabelsOn = False
  res1@gsnLeftString = ""
  res1@gsnRightString = ""
  size = 0.23
  res1@vpWidthF = size
  res1@vpHeightF = size
  res1@cnInfoLabelOn = False
  res1@lbLabelBarOn = False
  res1@cnFillMode = "RasterFill"
  res1@cnLevelSelectionMode = "ExplicitLevels"
  res1@cnMissingValFillColor = "white"
  res1@tmXBOn = False
  res1@tmXTOn = False
  res1@tmYLOn = False
  res1@tmYROn = False

  res2 = res1
  res2@cnFillOn = False
  res2@cnLinesOn = True
  res2@cnLineLabelsOn = False
  res2@cnInfoLabelOn = False
  res2@cnLineThicknessF = 0.5
  res2@cnLevels = (/ 0, 1 /) 

  ;--- set levels, variable specific ---
  res1@cnFillPalette = "NCV_banded"
  if (any(vname .eq. (/ "soilt1", "tg3"/))) then
    data = data-273.15 ; K -> degC
    res1@cnLevels = fspan(-30, 30, 61)
  end if
  if (any(vname .eq. (/ "evap"/))) then
    res1@cnLevels = fspan(-50, 200, 25)
  end if
  if (any(vname .eq. (/ "hflx"/))) then
    res1@cnLevels = fspan(-100, 300, 41)
  end if

  ;--- text resource ---
  txres = True
  txres@gsnFrame = False
  txres@txJust = "CenterCenter"

  ;--- box resource ---
  gnres = True
  gnres@gsFillColor = -1
  gnres@gsEdgesOn = True
  gnres@gsEdgeColor = "grey"

  ;--- tile order --- 
  ;   2
  ; 4 0 1 3
  ;   5
  xpos = tofloat((/ 1, 2, 1, 3, 0, 1 /))*size+0.05
  ypos = 0.9-tofloat((/ 1, 1, 0, 1, 1, 2 /))*size

  do k = 0, ntiles-1
    ;--- set position ---
    res1@vpXF = xpos(k)
    res1@vpYF = ypos(k) 
    res2@vpXF = xpos(k)
    res2@vpYF = ypos(k) 
    print(k+" "+res1@vpXF+" "+res1@vpYF)

    ;--- add labelbar ---
    if (k .eq. 1) then
      res1@lbLabelBarOn = True
      res1@pmLabelBarWidthF = size*1.8
      res1@pmLabelBarHeightF = 0.03
      res1@pmLabelBarParallelPosF = 0.97
    else
      res1@lbLabelBarOn = False
    end if

    ;--- add plot label ---
    res1@gsnStringFont = 22
    res1@gsnStringFontHeightF = 0.013
    res1@gsnLeftString = ""
    if (k .eq. 4) then
      res1@gsnLeftString = vname
    end if

    ;--- plot ---
    plot1 = gsn_csm_contour(wks, data(k,:,:), res1)
    draw(plot1)

    ;--- add box around each plot ---
    x_box = (/ 0.0, size, size, 0.0, 0.0 /)+xpos(k)
    y_box = (/ 0.0, 0.0, size, size, 0.0 /)+ypos(k)-size
    gsn_polygon_ndc(wks, x_box, y_box, gnres)

    ;--- add tile number ---
    txres@txPosXF = xpos(k)+0.01
    txres@txPosYF = ypos(k)-0.01
    txres@txAngleF = 0
    txres@txFont = 22
    txres@txFontHeightF = 0.015
    text = gsn_create_text(wks, sprinti("%d", k+1), txres)
    draw(text)

    ;--- add labelbar string ---
    if (k .eq. 3) then
      txres@txPosXF = xpos(k)+0.01
      txres@txPosYF = ypos(k)-size-0.04
      txres@txFontHeightF = 0.011
      text = gsn_create_text(wks, lstr, txres)
      draw(text)
    end if
  end do

 ;--- add date string ---
 txres@txPosXF = xpos(3)+size
 txres@txPosYF = ypos(3)+0.01
 txres@txJust = "centerright"
 text = gsn_create_text(wks, time_str, txres)
 draw(text)

 ;--- create plot ---
 frame(wks)
end

plot_cubed_ts.ncl

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
begin
  ;--- arguments ---
  ntiles = 6
  vname = "evap"
  ;vname = "hflx"
  ;vname = "tskin"

  ;--- coordinates of predefined land grid point for each tile ---
  ;--- different set of grid points can be defined in here --- 
  xc = (/ 79, 46, 78, 78, 10, 3 /)
  yc = (/ 47, 83, 12, 12, 42, 77 /)

  ;--- turn off missing value check ---
  setfileoption("nc", "MissingToFillValue", False)

  ;--- loop over tiles and prepare data ---
  do i = 0, ntiles-1
    ;--- open file ---
    lsA = systemfunc ("ls ufs.cpld.lnd.out.*.tile"+sprinti("%d", (i+1))+".nc")
    ncA = addfiles(lsA, "r")
    ListSetType(ncA, "cat")

    ;--- get time ---
    time = ncA[:]->time
    ntime = dimsizes(time)
    stime = cd_string(time, "%d %c - %H:%M")

    ;--- read fields ---
    dum = ncA[:]->$vname$(:,yc(i),xc(i))

    ;--- create data array ---
    if (i .eq. 0) then
      data = new((/ ntiles, ntime /), "double")
    end if

    data(i,:) = (/ dum /)

    ;--- delete temporary variables ---
    delete([/ ncA, dum /])
  end do

  ;--- plot data ---
  wks  = gsn_open_wks("newpdf", "plot_cubed_ts_"+vname)
  gsn_define_colormap(wks, "matlab_jet")

  res = True
  res@gsnDraw = False
  res@gsnFrame = False
  res@vpHeightF = 0.2
  res@vpWidthF = 0.4

  res@tmXTOn = False
  res@tmXTLabelsOn = False
  res@tmYROn = False
  res@tmYRLabelsOn = False
  res@tmYLMinorOn = False
  res@tmXBLabelAngleF = 90.0
  res@tmXBLabelFontHeightF = 0.01
  res@tmYLLabelFontHeightF = 0.01
  colors = (/ "black", "red", "blue", "green", "gold3", "darkorchid" /)
  res@xyLineColors = colors
  zeros = new((/ntiles/), "integer")
  zeros = 0
  res@xyDashPatterns = zeros
  res@xyLineThicknessF = 2.0
  res@gsnYRefLine = 0.0 
  res@gsnYRefLineDashPatterns = 16
  res@tiXAxisFontHeightF = 0.012
  res@tiYAxisFontHeightF = 0.012
  res@tiYAxisString = "Data"
  res@tiXAxisString = "Time"

  ;--- resource for legend ---
  res_text = True
  res_text@txFontColor = "black"
  res_text@txFontHeightF = 0.01
  res_text@txJust = "CenterLeft"

  ;--- plot time-series to different panel ---
  do i = 0, ntiles-1  
    res@vpXF = 0.1+res@vpWidthF*tofloat(i/3)
    res@vpYF = 0.9-res@vpHeightF*tofloat(i%3)
    if (i/3 .eq. 1) then
      res@vpXF = res@vpXF+0.1
    end if
    if (res@vpYF .ne. 0.9) then
      res@vpYF = res@vpYF-tofloat(i%3)*0.1
    end if
    print(res@vpXF+" "+res@vpYF)

    ;--- plot data ---
    plot = gsn_csm_xy(wks, ispan(0, ntime-1, 1), data(i,:), res)
    draw(plot)

    ;--- add legend ---
    gsn_text_ndc(wks, "Tile "+sprinti("%d", i+1), res@vpXF+0.02, res@vpYF-0.02, res_text)

  end do
  
  frame(wks)
end