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("/home/kunal/AFS_201504_21_to_30.nc","r")
  data = a->SWGNTCLR(0,:,:)   ; lat x lon

;---Region of interest
  minlat = 26
  maxlat = 32
  minlon = 76
  maxlon = 85

;---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 = (/"Uttaranchal"/)

  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 data after masking
  res@tiMainString      = "Plot"
  plot_mask     = gsn_csm_contour_map(wks,data_mask,res)

;---Add shapefile outlines to both plots
  lnres = True

  id1 = gsn_add_shapefile_polylines(wks,plot_mask,shp_filename_ind,lnres)
  id2 = 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_mask/),(/1,1/),pres)


end