;----------------------------------
; CALIPSO Satellite TRACK OVERLAYS for the month of January
;----------------------------------

load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

;----------------------------------
; SET DIRECTORIES
;----------------------------------

  diri = "/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/2007_to_2017_January_PalmerLTERDATA_for_Overlay_and_VFM/2017_January_PalmerLTERDATA_Overlay_Data/" ; input directory
  file_name = systemfunc("cd "+diri+" ; ls CAL*.hdf")
  printVarSummary(file_name)
  
  fig_diri = "/Users/srishtidasarathy/Documents/Bowman/CALIPSO_Data_Processing/2007_to_2017_January_PalmerLTER_Overlay_Figures/" ;directory that image gets saved in
  fig_save = fig_diri + "CALIPSO_Overlay_2017_January"; will save this image in the directory specified (fig_diri)

;----------------------------------
; OPEN WORKSTATION
;----------------------------------

  wks = gsn_open_wks("png", fig_save) ; open workstation as the first step. This doesn't need to be looped through
  n = dimsizes(file_name) ;calls out the number of files

;----------------------------------
; PLOT RESOURCES
;----------------------------------

   mpres = True
   mpres@gsnDraw = False; False so that it doesn't create a new image everytime and name it something like ".00002" etc.
   mpres@gsnFrame = False;
   
   ;Use lat/lon coordinates to limit area. I am trying to zoom in on WAP
   mpres@mpLimitMode = "LatLon"
   mpres@mpMinLatF  = -71
   mpres@mpMinLonF  = -80
   mpres@mpMaxLatF = -58
   mpres@mpMaxLonF = -58
   
   mpres@tiMainString    = "CALIPSO Tracks for January 2017, Western Antarctic Peninsula"   ; add titles
   mpres@tiXAxisString   = "Longitude"  
   mpres@tiYAxisString   = "Latitude"
   
   mpres@tiMainJust           = "Center"
   mpres@tiMainPosition       = "Center"
   mpres@tiMainFontAspectF = 1.8
   mpres@tiMainFontHeightF = 0.02
   mpres@tiDeltaF = 1.8

   mpres@tiXAxisJust          = "CenterRight"
   mpres@tiXAxisPosition      = "Center"

   mpres@tiYAxisJust          = "Center"
   mpres@tiYAxisSide          = "Left"
   mpres@tiYAxisPosition      = "Center"
   
   mpres@tiYAxisFontHeightF = 0.0180
   mpres@tiYAxisOffsetXF    = -0.06
   mpres@tiXAxisFontHeightF = 0.018
   mpres@tiXAxisOffsetYF    = -0.045
   mpres@tiXAxisOffsetXF    = 0.05
     
  map = gsn_csm_map(wks,mpres) ; outside of the loop because I want one map object
  lnid = new(n,"graphic") ; n is dimension size of number of files in the directory. assigns these to graphical object.
  
;----------------------------------
; LOOP START
;----------------------------------

  do i = 0, n-1
    print(i)
    hdf4_file = addfile(diri+file_name(i), "r") 
    dat = str_get_cols(file_name(i),31,40)
  
    plres = True
    plres@gsLineColor = "red"
   
    latitude = hdf4_file->Latitude
    longitude = hdf4_file->Longitude
  
    lnid(i) = gsn_add_polyline(wks,map,longitude(:,0),latitude(:,0),plres) ; only the tracks need to be looped through
 
    ;With each iteration of the loop the number of latitude and longitude values changes. Must have a delete item here.
    delete(latitude)  
    delete(longitude)
;------------------------   
   end do
   draw(map)
   frame(wks)
   
end