;***************************************************************
; binning_2.ncl
;
; Concepts illustrated:
;   - Create an array that spans the desired area
;   - Read data [ here, create bogus data]
;   - Loop over data and count instances of occurence
;   - Plot the data
;
;**************************************************************************
;--- Create desired grid. 
;**************************************************************************
  latS =  0
  latN =  50
  lonW =  90
  lonE =  180

  dlat =  5.0                    
  dlon =  5.0

  nlat = toint((latN-latS)/dlat) + 1
  mlon = toint((lonE-lonW)/dlon) + 1

  lat  = fspan(latS, latN, nlat)
  lon  = fspan(lonW, lonE, mlon)

  lat@units = "degrees_north"
  lon@units = "degrees_east"

  count     = new( (/nlat,mlon/), "float", 1e20)   ; [lat | 11] x [lon | 25]
  count!0   = "lat"
  count!1   = "lon"
  count&lat =  lat
  count&lon =  lon
  printVarSummary(count)

  valavg    = count

;********************************************************************
;--- Read data from csv file
;********************************************************************
  
  filename = "par_jtwc_above_ts_1979-1993.csv"
  lines    = asciiread(filename,-1,"string")

  data = lines(1:)
  clat = tofloat(str_get_field(data,7,","))    
  clon = tofloat(str_get_field(data,8,","))
  cval = tofloat(str_get_field(data,9,","))    ; vmax

  printVarSummary(clat)
  printMinMax(clat,0)
  printVarSummary(clon)
  printMinMax(clon,0)
  printVarSummary(cval)
  printMinMax(cval,0)                          ; min=35   max=145
  print("nval="+nval)

  vMax = 35
  nval = num(cval.ge.vMax)                     ; num(.not.ismissing(cval))
  print("nval="+nval)

;********************************************************************
;--- Bin count and sum; This assumes a simple rectilinear grid
;********************************************************************
  count  = 0
  valavg = 0
 
  cval@_FillValue = -999.0
  cval = where(cval.ge.vMax, cval, cval@_FillValue)
  printMinMax(cval,0)

  npts = dimsizes(cval)

  do n=0,npts-1
     if (clat(n).ge.latS .and. clat(n).le.latN .and.  \
         clon(n).ge.lonW .and. clon(n).le.lonE .and.  \
         .not.ismissing(cval(n)) ) then               ; cval(n).ge.vMax)  then

         jl = toint((clat(n)-latS)/dlat) 
         il = toint((clon(n)-lonW)/dlon) 
         count(jl,il)  = count(jl,il)  + 1
         valavg(jl,il) = valavg(jl,il) + cval(n)

     end if
  end do

  count@long_name = "Occurrence Count: vMax="+vMax
  count@units     = ""

  printVarSummary(count)
  print("count: min="+min(count)+"   max="+max(count))

  count  = where(count.eq.0, count@_FillValue,count)   ; don't divide by 0         
;********************************************************************
;--- Average
;********************************************************************
  
  valavg = valavg/count
  valavg@long_name = "Average"
  valavg@units     = "kts"

  printVarSummary(valavg)
  print("valavg: min="+min(valavg)+"   max="+max(valavg))

;********************************************************************
;--- Bin frequency (%)
;********************************************************************

  freq   = count
  nFilt  = num(.not.ismissing(cval))

  freq   = (count/nFilt)*100
  freq@long_name = "frequency"
  freq@units = "%"

  printVarSummary(freq)
  print("freq: min="+min(freq)+"   max="+max(freq))

  test_lat =  20.0
  test_lon = 130.0  
  print("---")
  print("---> ("+test_lat+","+test_lon+")")
  print("count ="+count({test_lat},{test_lon}))
  print("valavg="+valavg({test_lat},{test_lon}))
  print("  freq="+  freq({test_lat},{test_lon}))
  print("---")

;************************************************
; create plot
;************************************************
  freq  = where(freq .eq.0,  freq@_FillValue, freq)   ; 

  plot  = new (3, "graphic")

  wks = gsn_open_wks("png","binning")            ; send graphics to PNG file

  res                       = True     ; plot mods desired
  res@gsnAddCyclic          = False    
  res@gsnDraw               = False    
  res@gsnFrame              = False    

  res@cnFillOn              = True     ; turn on color fill
  res@cnFillPalette         = "BlAqGrYeOrRe"     ; set color map
  res@cnFillMode            = "RasterFill"       ; Raster Mode
  res@cnLinesOn             = False    ; turn of contour lines

  res@lbOrientation         = "vertical"

  res@mpMinLatF             = latS
  res@mpMaxLatF             = latN
  res@mpMinLonF             = lonW
  res@mpMaxLonF             = lonE
  res@mpCenterLonF          = (lonE+lonW)*0.5
  res@mpGridAndLimbOn       = True  
  res@mpGridLineDashPattern = 2             ; Dashed lines
  res@mpGridLatSpacingF     = 5.0
  res@mpGridLonSpacingF     = 10.0

  res@cnLevelSpacingF       = 5.0      ; contour spacing
  res@gsnCenterString       = "Occurrence Count"
  plot(0) = gsn_csm_contour_map_ce(wks,count, res)
  
  res@cnLevelSpacingF       = 0.5     ; contour spacing
  res@gsnCenterString       = "Frequency (%)"
  plot(1) = gsn_csm_contour_map_ce(wks,freq , res)
  
  res@cnLevelSpacingF       =  5.0     ; contour spacing
  res@gsnCenterString       = "Average"
  plot(2) = gsn_csm_contour_map_ce(wks,valavg,res)

  resP                      = True               ; modify the panel plot
  resP@gsnMaximize          = True
;;resP@gsnPanelMainString   = "A common title"
  gsn_panel(wks,plot,(/3,1/),resP)               ; now draw as one plot