;----------------------------------------------------------------------
; 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
;---Read var1 
         a = addfile("/home/kunal/DAY_30.nc","r")
  filename = "/home/kunal/IND_adm/india_state.shp"
   
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"


;---Draw the "outside" points.
  mkres@gsMarkerColor = "ForestGreen"
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
 ; maxdata = 1500
  data = a->bcfire(0,:,:)

  
  
  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/)

  lonb = (/80.10,80.98,81.57,81.98,83.01,84.00,83.99,83.01,83.01,80.98,81.00,80.80,80.10 /)
 
  latb = (/28.91,28.43,28.01,27.98,27.47,27.38,28.64,28.99,29.40,29.56,29.99,30.00,28.91/)

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

;---Create a lat,lon pair of each original lat,lon point in our dataset
  lat1d = ndtooned(conform(data,data&latitude,0))
  lon1d = ndtooned(conform(data,data&longitude,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)

  inoutb   = gc_inout(lat1d,lon1d,latb,lonb) 
  indb_in  = ind(inoutb)                         ; Get all indexes where inout is True (inside the polygon)
  indb_out = ind(.not.inoutb)                    ; Get all indexes where inout is False (outside the polygon)


;---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
  
  datab1d            = ndtooned(data)              ; Reshape to 1D array so we can apply inout mask
  datab1d@_FillValue = -9999
  datab1d = where(inoutb,datab1d,datab1d@_FillValue)

  datab_mask = reshape(datab1d,dimsizes(data))      ; Convert back to 2D array and copy metadata
  copy_VarMeta(data,datab_mask)    


;---Write to a file
  write_to_ascii_file("2004_DAY_30_UK_bc.txt",lat1d,lon1d,data1d,ind_in)
  write_to_ascii_file("2004_DAY_30_NEP_bc.txt",lat1d,lon1d,datab1d,indb_in)

end