load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "./shapefile_utils.ncl"

begin
;---Read data
  a    = addfile("dataforfigure2.nc","r")
  data = a->MYD08_D3_6_Deep_Blue_Aerosol_Optical_Depth_550_Land_Mean(0,:,:)   ; lat x lon

;---Region of interest
  minlat = 23
  maxlat = 34
  minlon = 72
  maxlon = 89

;---These are the two shapefiles we'll use for masking.
  shp_filename_ind = "IND_adm/IND_adm1.shp"
  shp_filename_npl = "NPL_adm/NPL_adm0.shp"

;---Mask the data based on Nepal outline in shapefile
  opt             = True
  opt@minlat      = minlat    ; this makes the masking go faster
  opt@maxlat      = maxlat
  opt@minlon      = minlon
  opt@maxlon      = maxlon
  opt@debug       = True
  data_mask_npl = shapefile_mask_data(data,shp_filename_npl,opt)

;---Mask the data based on Northern India outlines in shapefile
  opt@shape_var   = "NAME_1"
  opt@shape_names = (/"Uttar Pradesh","Himachal Pradesh","Uttaranchal","Punjab","Haryana"/)

  data_mask = shapefile_mask_data(data,shp_filename_ind,opt)

;---Combine the two sets of data back into "data_mask"
  data_mask = where(.not.ismissing(data_mask_npl),data_mask_npl,data_mask)

;---Debug prints
  printMinMax(data,0)
  printMinMax(data_mask,0)

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

  res                   = True

  res@gsnDraw           = False    ; will panel later
  res@gsnFrame          = False

  res@gsnAddCyclic      = False    ; turn off longitude cyclic point

  res@cnFillOn          = True     ; turn on contour fill
  res@cnFillPalette     = "WhBlGrYeRe"
  res@cnLinesOn         = False    ; turn off contour lines
  res@cnLineLabelsOn    = False    ; turn off line labels
  res@cnFillMode        = "RasterFill"
  res@cnRasterSmoothingOn = True
  res@lbLabelBarOn      = False    ; will add labelbar in panel

;---Fix the contor levels for all plots
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/0.3,0.5,0.7,0.9,1.,2./)
 
;---Zoom in on region of interest
  res@mpMinLatF         = minlat
  res@mpMaxLatF         = maxlat
  res@mpMinLonF         = minlon
  res@mpMaxLonF         = maxlon

  res@mpFillOn          = False     ; will add shapefile outlines later
  res@mpOutlineOn       = False

;---Titles and tickmarks
  res@gsnLeftString     = ""
  res@gsnRightString    = ""
  res@tiMainString      = "original data"
  res@pmTitleZone       = 1
  res@pmTickMarkDisplayMode = "Always"   ; better tickmarks

;---Plot of original data before masking
  plot_orig = gsn_csm_contour_map(wks,data,res)

;---Plot of data after masking
  res@tiMainString      = "masked data"
  plot_mask     = gsn_csm_contour_map(wks,data_mask,res)

;---Add shapefile outlines to both plots
  lnres = True
;  lnres@gsLineThicknessF = 2.0
  id1 = gsn_add_shapefile_polylines(wks,plot_orig,shp_filename_ind,lnres)
  id4 = gsn_add_shapefile_polylines(wks,plot_mask,shp_filename_ind,lnres)

  id2 = gsn_add_shapefile_polylines(wks,plot_orig,shp_filename_npl,lnres)
  id3 = gsn_add_shapefile_polylines(wks,plot_mask,shp_filename_npl,lnres)

;---Panel original data and masked data
  pres = True
  pres@gsnMaximize = True
  pres@gsnPanelLabelBar = True
  gsn_panel(wks,(/plot_orig,plot_mask/),(/2,1/),pres)

;---Add filled dots at grid locations and panel
  mkres = True
  mkres@gsMarkerSizeF = 3
  mkres@gsMarkerIndex = 16
  mkres@gsnCoordsAttach = True
  gsn_coordinates(wks,plot_orig,data,mkres)
  gsn_coordinates(wks,plot_mask,data_mask,mkres)
  gsn_panel(wks,(/plot_orig,plot_mask/),(/2,1/),pres)
end