;----------------------------------------------------------------------
; station_1.ncl
;----------------------------------------------------------------------
; Concepts illustrated:
;   - Contouring one-dimensional X, Y, Z data
;   - Reading an ASCII file with several columns of data
;   - Drawing lat/lon locations as filled dots using gsn_coordinates
;   - Controlling which contour lines get drawn
;   - Using opacity to emphasize or subdue overlain features
;   - Reversing a color map
;----------------------------------------------------------------------
; This example reads in station data represented by
; 1D arrays, and generates a filled contour plot over a map.
;
; It uses a newer method (NCL V6.4.0 and later) for specifying the
; lat/lon information via special lat1d/lon1d attributes attached
; to the data to be plotted. See station_old_1.ncl if you have an 
; older version of NCL.
;----------------------------------------------------------------------
; 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"

begin
;
; Data is stored in four columns:  station_name  lat  lon  pwv
; Read in each line as a string, and use "str_get_field" to 
; read in the fields of interest.
;
  fname = "2010_pcp.pnt"
  lines = asciiread(fname,-1,"string")
;
; Use "str_get_field" to indicate which fields to read in.
; Each field is separated by an arbitrary number of spaces.
;
  pwv_dec       = tofloat(str_get_field(lines(1:),16," "))
  pwv_dec@lat1d = tofloat(str_get_field(lines(1:),3," "))  ; Attach as lat1d, lon1d
  pwv_dec@lon1d = tofloat(str_get_field(lines(1:),2," "))  ; for plotting later

;  print(pwv_dec)

  wks = gsn_open_wks("x11","station_precip")       ; send graphics to PNG file

  cmap = read_colormap_file("WhViBlGrYeOrRe") ; read color map
  cmap = cmap(::-1,:)                         ; reverse the color map

  res                             = True
  res@gsnMaximize                 = True
  res@gsnFrame                    = False   ; Want to draw markers later.
  res@gsnDraw                     = False

  res@cnLineLabelPlacementMode    = "Constant"
  res@cnLineLabelFontColor        = "Gray15"
  res@cnLineDashSegLenF           = 0.3

  res@cnLevelSelectionMode        = "ManualLevels"
  res@cnMinLevelValF              = 15         	; 15.25
  res@cnMaxLevelValF              = 50		; 49.75
  res@cnLevelSpacingF             = 0.25

  res@cnFillOn                    = True
  res@cnFillPalette               = cmap(2:88,:)
  res@cnFillOpacityF              = 0.75 
  res@cnLinesOn                   = True
  res@cnLineLabelsOn              = True
  res@cnLevelFlags		  = new(139,"string")
  res@cnLevelFlags(:)		  = "NoLine"
  res@cnLevelFlags(0::20)         = "LineAndLabel"

  res@lbOrientation               = "vertical"
  res@lbBoxLinesOn                = False

;---Zoom in on map area of interest
  res@mpMinLatF                   = min(pwv_dec@lat1d)
  res@mpMinLonF                   = min(pwv_dec@lon1d)
  res@mpMaxLatF                   = max(pwv_dec@lat1d)
  res@mpMaxLonF                   = max(pwv_dec@lon1d)

  res@mpFillOn                    = False
  res@mpOutlineDrawOrder          = "PostDraw"
  res@mpFillDrawOrder             = "PreDraw"

  res@pmTickMarkDisplayMode       = "Always"   ; nicer map tickmarks

;---Create the plot; it won't get drawn because gsnFrame was set to False.
;  res@tiMainString = "GPS PWV (18Z)"
  plot = gsn_csm_contour_map(wks,pwv_dec,res)

;---Draw markers on the plot in the lat/lon locations.
  mkres               = True
  mkres@gsMarkerIndex = 16     ; Filled circle
  gsn_coordinates(wks,plot,pwv_dec,mkres)
end