;----------------------------------------------------------------------
; Generate some dummy data given a lat/lon region of interest, and
; the data min/max.
;----------------------------------------------------------------------
function generate_dummy_data(mindata,maxdata,minlat,maxlat,minlon,maxlon)
local nlat, nlon, lat, lon
begin
  nlat = 30
  nlon = 30
  data = generate_2d_array(10, 12, mindata, maxdata, 0, (/nlat,nlon/))
  lat = fspan(minlat,maxlat,nlat)
  lon = fspan(minlon,maxlon,nlon)
  lat!0 = "lat"
  lon!0 = "lon"
  lat@units = "degrees_north"
  lon@units = "degrees_east"
  data!0 = "lat"
  data!1 = "lon"
  data&lat = lat
  data&lon = lon
  return(data)
end

;----------------------------------------------------------------------
; This procedure adds the latv, lonv polygon to the plot. 
;----------------------------------------------------------------------
procedure add_polygon_to_plot(wks,plot,latv,lonv)
local lnres
begin
  lnres = True
  lnres@gsLineThicknessF   = 5
  lnres@gsLineDashPattern  = 11
  plot@polygon = gsn_add_polyline(wks,plot,lonv,latv,lnres)
end

;----------------------------------------------------------------------
; This procedure adds the location of each lat, lon location to the
; given plot as a filled dot in one of two colors, depending on whether
; it is inside or outside the polygon.
;----------------------------------------------------------------------
procedure add_markers_to_plot(wks,plot,lat1d,lon1d,idx_in,idx_out)
local mkres
begin
  mkres = True
  mkres@gsMarkerIndex = 16
  mkres@gsMarkerSizeF = 5

;---Draw the "inside" points.
  mkres@gsMarkerColor = "Black"
  plot@markers_in = gsn_add_polymarker(wks,plot,lon1d(idx_in),lat1d(idx_in),mkres)

;---Draw the "outside" points.
  mkres@gsMarkerColor = "ForestGreen"
  plot@markers_out = gsn_add_polymarker(wks,plot,lon1d(idx_out),lat1d(idx_out),mkres)
end

;----------------------------------------------------------------------
; Create a filled contour plot and return it (don't draw it).
;----------------------------------------------------------------------
function create_plot(wks,data,mindata,maxdata,minlat,maxlat,minlon,\
                     maxlon,title)
local res
begin
  res                       = True

  res@gsnMaximize           = True
  res@gsnDraw               = False       ; Don't draw plot
  res@gsnFrame              = False       ; Don't advance frame

  res@cnFillOn              = True        ; Turn on contour fill
  res@cnFillPalette         = "BlRe"
  res@cnFillOpacityF        = 0.9         ; Make slightly transparent
  res@cnLinesOn             = False       ; Turn off contour lines

  res@cnLevelSelectionMode  = "ManualLevels"
  res@cnMinLevelValF        = mindata 
  res@cnMaxLevelValF        = maxdata 
  res@cnLevelSpacingF       = 0.01

  res@gsnAddCyclic          = False

  res@mpOutlineBoundarySets = "AllBoundaries"
  res@mpDataSetName         = "Earth..3"
  res@mpOutlineOn           = True
  res@mpMinLatF             = minlat
  res@mpMaxLatF             = maxlat
  res@mpMinLonF             = minlon
  res@mpMaxLonF             = maxlon

  res@pmTickMarkDisplayMode = "Always"

  res@tiMainString          = title

;---Create the plot and return it
  plot = gsn_csm_contour_map(wks,data,res)

  return(plot)
end

;----------------------------------------------------------------------
; Given indexes into 1D arrays, write the indexed locations to
; an ASCII file. Use write_table for formatted output.
;
; Also print to screen using same formatted strings.
;----------------------------------------------------------------------
procedure  write_to_ascii_file(filename,lat,lon,data,idx)
local alist, header
begin
  header = "   lat     lon    data"
  alist = [/lat(idx),lon(idx),data(idx)/]
  write_table(filename, "w", [/header/], "%s")
  write_table(filename, "a", alist, "%7.2f %7.2f %7.3f")

;---Print same data to the screen
  print_table([/header/], "%s")
  print_table(alist, "%7.2f %7.2f %7.3f")
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
  minlat = 26
  maxlat = 32
  minlon = 76
  maxlon = 85
  mindata = -0.04
  maxdata = 0.04
  data = generate_dummy_data(mindata,maxdata,minlat,maxlat,minlon,maxlon)

;---Polygon of interest
  lonv = (/77.98,77.96,79.01,80.02,80.01,79.64,79.63,79.01,77.98/)
  latv = (/31.00,29.70,29.15,28.82,30.05,30.01,30.61,30.64,31.00/)

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

;---Create plot, add latv,lonv polygon, and draw
  plot = create_plot(wks,data,mindata,maxdata,minlat,maxlat,minlon,maxlon,"Original data")
  add_polygon_to_plot(wks,plot,latv,lonv)
  draw(plot)
  frame(wks)

;---Create a lat,lon pair of each original lat,lon point in our dataset
  lat1d = ndtooned(conform(data,data&lat,0))
  lon1d = ndtooned(conform(data,data&lon,1))

;---Check whether each lat,lon pair is in/outside of latv,lonv polygon
  inout   = gc_inout(lat1d,lon1d,latv,lonv)    
  ind_in  = ind(inout)                         ; Get all indexes where inout is True (inside the polygon)
  ind_out = ind(.not.inout)                    ; Get all indexes where inout is False (outside the polygon)

;---Add lat/lon markers to plot and draw again
  add_markers_to_plot(wks,plot,lat1d,lon1d,ind_in,ind_out)
  draw(plot)  
  frame(wks)

;---Mask original data based on inout array
  data1d            = ndtooned(data)              ; Reshape to 1D array so we can apply inout mask
  data1d@_FillValue = -9999
  data1d = where(inout,data1d,data1d@_FillValue)

  data_mask = reshape(data1d,dimsizes(data))      ; Convert back to 2D array and copy metadata
  copy_VarMeta(data,data_mask)                    ; from original variable
  
;---Draw plot of masked data
  plot_mask = create_plot(wks,data_mask,mindata,maxdata,minlat,maxlat,minlon,maxlon,"Masked data")
  add_polygon_to_plot(wks,plot_mask,latv,lonv)
  add_markers_to_plot(wks,plot_mask,lat1d,lon1d,ind_in,ind_out)
  draw(plot_mask)  
  frame(wks)

;---Write to a file
  write_to_ascii_file("file.txt",lat1d,lon1d,data1d,ind_in)

end