;----------------------------------------------------------------------
; Procedure for adding a labelbar at a given NDC location, given
; the levels and colors to use.
;----------------------------------------------------------------------
undef("add_labelbar")
procedure add_labelbar(wks,plot,levels,colors)
local lbres, labels
begin
  nlevels = dimsizes(levels)
;----------------------------------------------------------------------
; Draw a labelbar
;----------------------------------------------------------------------
  lbres                    = True
  lbres@vpWidthF           = 0.80             ; width
  lbres@vpHeightF          = 0.10             ; height
  lbres@lbPerimOn          = False            ; Turn off perimeter.
  lbres@lbOrientation      = "Horizontal"     ; Default is vertical.
  lbres@lbLabelAlignment   = "InteriorEdges"  ; Default is "BoxCenters"
  lbres@lbFillColors       = colors           ; Colors for boxes.
  lbres@lbMonoFillPattern  = True             ; Fill them all solid.
  lbres@lbLabelFontHeightF = 0.012            ; label font height

  labels = sprintf("%4.2f",levels)

  lbid = gsn_create_labelbar(wks,nlevels+1,labels,lbres)
   
;
; Now, create some annotation resources indicating how we want to
; attach the labelbar to the plot. Here, we are using the top right
; corner of the labelbar as the point which we are going to position
; it, and then we use amParallelPosF and amOrthogonalPosF to indicate
; where we want to place it.
;
; amParallelPosF/amOrthogonalPosF                                       ;
;    0.0/ 0.0  -  annotation in dead center of plot
;    0.5/ 0.5  -  annotation at bottom right of plot
;    0.5/-0.5  -  annotation at top right of plot
;   -0.5/-0.5  -  annotation at top left of plot
;   -0.5/ 0.5  -  annotation at bottom left of plot
;                                                                
  amres                  = True
  amres@amJust           = "TopCenter"
  amres@amParallelPosF   =  0.0   ; keep labelbar centered
  amres@amOrthogonalPosF =  0.6   ; move down and outside of plot

;
; Give both annotation id and labelbar id unique names.
;
; Attaching them to plot with unique names ensures that
; labelbar "lives" outside this procedure.
;
  tmpid1 = "anno"+unique_string("id")
  tmpid2 = "lbar"+unique_string("id")

  plot@$tmpid1$ = gsn_add_annotation(plot,lbid,amres)
  plot@$tmpid2$ = lbid
end
   

    levels =  ispan(5,50,5)
    nlevels = dimsizes(levels)
    print(nlevels)
    colors  = span_color_rgba("NCV_jet",nlevels+1)

    R = (/16,0,39,7,47,5,7,30,0,17,0,7,3,0,0,47,0,0,46,32,0,46,3,7,0,0,7,0,25,0,46,4,0,7,7,0,0,0,47,46,47,7,0,3/)
    lat = (/27.156,27.883,25.44,30.383,31.71,26.681,27.567,28.422,24.133,25.233,23.233,30.733,26.333,26.167,24.917,28.585,25.933,26.3,24.744,26.74,27.383,39.179,26.533,25.491,26.404,29.7,22.533,23.4,25.45,26.867,26.761,30.867,27.233,25.033,29.017,26.667,26.117,29.617,22.655,30.333,25.591,25.267,26.983,25.3/)
    lon = (/77.961,78.067,81.734,76.767,74.797,88.329,81.6,79.451,88.267,86.95,87.85,76.883,89.497,85.9,84.183,77.206,80.833,87.267,84.951,83.45,80.167,75.755,88.717,78.558,80.41,77.033,88.333,88.517,82.867,80.933,80.889,75.933,79.05,88.133,77.717,84.917,85.4,78.383,88.447,76.467,85.088,87.467,84.85,83.017/)
    npts = 44
 
    num_distinct_markers = nlevels+1        ; number of distinct markers
    lat_new = new((/num_distinct_markers,npts/),float,-999)
    lon_new = new((/num_distinct_markers,npts/),float,-999) 

    do i = 0, num_distinct_markers-1
    if (i.eq.0) then
      indexes = ind(R.lt.min(levels))
    end if
    if (i.eq.num_distinct_markers-1) then
      indexes = ind(R.ge.max(levels))
    end if
    if (i.gt.0.and.i.lt.num_distinct_markers-1) then       
      indexes = ind(R.ge.levels(i-1).and.R.lt.levels(i))
    end if
   
    if (.not.any(ismissing(indexes))) then 
      npts_range = dimsizes(indexes)   ; # of points in this range.
      lat_new(i,0:npts_range-1) = lat(indexes)
      lon_new(i,0:npts_range-1) = lon(indexes)
    end if
    delete(indexes)            ; Necessary b/c "indexes" may be a different
                               ; size next time.
  end do 

  wks = gsn_open_wks("png","polyg")              ; send graphics to PNG file
  
  ;---Set up some map resources.
  mpres                       = True
  mpres@gsnMaximize           = True             ; Maximize plot in frame.
  mpres@gsnDraw               = False            ; Will draw later
  mpres@gsnFrame              = False            ; Don't advance the frame
  mpres@pmTickMarkDisplayMode = "Always"

  ;---Zoom in on United States.
  mpres@mpMinLatF             = 20.
  mpres@mpMaxLatF             = 37.
  mpres@mpMinLonF             = 68.
  mpres@mpMaxLonF             = 90.

 ; mpres@tiMainString          = "Dummy station data colored and~C~sized according to range of values"
  map = gsn_csm_map(wks,mpres)

  ;--Create logical variables to hold the marker resources.
  gsres               = True
  gsres@gsMarkerIndex = 16          ; Use filled dots for markers.


  ; Loop through each grouping of markers, and draw them one set at
; a time, assigning the proper color and size with gsn_marker.
;
  base_size = 0.01
  pmid = new(num_distinct_markers,graphic)
  do i = 0, num_distinct_markers-1
    if (.not.ismissing(lat_new(i,0)))
      gsres@gsMarkerColor      = colors(i,:)
      gsres@gsMarkerSizeF      = base_size * (i+1)/3.
      gsres@gsMarkerThicknessF = 0.7*(i+1)
      pmid(i) = gsn_add_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres)
    end if
  end do

filenames = "IND_WHOLE" + ".shp"
;;---Attach two sets of polylines to the same map.
;  lnres             = True
;  lnres@gsLineColor = "black"
;  lnres@gsLineThicknessF=1.8
;  poly1 = gsn_add_shapefile_polylines(wks,map,filenames,lnres)


;---Draw labelbar and advance frame.
  add_labelbar(wks,map,levels,colors)
  draw(map)
  frame(wks)