;***************************************************************
; binning_2.ncl
;
; Concepts illustrated:
;   - Reading multiple fileslwith data from MODIS "MCD19A2"
;   - Binning and summing the data from each pass [swath]
;   - Averaging binned data
;   - Plotting
;   - Creating a netCDF file
;
;***************************************************************
; Utility functions to get (1) date information (2) facilitate netCDF 
;********************************************************************
; MCD19A2.AYYYYDDD.hHHvVV.CCC.YYYYDDDHHMMSS.hdf
; YYYYDDD = Year and Day of Year of acquisition
; hHH = Horizontal tile number (0-35)
; vVV = Vertical tile number (0-17)
; CCC = Collection number
; YYYYDDDHHMMSS = Production Date and Time
; Example:
; MCD19A2.A2019161.h32v11.006.2019163201550.hdf
;********************************************************************

undef("parseFileName")
function parseFileName ( fNam[*]:string )
; parse the file name(s) and extract data information
; https://www.ncl.ucar.edu/Document/Functions/Built-in/str_get_cols.shtml
; 012345678901234567890123456789012345678901234
;           1         2         3         4
; MCD19A2.AYYYYDDD.hHHvVV.CCC.YYYYDDDHHMMSS.hdf
;          YYYYDDD  tile
local n, nNam, time_info, YYYY, DDD, yyyyddd, yyyymmdd, mmdd, mm, dd  
begin
   nNam = dimsizes( fNam )

   YYYY      = toint(str_get_cols(fNam, 9,12))     ; YYYY 
   DDD       = toint(str_get_cols(fNam,13,15))     ; DDD

   yyyyddd   = YYYY*1000 + DDD
   yyyymmdd  = yyyyddd_to_yyyymmdd(yyyyddd) 
   mmdd      = yyyymmdd-(yyyymmdd/10000)*10000
   mm        = mmdd/100                            
   dd        = mmdd-mm*100

   time_info    = new( (/nNam,6/), integer,"No_FillValue")

   do n=0,nNam-1
      time_info(n,0) = yyyyddd 
      time_info(n,1) = YYYY      
      time_info(n,2) = DDD
      time_info(n,3) = yyyymmdd
      time_info(n,4) = mm           
      time_info(n,5) = dd           
   end do
   
   return (time_info)   
end

;**************************************************************
;                     MAIN
;**************************************************************

  vNam  = "Optical_Depth_047_grid1km"
  PLOT  = True

;*****************************************************************
; binned grid definition     
; The grid must be constant spacing in lat and in lon.
; The spacing in the lat and lon directions may be different.
; so: 2x2, 3x5, etc
;*****************************************************************


  nlat  = 10800                     ; grid coordinates
  mlon  = 10800      
;  nlat  = 3600
;  mlon  = 7200
  lat   = latGlobeFo(nlat,"lat","latitude","degrees_north")
  lon   = lonGlobeFo(mlon,"lon","longitude","degrees_east")
  lon   = (/ lon - 180. /)  ; subtract 180 from all values 
  lon&lon = lon             ; update coordinates

;*****************************************************************
; Variables to hold binned quantities
;*****************************************************************

  ; creating a new array to place numerous days/months/years together
;  X01 = new( (/4,nlat,mlon/),"integer")
;  X02 = new((/4,1200,1200/),"integer")
  GBIN  = new ( (/nlat,mlon/), float ) 
  GKNT  = new ( (/nlat,mlon/), integer ) 
  GBIN2  = new ( (/nlat,mlon/), float )
  GKNT2  = new ( (/nlat,mlon/), integer )
  GBIN3  = new ( (/nlat,mlon/), float )
  GKNT3  = new ( (/nlat,mlon/), integer )
;  X01   = new ( (/nlat,mlon/), integer )

  GBIN  = 0.0                      ; initialization
  GKNT  = 0
  GBIN2  = 0.0                      ; initialization
  GKNT2  = 0
  GBIN3  = 0.0                      ; initialization
  GKNT3  = 0

  X01   = 0

;*****************************************************************
; File info                          
;     HDFEOSVersion : HDFEOS_V2.19
;     StructMetadata_0 : GROUP=SwathStructure
;*****************************************************************

  diri = "./"
;  fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2017259*.hdf")
;  fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2017*.hdf")
;************ multiple year averaging ****************************
  fili = systemfunc("cd "+diri+" ; ls MCD19A2.A200*.hdf")
  fili1 = systemfunc("cd "+diri+" ; ls MCD19A2.A201*.hdf")
  nfil = dimsizes( fili )
  nfil1 = dimsizes( fili1 )
  print("nfil="+nfil)
  print("nfil1="+nfil1)
  f0   = addfile(diri+fili(0)+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays
  f1   = addfile(diri+fili1(0)+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays 
  x    = short2flt( f0->$vNam$   )   ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km )
  y    = short2flt( f1->$vNam$   )   ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km 
printVarSummary(x)
printMinMax(x,0)
print("=============================")
  nmsg = num(ismissing(x))
  ngood= num(.not.ismissing(x))
  npts = nmsg+ngood
  npc  = (ngood/tofloat(npts))*100
print("nmsg ="+nmsg)
print("ngood="+ngood)
print("npc (% good)="+npc)
print("=============================")
  nx   = num(x.gt.-100 .and. x.lt.5000)   ; #  of values bwtweem -100 and 5000 after scale fac -.01 and 5
; *****************************************************************
  nx   = num(x.ge. -0.1 .and. x.le.5.0)   ; values after unpacking
                                          ; within allowed range
print("nx="+nx+":  this should match ngood="+ngood) ; ie: there are no out of range values
;  X01  = where(x.ge. -0.1 .and. x.le.5.0, 1, 0)  
printVarSummary(X01)                               
  n1n   = num(X01.eq.1)
  n1s   = sum(X01)
;  X02   = nls
print("n1n="+n1n+"   n1s="+n1s)
print("=============================")
;*****************************************************************
; Loop over the files: Sum the quantities
; *NO* data filtering is performed: No range checking nor use of AOD_QA_grid1km
;*****************************************************************
  tStrt = systemfunc("date")         ; time the loop (wall clock)

  do nf=0,nfil-1
     print(nf+"   "+fili(nf))
     f       = addfile(diri+fili(nf)+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays 
                                     ; read data
     lat2d   = f->GridLat_grid1km 
     lon2d   = f->GridLon_grid1km 
     x       = short2flt( f->$vNam$   )   ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km )
 
  printMinMax(x,0)
     dimx    = dimsizes(x)
     norbit  = dimx(0)

    do no=0,norbit-1                 ; loop over all orbits within the current file  
;    do no=3,3                 ; loop over all orbits within the current file  
        bin_sum(GBIN,GKNT,lon,lat   \
               ,ndtooned(lon2d), ndtooned(lat2d),ndtooned(x(no,:,:)) )
   end do
     delete(  x  )                   ; size may change
;     delete( X01 )                   ; size may change
  end do
 
  wallClockElapseTime(tStrt, "Main Sum Loop: nlat="+nlat+", mlon="+mlon, 0)

;*****************************************************************
; Perform averaging
;*****************************************************************
                                    ; avoid division by 0.0
  GKNT     = where(GKNT.eq.0 , GKNT@_FillValue, GKNT)
  GBIN     = GBIN/GKNT
 
  GBIN!0   = "lat"
  GBIN!1   = "lon"
  GBIN&lat =  lat
  GBIN&lon =  lon

  copy_VarCoords(GBIN, GKNT)        ; copy coords

;  if (isfilevaratt(x, vNam, "long_name")) then
 if (isfilevaratt(f, vNam, "long_name")) then
      GBIN@long_name = "BINNED AVERAGE: "+vNam
      GKNT@long_name = "BINNED COUNT: "+vNam
;      X01@long_name  = "BINNED COUNT: "+vNam
  end if

 if (isfilevaratt(f, vNam, "units")) then
     GBIN@units     = f->$vNam$@units
  end if
printVarSummary(GKNT)
;************************************************************************
; Process the second set of data on map
;************************************************************************
;*****************************************************************
; Loop over the files: Sum the quantities
; *NO* data filtering is performed: No range checking nor use of AOD_QA_grid1km
;*****************************************************************
  tStrt = systemfunc("date")         ; time the loop (wall clock)

  do nf1=0,nfil1-1
     print(nf1+"   "+fili1(nf1))
     f1       = addfile(diri+fili1(nf1)+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays 
                                     ; read data
     lat2d   = f1->GridLat_grid1km 
     lon2d   = f1->GridLon_grid1km 
     y       = short2flt( f1->$vNam$   )   ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km )
 
  printMinMax(y,0)
     dimy    = dimsizes(y)
     norbit  = dimy(0)

    do no=0,norbit-1                 ; loop over all orbits within the current file  
;    do no=3,3                 ; loop over all orbits within the current file  
        bin_sum(GBIN2,GKNT2,lon,lat   \
               ,ndtooned(lon2d), ndtooned(lat2d),ndtooned(y(no,:,:)) )
   end do
     delete(  y  )                   ; size may change
;     delete( X01 )                   ; size may change
  end do
 
  wallClockElapseTime(tStrt, "Main Sum Loop: nlat="+nlat+", mlon="+mlon, 0)

;*****************************************************************
; Perform averaging
;*****************************************************************
                                    ; avoid division by 0.0
  GKNT2     = where(GKNT2.eq.0 , GKNT2@_FillValue, GKNT2)
  GBIN2     = GBIN2/GKNT2
 
  GBIN2!0   = "lat"
  GBIN2!1   = "lon"
  GBIN2&lat =  lat
  GBIN2&lon =  lon

  copy_VarCoords(GBIN2, GKNT2)        ; copy coords

;  if (isfilevaratt(x, vNam, "long_name")) then
 if (isfilevaratt(f1, vNam, "long_name")) then
      GBIN2@long_name = "BINNED AVERAGE: "+vNam
      GKNT2@long_name = "BINNED COUNT: "+vNam
;      X01@long_name  = "BINNED COUNT: "+vNam
  end if

 if (isfilevaratt(f1, vNam, "units")) then
     GBIN2@units     = f1->$vNam$@units
  end if
printVarSummary(GKNT2)
; % difference
;   GBIN3 = ((GBIN2-GBIN)/GBIN)*100
;  RATIO
   GBIN3= GBIN2/GBIN
   GKNT3 = GKNT2/GKNT

;*****************************************************************
; time/date
;*****************************************************************
  dateinfo= parseFileName(fili(0))

  yyyyddd = dateinfo(0,0)    ; this is the date info for the 1st file
  yyyy    = dateinfo(0,1)
  ddd     = dateinfo(0,2)
  yyyymmdd= dateinfo(0,3)
  month   = dateinfo(0,4)
  day     = dateinfo(0,5)
;*****************************************************************
; PLOTS   
;*****************************************************************

  if (PLOT) then

      setvalues NhlGetWorkspaceObjectId()
                "wsMaximumSize" : 1000000000
      end setvalues   

      plot    = new ( 2, "graphic")

     
      wks     = gsn_open_wks("png", "MCD19A2_Ratio-1x1km_-20to20"); send graphics to PNG file
      
      res                     = True          ; Use plot options
      res@gsnDraw             = False
      res@gsnFrame            = False

      res@cnFillOn            = True          ; Turn on color fill
;     res@cnFillPalette       = "wh-bl-gr-ye-re"  ; set color map
      res@cnFillPalette       = "BlAqGrYeOrReVi200"  ; set color map
      res@cnFillMode          = "RasterFill"  ; Turn on raster color
;  res@lbLabelAutoStride=True ; ensure labels do not overlap
;.........................................................................
;  res@cnLevelSelectionMode = "ManualLevels"
;  res@cnMinLevelValF = -50.
;  res@cnMaxLevelValF = 50.
;  res@cnLevelSpacingF = 10.
;   plot(0) = gsn_csm_contour_map(wks,GBIN3, res)
;.........................................................................
; res@cnLevelSelectionMode = "AutoLevels"
;  res@cnMinLevelValF = 0.
;  res@cnMaxLevelValF = 1.
;  res@cnLevelSpacingF = .1
;  plot(1) = gsn_csm_contour_map(wks,GKNT3, res)
;..............................................................................

     res@cnLinesOn           = False         ; Turn off contourn lines
  res@mpDataSetName         = "Earth..4"   ; This database contains
                                           ; divisions for other countries.
  res@mpDataBaseVersion = "MediumRes"

; The following 5 lines are used to create a zoomed image.
; Set limits of map, based on the min/max of the dataset latitude/longitude
;  res@mpLimitMode = "LatLon"
;  res@mpMinLatF        = min(data@lat2d)
;  res@mpMaxLatF        = max(data@lat2d)
;  res@mpMinLonF        = min(data@lon2d)
;  res@mpMaxLonF        = max(data@lon2d)
;***************** Full NYS Plot ***************
  res@mpLimitMode = "LatLon"
  res@mpMinLatF = 40
  res@mpMaxLatF = 46
  res@mpMinLonF = -80.
  res@mpMaxLonF = -71.
;************** NYC Plot Only *****************
;  res@mpLimitMode = "LatLon"
;  res@mpMinLatF = 40
;  res@mpMaxLatF = 41
;  res@mpMinLonF = -74.5
;  res@mpMaxLonF = -73.5

;  res@mpUSStateLineColor = "gray"
  res@mpUSStateLineDashPattern = 7
  res@mpUSStateLineThicknessF = 1
;  res@mpNationalLineThicknessF= 2.         ; double the thickness of national boundaries
  res@mpNationalLineThicknessF= 1.         ; double the thickness of national boundaries
  res@mpOutlineBoundarySets = "National"                   ; draw national boundaries
  res@mpOutlineSpecifiers = (/"Canada : Provinces","United States : States"/)  ; draw US States, Canadian provinces
;  res@mpOutlineSpecifiers = (/"China:states","Mongolia","India","Nepal","Bhutan","Bangladesh","Myanmar", \
;                                    "Thailand","Cambodia","Vietnam","Taiwan","Japan"/)
      res@mpCenterLonF        = 180           ; set map center at 180
      res@lbOrientation       = "vertical"

;.........................................................................
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF = 0.0
  res@cnMaxLevelValF = 1.3
  res@cnLevelSpacingF = .1
   plot(0) = gsn_csm_contour_map(wks,GBIN3, res)
;.........................................................................
 res@cnLevelSelectionMode = "AutoLevels"
;  res@cnMinLevelValF = 0.
;  res@cnMaxLevelValF = 1.
;  res@cnLevelSpacingF = .1
  plot(1) = gsn_csm_contour_map(wks,GKNT3, res)
;..............................................................................



;      plot(0) = gsn_csm_contour_map(wks,GBIN3, res)
;      plot(1) = gsn_csm_contour_map(wks,GKNT3, res)
;      plot(2) = gsn_csm_contour_map(wks, X01, res)
;*****************************************************************
; Panel   
;*****************************************************************
      dlat    = lat(1)-lat(0)

      resP                     = True
;      resP@gsnPanelMainString  = "MCD19A2: "+yyyy+"/"+month+"/"+day+ \
;                                 ": "+nlat+"x"+mlon+" dlat="+sprintf("%5.3f", dlat)
      resP@gsnPanelMainString  = "MCD19A2:02-06v14-19: "+nlat+"x"+mlon+" dlat="+sprintf("%5.3f", dlat)
      resP@gsnMaximize         = True          ; max size 
     ;resP@gsnPaperOrientation = "portrait"
     gsn_panel(wks,plot,(/2,1/),    resP)

  end if