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"

;----------------------------------------------------------------------
; Procedure to add shapefile outlines to the given plot.
;----------------------------------------------------------------------
procedure add_shp_outlines(wks,plot,shp_filename)
local lnres
begin
;---Resources for polyline
  lnres                  = True
  lnres@gsLineColor      = "Red"
  lnres@gsLineThicknessF = 3.0            ; 3x thickness

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

;----------------------------------------------------------------------
; Function to create ummy rectilinear data over a map area of interest.
;----------------------------------------------------------------------
function create_dummy_array(minlat,maxlat,minlon,maxlon,nlon,nlat)
local lat1d, lon1d
begin
;---Generate some dummy data over map
  data            = generate_2d_array(10, 10, 0., 100., 0, (/nlat,nlon/))
  data@_FillValue = -9999

;---Add lat/lon coordinate array information.
  lat1d       = fspan(minlat,maxlat,nlat)
  lon1d       = fspan(minlon,maxlon,nlon)
  lat1d@units = "degrees_north"
  lon1d@units = "degrees_east"
  data!0      = "lat"
  data!1      = "lon"
  data&lat    = lat1d
  data&lon    = lon1d

  return(data)
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------

begin
;---Create dummy array
  GRID_TYPE = "finer"     ; "coarser"

  minlat = 20
  maxlat = 60
  minlon = -125
  maxlon = -60
  if(GRID_TYPE.eq."finer") then
    nlat = 64
    nlon = 128
  else
    nlat = 32
    nlon = 64
  end if
  data = create_dummy_array(minlat,maxlat,minlon,maxlon,nlat,nlon)

;---Start the graphics
  wks  = gsn_open_wks("png",GRID_TYPE)

  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@cnFillOn            = True
  res@cnFillPalette       = "MPL_Set3"
; res@cnFillOpacityF      = 0.5            ; Fade the plot slightly
  res@cnLinesOn           = False
  res@cnLineLabelsOn      = False

  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = ispan(6,96,4)
  res@lbLabelBarOn         = False

  res@gsnAddCyclic        = False

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

;---Zoom in on area of interest
  res@mpMinLatF           = minlat
  res@mpMaxLatF           = maxlat
  res@mpMinLonF           = minlon
  res@mpMaxLonF           = maxlon

  res@tiMainString        = "Original data"
  plot_orig = gsn_csm_contour_map(wks,data,res)        ; Create plot of original data

;---Use shapefile to mask this dummy array based on shapefile outline
  dir          = ncargpath("data") + "/shp/"
  shp_filename = dir + "mrb.shp"
  data_mask    = shapefile_mask_data(data,shp_filename,True)


  res@tiMainString = "Data with mask applied (" + GRID_TYPE + " grid)"
  plot_mask = gsn_csm_contour_map(wks,data_mask,res)   ; Create plot of masked data

;---Add shapefile outlines to both plots
  add_shp_outlines(wks,plot_orig,shp_filename)
  add_shp_outlines(wks,plot_mask,shp_filename)

;---Add dots at the lat/lon grid locations
  mkres                 = True
  mkres@gsMarkerIndex   = 16         ; Filled dots
  mkres@gsMarkerSizeF   = 0.001      ; Make them small
  mkres@gsnCoordsAttach = True
  gsn_coordinates(wks,plot_orig,data,mkres)
  gsn_coordinates(wks,plot_mask,data_mask,mkres)

  pres                  = True
  pres@gsnMaximize      = True
  pres@gsnPanelLabelBar = True
  gsn_panel(wks,(/plot_orig,plot_mask/),(/2,1/),pres)
end