;***************************************************************
; 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
;
;***************************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"    

;**************************************************************************
;--- Create desired grid. Here, 2x2 but can be (say) 1x3 if dlat=1, dlon= 3 
;**************************************************************************
  latS =  0
  latN =  50
  lonW =  60
  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) 
  count!0   = "lat"
  count!1   = "lon"
  count&lat =  lat
  count&lon =  lon

  valavg    = count

;********************************************************************
;--- Read data from csv file
;********************************************************************
  
  filename = "test.txt"
  lines    = asciiread(filename,-1,"string")   ; Read data as lines of strings
  nrows    = dimsizes(lines)                   ; This includes the header line
  
  delim = " "
  header = str_split_csv(lines(0),delim,0)  ; SN, CY, Year, Month, Day, Hour, Lat, Lon, VMax, Cat 
  
;---Read data by selecting column and converting to int or float, or leaving as a string.
  sn     = toint(str_get_field(lines(1:),1,delim))
  cy     = toint(str_get_field(lines(1:),2,delim))
  year   = toint(str_get_field(lines(1:),3,delim))
  month  = toint(str_get_field(lines(1:),4,delim))
  day    = toint(str_get_field(lines(1:),5,delim))
  hour   = toint(str_get_field(lines(1:),6,delim))
  clat    = tofloat(str_get_field(lines(1:),7,delim))
  clon    = tofloat(str_get_field(lines(1:),8,delim))
  vmax   = tofloat(str_get_field(lines(1:),9,delim))
  cat    =         str_get_field(lines(1:),10,delim)
  
  clat@units = "degrees_north"
  clon@units = "degrees_east"
  printMinMax(vmax,0)
  printMinMax(lat,0)
  printMinMax(lon,0) 

  ;vmax = tofloat(str_get_field(data,9,","))
  
  ;filt = mask(data,(vmax.ge.35),True)
  
  clon = where(clon.lt.lonW, lonW, clon)  
  clon = where(clon.gt.lonE, lonE, clon)
  clat = where(clat.lt.latS, latS, clat)
  clat = where(clat.gt.latN, latN, clat)

;********************************************************************
;--- Bin count and sum; This assumes a simple rectilinear grid
;********************************************************************
  count  = 0
  valavg = 0
  

  sn_unique = get_unique_values(sn)
  snpts       = dimsizes(sn_unique)
  do a=0,snpts-1
         ii := ind(sn.eq.sn_unique(a)) 
  end do

  npts = dimsizes(ii)
  nlat=lat(ii)
  nlon=lon(ii)
  nvmax=vmax(ii)

  do n=0,npts-1  
           if (nlat(n).ge.latS .and. nlat(n).le.latN .and.  \
                 nlon(n).ge.lonW .and. nlon(n).le.lonE .and.  \
                .not.ismissing(nvmax(n)) ) then

                 jl = toint((nlat(n)-latS)/dlat) 
                 il = toint((nlon(n)-lonW)/dlon) 
                 count(jl,il)  = count(jl,il)  + 1
                 valavg(jl,il) = valavg(jl,il) + nvmax(n)
              end if
  end do

 ;count@long_name = "Occurrence Count"
 ;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 = "..."
 ;valavg@units     = "..."

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

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

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

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

;************************************************
; 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       = 5.0     ; 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