;----------------------------------------------------------------------
; This script reads surface temperature off a NetCDF file,
; averages the data based on climate divisions, and plots filled
; polygons based on the average value.
;
; The climate divisions shapefile was downloaded from NOAA:
;
;     http://www.esrl.noaa.gov/psd/data/usclimdivs/boundaries.html
;----------------------------------------------------------------------
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 "./shapefile_mask_data.ncl"
load "./shapefile_average_data.ncl"

;----------------------------------------------------------------------
; Adds dots to the given plot at data's lat/lon grid locations.
;----------------------------------------------------------------------
undef("add_latlon_points")
procedure add_latlon_points(wks,plot,data)
local mkres
begin
  mkres                 = True
  mkres@gsMarkerIndex   = 16         ; Filled dots
  mkres@gsMarkerSizeF   = 0.002      ; Make them smaller
  mkres@gsMarkerColor   = "black"
  mkres@gsnCoordsAttach = True
  gsn_coordinates(wks,plot,data,mkres)
end

;----------------------------------------------------------------------
; Add shapefile outlines to the given plot.
;----------------------------------------------------------------------
undef("add_shapefile_outlines")
procedure add_shapefile_outlines(wks,plot,shp_filename)
local lnres
begin
;---Resources for polyline
  lnres                  = True
  lnres@gsLineColor      = "Navyblue"
; lnres@gsLineThicknessF = 3.0            ; 3x thickness

  plot@lines = gsn_add_shapefile_polylines(wks, plot, shp_filename, lnres)
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
;---Estimated map area to focus on (United States)
  minlat = 22
  maxlat = 51
  minlon = -126
  maxlon = -65

;---Read data off NetCDF file
  fname    = "ts_Amon_CESM1-CAM5_historical_r1i1p1_185001-200512.nc"
  f        = addfile(fname,"r")
  ts       = f->ts(0,:,:)        ; Original ts: 1872 x 192 x 288   (time x lat x lon)
  ts       = ts-273.15           ; Convert from Kelvin -> Celsius.
  ts@units = "degC"

  ts = lonFlip(ts)      ; 0->360 to -180->180
  printVarSummary(ts)

;
; Use shapefile to mask the dummy array based on shapefile outline
; This will make the data averaging go a little faster.
;
  shp_filename = "GIS.OFFICIAL_CLIM_DIVISIONS.shp"
  ts_mask      = shapefile_mask_data(ts({minlat:maxlat},{minlon:maxlon}),shp_filename,True)

;---Start the graphics
  wks  = gsn_open_wks("x11","shapefiles_data_recti")

  res                     = True

  res@gsnMaximize         = True           ; maximize plot in frame
  res@gsnDraw             = False          ; don't draw plot yet
  res@gsnFrame            = False          ; don't advance frame yet

  res@mpDataBaseVersion   = "MediumRes"    ; slightly better resolution
  res@mpFillOn            = False
  res@mpOutlineOn         = False

;---Zoom in on area of interest
  res@mpMinLatF           = minlat
  res@mpMaxLatF           = maxlat
  res@mpMinLonF           = minlon
  res@mpMaxLonF           = maxlon
;
; Create a map plot so we can add color-coded climate (by average) 
; division polygons later.
;
  res@gsnLeftString      = "Average " + ts@long_name
  res@gsnRightString     = ts@units

  plot_map = gsn_csm_map(wks,res)    

;---Create filled contour plot of original data, for comparison.
  cnres                     = res
  cnres@cnFillOn            = True
  cnres@cnLinesOn           = False
  cnres@cnLineLabelsOn      = False
  cnres@cnLevelSpacingF     = 2.5
  cnres@lbLabelBarOn        = False
  cnres@gsnAddCyclic        = False
  cnres@gsnLeftString      = ts@long_name
  cnres@gsnRightString     = ts@units

  plot_data = gsn_csm_contour_map(wks,ts({minlat:maxlat},{minlon:maxlon}),cnres)

;---Get the colors and levels used for contour plot
  getvalues plot_data@contour
    "cnLevels"     : levels
    "cnFillColors" : colors
  end getvalues

;---Add color-coded polygons that represent average of data in a particular climate division
  plot_averages_using_shapefile_features(wks,plot_map,ts,shp_filename,levels,colors)

;---Add shapefile outlines to both plots
  add_shapefile_outlines(wks,plot_data,shp_filename)
  add_shapefile_outlines(wks,plot_map,shp_filename)

;---Add dots at lat/lon locations
;  add_latlon_points(wks,plot_data,ts)
  add_latlon_points(wks,plot_map,ts_mask)

;
; Panel both plots. All the stuff we've added to both plots,
; shapefile outlines, lat/lon dots, filled polygons, will
; be seen here.
;
  pres                  = True
  pres@gsnMaximize      = True
  pres@gsnPanelLabelBar = True
  gsn_panel(wks,(/plot_data,plot_map/),(/2,1/),pres)

end