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/wrf/WRFUserARW.ncl"
load "./shp_utils.ncl"

begin
;---Open WRF output file and read data
  filename     = "wrfout_d01_2002-01.nc"
  a            = addfile(filename,"r")
  ter          = a->HGT(0,:,:)             ; Read the variable to memory
  lat2d        = a->XLAT(0,:,:)
  lon2d        = a->XLONG(0,:,:)
  ter@lat2d    = lat2d
  ter@lon2d    = lon2d

;---Contour levels and colors to use
  color_map = "OceanLakeLandSnow"
  levels  = ispan(1,1100,20)
  nlevels = dimsizes(levels)
  colors  = span_color_rgba(color_map,nlevels+1)

  wks = gsn_open_wks("png","shapefiles")
  gsn_define_colormap(wks,color_map)  ; Change color map

  res                   = True               ; Contour options
  res = wrf_map_resources(a,res)             ; Add WRF map resources
  res@gsnDraw           = False
  res@gsnFrame          = False
  res@mpOutlineOn       = False  ; Turn off map outlines
  res@mpFillOn          = False  ; Turn off map fill
  res@cnFillOn          = True
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels          = levels
  res@cnFillPalette     = colors
  res@lbLabelBarOn      = False
  res@cnInfoLabelOn     = False
  res@cnLinesOn         = False
  res@cnLineLabelsOn    = False

  plot1 = gsn_csm_contour_map(wks,ter,res)       ; Create contour plot

;---Attach the shapefile polylines using files read off gadm.org/country.
  lnres                  = True
;  lnres@gsLineColor      = "gray"
  lnres@gsLineThicknessF = 0.5

  usa_shp_name = "USA_adm/USA_adm2.shp"
  usa_id1 = gsn_add_shapefile_polylines(wks,plot1,usa_shp_name,lnres)

  draw(plot1)      ; This draws the contour/map and shapefile outlines.
  frame(wks)       ; Advance the frame

;---Create just a map plot for the the filled polygons to be added to.
  mpres          = True
  mpres          = wrf_map_resources(a,mpres)
  mpres@gsnDraw  = False
  mpres@gsnFrame = False

  plot2 = gsn_csm_map(wks,mpres) 

;
; Calculate data averages over each county of interest and also
; attach colored polygons to existing map based on average.
;
  ter_avg = plot_averages_by_county(wks,plot2,ter,lat2d,lon2d,\
                                    usa_shp_name,levels,colors)
  printMinMax(ter,0)
  printMinMax(ter_avg,0)

  usa_id2 = gsn_add_shapefile_polylines(wks,plot2,usa_shp_name,lnres)

  draw(plot2)
  frame(wks)

  pres = True
  pres@gsnMaximize = True
  gsn_panel(wks,(/plot1,plot2/),(/1,2/),pres)
end