; Ozone site plotting for R values of DEC verses GEOS-Chem scatter plot values to display the spatial difference across NYS.
;
;
;
begin
;---------------------------
; Loading input data file
    strs = asciiread("ozonetest.txt",-1,"string")

	delim = " "
	nfields = str_fields_count(strs(0), delim)
	
;	print(nfields)

	fieldavg =  11  ;avged data
	R = tofloat(str_get_field(strs,10," "))
	print(R)


   lat       = tofloat(str_get_field(strs,2," "))
	print(lat)
;	Nlat = new((/51,1/),"float",-999.)
	Nlat = lat
;	print(Nlat) 
;      lat@_FillValue=integertoshort(-999)
   lon       = tofloat(str_get_field(strs,3," "))
	print(lon)
;	Nlon = new((/51,1/),"float",-999.)
	Nlon = lon
;	print(Nlon)
;       lon@_FillValue=integertoshort(-999)



;---------------------------
; Set up markers

	arr = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1./) 
;	colors = read_colormap_file("WhViBlGrYeOrRe") ; setting color markers
        colors = (/"10","30","38","48","56","66","74","94","450","600"/)
	labels = new(dimsizes(arr)+1,string)  ; Labels for legend.
;-------------------------------------------------------------

num_distinct_markers = dimsizes(arr)+1      ; number of distinct markers
  lat_new = new((/num_distinct_markers,dimsizes(R)/),float,-999)
  lon_new = new((/num_distinct_markers,dimsizes(R)/),float,-999)
;	print(arr)
;	print(lat_new)
;	print(lon_new)

;
; Group the points according to which range they fall in. At the
; same time, create the label that we will use later in the legend.
;       
  do i = 0, num_distinct_markers-1
    if (i.eq.0) then
      indexes = ind(R.lt.arr(i))
      labels(i) = "x < " + arr(i)
    end if
    if (i.eq.num_distinct_markers-1) then
      indexes = ind(R.ge.max(arr))
      labels(i) = "x >= " + max(arr)
    end if
    if (i.gt.0.and.i.lt.num_distinct_markers-1) then       
      indexes = ind(R.ge.arr(i-1).and.R.lt.arr(i))
      labels(i) = arr(i-1) + " <= x < " + arr(i)
    end if
;
; Now that we have the set of indexes whose values fall within 
; the given range, take the corresponding lat/lon values and store
; them, so later we can color this set of markers with the appropriate
; color.
;
    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


;============================================================
; Start the graphics.
   wks      = gsn_open_wks("png","O3-NYS-DEC-GEOS-Rvals")         ; send graphics to PNG file
;  gsn_define_colormap(wks,"gui_default")      ; define a different color map
  gsn_define_colormap(wks,"WhViBlGrYeOrRe")                
;
; Set up some map resources.
;
  mpres              = True
  mpres@gsnMaximize  = True             ; Maximize plot in frame.
  mpres@gsnFrame     = False            ; Don't advance the frame


  mpres@cnFillOn = True
  mpres@lbLabelBarOn = True
  mpres@lbAutoManage = False
  mpres@pmLabelBarDisplayMode = "Always"
  mpres@lbOrientation       = "vertical"

; We want county outlines.
   mpres@mpDataBaseVersion     = "MediumRes"
   mpres@mpDataSetName         = "Earth..2"
   mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates"
   mpres@mpOutlineSpecifiers   = (/"Land","New York:counties"/)


; Zoom in on NY, U.S.
   mpres@mpLimitMode           = "LatLon"
   mpres@mpMinLatF             = 40
   mpres@mpMaxLatF             = 46
   mpres@mpMinLonF             = -80.0
   mpres@mpMaxLonF             = -71.0

; Turn on tickmarks.
   mpres@gsnTickMarksOn        = True

  mpres@mpFillColors = (/"transparent","transparent", \
  "lightgray","transparent"/)    ;assign light gray to land masses

; Main title
  mpres@tiMainString          = "DEC Ozone Site R values"
  map = gsn_csm_map(wks,mpres)

;
;  Create logical variables to hold the marker and text resources.
;  These markers are different than the XY markers, because they are not
;  associated with an XY plot. You can put these markers on any plot.
;
  gsres               = True
  gsres@gsMarkerIndex = 16          ; Use filled dots for markers.

  txres               = True
  txres@txFontHeightF = 0.015

;
; Loop through each grouping of markers, and draw them one set at
; a time, assigning the proper color and size with gsn_marker.
;
; At the same time, draw a legend showing the meaning of the markers.
;

;  xleg = (/0.1,0.1,0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.8,0.9,0.9/)   ; Location of
;  xtxt = (/0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.7,0.9,0.9,1.0,1.0/)   ; legend markers
;  xleg = (/0.1,0.1,0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.8,0.9,0.9/)   ; Location of
;  xtxt = (/0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.7,0.9,0.9,1.0,1.0/)   ; legend markers
;  yleg = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/)   ; and text
;  ytxt = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/)   ; strings.

  do i = 0, num_distinct_markers-1
   if (.not.ismissing(lat_new(i,0)))
      gsres@gsMarkerColor      = colors(i)
      gsres@gsMarkerThicknessF = 0.7*(i+1)
;      gsres@gsMarkerThicknessF = 0.7
      gsn_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres)
;
; Add marker and text for the legend.

;
;      gsn_polymarker_ndc(wks,          xleg(i),yleg(i),gsres)
;      gsn_text_ndc      (wks,labels(i),xtxt(i),ytxt(i),txres)
    end if
  end do

  frame(wks)    ; Advance the frame.
end