;*************************************************
; rdm2grid_1.ncl
;*************************************************
;
; Concepts illustrated:
;   - Reading SLP from netCDF
;   - Using array syntax[ ::-1 ] so that lat is South-to-North
;   - Generating random index values for demo
;   - Using natgrid to interpolate to grid
;   - Plotting original and reconstructed grids

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

undef("natgrid_Wrap")
function natgrid_Wrap (xi[*]:numeric,yi[*]:numeric, fi:numeric \
                      ,xo[*]:numeric,yo[*]:numeric)

local fo, dimfi, nDimi, dimfo, nDimo
begin
  fo = natgrid (xi,yi,fi, xo,yo)            ; perform interpolation 
                                             
  dimfi = dimsizes(fi)
  nDimi = dimsizes(dimsizes(fi))            ; rank: # of dimensions
  dimfo = dimsizes(fo)
  nDimo = dimsizes(dimsizes(fo))            ; rank: # of dimensions

  copy_VarAtts (fi, fo)                     ; copy variable attributes
  copy_VarCoords_2 (fi, fo)                 ; copy left coords if present

  if (isdimnamed(xo,0)) then
      fo!(nDimo-2) = xo!0
      fo&$xo!0$    = xo  
  else
      fo!(nDimo-2) = "x"  
      fo&x         = xo  
  end if

  if (isdimnamed(yo,0)) then
      fo!(nDimo-1) = yo!0
      fo&$yo!0$    = yo  
  else
      fo!(nDimo-1) = "y"  
      fo&y         = yo  
  end if

  return(fo)        ; NOTE: the last two are ([...,]x,y)
                    ; may have to reverse upon return
end

begin

  NOBS  = 500      

;------------------------------------------------------------------
; Read data
;------------------------------------------------------------------

  diri = "/Users/shea/Data/CDC/"  
  fili = "slp_1988mm.nc"
  f    = addfile(diri+fili, "r")
  Z    = f->slp(0,::-1,:)      ; natgrid requires monotonic increasing y [lat]

;------------------------------------------------------------------
; Generate N random values from the above
;------------------------------------------------------------------

  dimz  = dimsizes( Z )
  nlat  = dimz(0)
  mlon  = dimz(1)

  N     = nlat*mlon
  iZ    = generate_unique_indices( N )
 ;iZ    = floattoint(random_uniform (0, N-1, N))  
  trip  = grid2triple (Z&lon, Z&lat, Z) 

  dlon  = Z&lon(2)-Z&lon(1)
  reps  = random_uniform (-dlon, dlon, N)  ; minor location perturbations

  IOBS  = iZ(0:NOBS)   ; convenience
                       ; clarity only ... create explicit variables  
  rlon  = trip(0,IOBS) + reps(IOBS)  
  rlat  = trip(1,IOBS) + reps(IOBS(::-1))
  rZ    = trip(2,IOBS)

  rlon@units = "degrees_east"
  rlat@units = "degrees_north"

;------------------------------------------------------------------
; natgrid
;------------------------------------------------------------------
  zgrd  = natgrid_Wrap (rlon, rlat, rZ, Z&lon,Z&lat)
  znat  = zgrd(lat|:,lon|:)

;------------------------------------------------------------------
; plot
;------------------------------------------------------------------
  maxLev   = 16
  mnmxint  = nice_mnmxintvl( min(znat), max(znat), maxLev, True)

  dimznat  = dimsizes(znat)
  kplt     = dimznat(0)

  wks      = gsn_open_wks("ps","rdm2grid")        
  gsn_define_colormap(wks,"amwg")
  plot     = new(2,graphic)           
  
  res                     = True
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnFillOn            = True            ; turn on color
 ;res@cnLinesOn           = False           ; turn off contour lines
  res@cnLineLabelsOn      = False           ; turn off contour line labels

  res@gsnSpreadColors     = True            ; spread out color table
  res@lbLabelBarOn        = False           ; turn off individual cb's
  res@mpCenterLonF        = 180.

  res@cnLevelSelectionMode =  "ManualLevels"   
  res@cnMinLevelValF       = mnmxint(0) 
  res@cnMaxLevelValF       = mnmxint(1) 
  res@cnLevelSpacingF      = mnmxint(2) 

  pmres = True
  pmres@gsMarkerIndex      = 17             ; Filled circle
  pmres@gsMarkerSizeF      = 0.0125

  res@gsnCenterString = "Original grid"
  plot(0) = gsn_csm_contour_map_ce(wks,Z,res) 
   
  res@gsnCenterString = "NATGRID: N="+NOBS
  plot(1) = gsn_csm_contour_map_ce(wks,znat,res) 
          ; add polymarkers
  dum1    = gsn_add_polymarker(wks, plot(1), rlon, rlat, pmres)

;************************************************
; create panel
;************************************************
  resP                     = True               ; modify the panel plot
 ;resP@txString            = title
  resP@gsnMaximize         = True               ; make large
  resP@gsnPanelLabelBar    = True               ; add common colorbar
  resP@lbLabelAutoStride   = True               ; Let NCL decide spacing
 ;resP@lbLabelStride       = 2                  ; force every other label 
 ;resP@lbLabelFontHeightF  = 0.0125             ; make labels smaller [0.2 default]

  gsn_panel(wks,plot,(/2,1/),resP)              ; now draw as one plot
end