;----------------------------------------------------------------------
; This procedure queries a given contour plot for its levels
; and colors, and then it addes filled markers that are 
; colored in based on the levels. The colors should match
; with the colors used in the filled contour plot.
;----------------------------------------------------------------------
procedure add_markers(wks:graphic,plot:graphic,data[*]:numeric,\
                      x[*]:numeric,y[*]:numeric)
local mkres, hollow_markers, filled_markers,levels, nlevels, colors, ncolors, i, ii
begin
;---Get the levels and colors used in the existing contour plot
    getvalues plot
      "cnLevels"     : levels
      "cnFillColors" : colors
    end getvalues

    nlevels = dimsizes(levels)
    ncolors = dimsizes(colors)

;---Create arrays to hold two sets of markers: filled and hollow
    filled_markers = new(ncolors,graphic)
    hollow_markers = new(ncolors,graphic)

;---Set some marker resources
    mkres  = True
    mkres@gsMarkerIndex      = 16           ; filled dot
    mkres@gsMarkerSizeF      = 0.009
    mkres@gsMarkerThicknessF = 1.0

;
; Loop through each level and find all the indexes into the data 
; where the points fall in that level.
;
    do i=0,nlevels
      if(i.eq.0) then
        ii := ind(data.lt.levels(0))                          ; first level
      else if(i.eq.nlevels) then 
        ii := ind(data.ge.levels(nlevels-1))                  ; middle levels
      else
        ii := ind(levels(i-1).le.data.and.data.lt.levels(i))   ; last level
      end if
      end if
      if(.not.any(ismissing(ii))) then
;---Filled marker
        mkres@gsMarkerIndex      = 16
        mkres@gsMarkerColor      := colors(i)
        mkres@gsMarkerSizeF      = 0.009
        mkres@gsMarkerThicknessF = 1.0
        filled_markers(i) = gsn_add_polymarker(wks,plot,x(ii),y(ii),mkres)

;---Hollow marker (makes the markers "pop" a little)
        mkres@gsMarkerIndex      = 4
        mkres@gsMarkerColor      := "black"
        mkres@gsMarkerSizeF      = 0.01
        mkres@gsMarkerThicknessF = 2.0
        hollow_markers(i) = gsn_add_polymarker(wks,plot,x(ii),y(ii),mkres)
      end if
    end do
    plot@$unique_string("hmarkers")$ = hollow_markers    ; make sure these marker ids live outside this procedure
    plot@$unique_string("fmarkers")$ = filled_markers
    return
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
      fname                    = "data.txt"
      lines                    = asciiread(fname,-1,"string")
      lines@_FillValue         = -999
      lines                    = str_sub_str(lines,",-999,",",lines@_FillValue,")
      Depthm                   = tofloat(str_get_field(lines(1:),1," "))
      T_degC                   = tofloat(str_get_field(lines(1:),2," "))
      Station                  = tofloat(str_get_field(lines(1:),6," "))

      wks = gsn_open_wks("png","plot_with_markers")

      res = True
      res@gsnMaximize          = True
      res@gsnDraw              = False
      res@gsnFrame             = False
     
      res@sfXArray             = Station
      res@sfYArray             = Depthm

      res@cnLevelSelectionMode = "ManualLevels"
      res@cnMinLevelValF       = 6.
      res@cnMaxLevelValF       = 14.
      res@cnLevelSpacingF      = 1.0 ; contour spacing
      res@cnFillOn             = True
      ; use a cmocean colormap
      res@cnFillPalette        = "cmocean_thermal"
      res@cnLinesOn            = False
      res@lbOrientation        = "Vertical"

      res@tiMainString         = "CalCOFI line 90 section"
      res@tiYAxisString        = "Depth"       ; y-axis title
      res@tiXAxisString        = "Station"       ; x-axis title

;----------------------------------------------------------------------
; Create filled contour plot; no axes are reversed yet.
;----------------------------------------------------------------------
      plot_norev  = gsn_csm_contour(wks,T_degC,res)

      draw(plot_norev)
      frame(wks)

;----------------------------------------------------------------------
; Create plot again, and then add colored markers so we can compare 
; against original contour plot. This is mainly for debug purposes.
;----------------------------------------------------------------------
      res@cnFillOpacityF        = 0.5        ; make contours partially transparent so we can see markers
      res@lbOverrideFillOpacity = True       ; don't make labelbar transparent like contours

      plot_norev_with_markers  = gsn_csm_contour(wks,T_degC,res)
      add_markers(wks,plot_norev_with_markers,T_degC,Station,Depthm)

      draw(plot_norev_with_markers)
      frame(wks)

;----------------------------------------------------------------------
; Create plot again, but with two axes reversed. This causes the
; contour plot to be incorrect. Make the bad contours transparent and 
; add the locations of the data as markers colored based on the colors
; and levels used in the contour plot.
;----------------------------------------------------------------------
      res@trYReverse            = True       ; reverses y-axis
      res@trXReverse            = True       ; reverses x-axis 
      res@cnFillOpacityF        = 0.0        ; make contours transparent

      plot_rev = gsn_csm_contour(wks,T_degC,res)
      add_markers(wks,plot_rev,T_degC,Station,Depthm)

      draw(plot_rev)
      frame(wks)

;---Panel all three plots for easier comparison
      pres                 = True
      res@gsnMaximize      = False
      pres@gsnPanelRowSpec = True
      gsn_panel(wks,(/plot_norev,plot_norev_with_markers,plot_rev/),(/2,1/),pres)
end