sub.draw_series.ncl

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; sub.draw_series.ncl
; Carl Schreck (cjschrec@ncsu.edu)
; October 2014
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Draw time series plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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 "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
load "$CJS_NCL_LIB/print_clock.ncl"
load "$CJS_NCL_LIB/lib.array.ncl"
load "$CJS_NCL_LIB/lib.cjs_graphics.ncl"


undef ( "draw_series" )
function draw_series( \
  io_wks    [1] : graphic, \
  i_prefix  [1] : string, \
  i_sensor  [*] : string, \
  i_seriesType [1] : string, \
  i_minTime [1] : numeric, \
  i_maxTime [1] : numeric, \
  i_res     [1] : logical \
)
begin ; draw_series
  ; These are some parameters that could be useful to have up top
  if( i_prefix.eq."olr" ) then
    buffer = 0
  else
    buffer = 365
  end if
  shortSmooth = 365
  longSmooth = shortSmooth*5
  inPath = "/home/carl/data/olr/compare/time_series/" + i_prefix + "."  \
      + i_sensor + ".20S-20N.nc"

  res = i_res
  res@gsnLeftString = i_prefix
  res@trXMinF = i_minTime
  res@trXMaxF = i_maxTime
  res@xyExplicitLegendLabels = i_sensor
  res@tiXAxisString = ""
  res@xyMarkLineMode = "Lines"
  res@xyMonoLineThickness = True
  res@vpHeightF = 0.4

  print_clock( "Reading times" )
  do sensorCounter = 0, dimsizes(i_sensor)-1
    inFile = addfile( inPath(sensorCounter), "r" )
    inTime := inFile->time
    inTime = floor(inTime)
    if( sensorCounter.eq.0 ) then
      allTimes = inTime
    else
      newAllTimes := array_append_record( allTimes, inTime, 0 )
      allTimes := newAllTimes
    end if
  end do
  xData = remove_duplicates( allTimes )
  xData@units = inTime@units
  xData!0 = "time"
  xData&time = xData

  print_clock("Reading data")
  yData = new( (/ dimsizes(i_sensor), dimsizes(xData) /), float )
  yData!1 = "time"
  yData&time = xData

  do sensorCounter = 0, dimsizes(i_sensor)-1
    inFile = addfile( inPath(sensorCounter), "r" )
    inTime := floor(inFile->time)
    startInd = buffer
    endInd = dimsizes(inTime)-1-buffer
    inData := inFile->$i_seriesType$(startInd:endInd)
    yData(sensorCounter,{inTime(startInd):inTime(endInd)}) = inData
  end do

  print_clock("Trending")
  trendData = new( dimsizes(yData), double )
  do sensorCounter = 0, dimsizes(i_sensor)-1
    rc = regline( xData({i_minTime:i_maxTime}),  \
        yData(sensorCounter,{i_minTime:i_maxTime}) )
    trendData(sensorCounter,:) = rc * xData + rc@yintercept
  end do

  print_clock("Smoothing")
  yData = runave_n_Wrap( yData, shortSmooth, 0, 1 )
  smoothData = runave_n_Wrap( yData, longSmooth, 0, 1 )

  print_clock("Plotting")
  resTick = True
  resTick@ttmFormat      = "%Y"
  resTick@ttmAxis        = "XB"
  month_axis_labels( i_minTime, 60, 12, 1, res, resTick )

  res@xyLineThicknessF = 2
  res@xyBoldColors = False
  rawPlot = cjs_draw_xy( io_wks, xData, yData, res )

  res@xyLineThicknessF = 1
  res@xyBoldColors = False
  trendPlot = cjs_draw_xy( io_wks, xData, trendData, res )
  overlay( rawPlot, trendPlot )

  res@xyLineThicknessF = 6
  res@xyBoldColors = True
  smoothPlot = cjs_draw_xy( io_wks, xData, smoothData, res )
  overlay( rawPlot, smoothPlot )

  retVal = rawPlot

  return(retVal)
end; draw_series