;*************************************************
; station_1.ncl
;*************************************************
;
; This example reads in station data represented by
; 1D arrays, and generates a filled contour plot
; via the triangular mesh capability in NCL.
;
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
;
; This is a tricky file for NCL to read, because it contains a header
; and then a mix of alpha characters and numeric values. Unfortunately,
; some of the station names contain numbers, so we can't just do a straight
; numeric read of the file. We have to parse out column one first.
;
; There are multiple ways you can read this file, and below
; demonstrates just one of those ways, using various conversion functions
; like stringtofloat and charactertostring.
;
  fname = "pw.dat"
;
; Read in each row as a single string, and then convert to a character
; array so we can parse it.
;
  data  = asciiread("pw.dat", -1, "string")    ; -1 means read all rows.
  cdata = stringtochar(data)
;
; The first row is just header information, so we can discard this.
; That is, we can start with the second row, which is represented
; by index '1'.
;
; The latitude values fall in columns   6-12 (indices 7:13)
; The longitude values fall in columns 13-21 (indices 14:22)
; The pwv data values fall in columns  22-31 (indices 23:)
;
  lat = stringtofloat(charactertostring(cdata(1:,7:13)))
  lon = stringtofloat(charactertostring(cdata(1:,14:22)))
  pwv = stringtofloat(charactertostring(cdata(1:,23:)))

;
; This second file is not so tricky.  The 2D lat/lon data is sorted with 
; lat values first, and then lon values.
;
  nlat     = 70
  nlon     = 70
  latlon2d = asciiread("stn_latlon.dat",(/2,nlat,nlon/),"float")
  lat2d    = latlon2d(0,:,:)
  lon2d    = latlon2d(1,:,:)

  delete(data)        ; Remove arrays we don't need anymore.
  delete(cdata)   
  delete(latlon2d)

  tlat1 = 60.0
  tlat2 = 30.0
  clon  = -98.5
  clat  =  36.3

  wks = gsn_open_wks("ps","station")
  gsn_define_colormap(wks,"WhViBlGrYeOrRe")
  gsn_reverse_colormap(wks)

  res                             = True
  res@gsnFrame                    = False   ; So we can draw markers
  res@gsnMaximize                 = True

  res@gsnSpreadColors             = True
  res@gsnSpreadColorEnd           = 90
  res@pmTickMarkDisplayMode       = "Always"

  res@trGridType                  = "TriangularMesh"    ; The default if you
                                                        ; have 1D data

  res@cnLineLabelPlacementMode    = "Constant"
  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@cnLinesOn                   = True
  res@cnLineLabelsOn              = True
  res@cnLevelFlags		  = new(139,"string")
  res@cnLevelFlags(:)		  = "NoLine"
  res@cnLevelFlags(0::20)         = "LineAndLabel"

  res@lbLabelStride               = 20  
  res@lbBoxLinesOn                = False
  res@tiMainString                = "GPS PWV (18Z)"

  res@sfXArray                    = lon
  res@sfYArray                    = lat

  res@mpProjection                = "LambertConformal"
  res@mpLambertParallel1F         = tlat1
  res@mpLambertParallel2F         = tlat2
  res@mpLambertMeridianF          = clon

  res@mpLimitMode                 = "Corners"
  res@mpLeftCornerLatF            = lat2d(20,18)
  res@mpLeftCornerLonF            = lon2d(20,18)
  res@mpRightCornerLatF           = lat2d(58,52)
  res@mpRightCornerLonF           = lon2d(58,52)

  res@mpFillOn                    = False
  res@mpOutlineDrawOrder          = "PostDraw"
  res@mpFillDrawOrder             = "PreDraw"
  res@mpOutlineBoundarySets       = "GeophysicalAndUSStates"
  res@mpUSStateLineColor          = "Gray10"
  res@mpUSStateLineDashPattern    = 2

;
; Draw the map (frame won't get advanced because gsnFrame was set to False).
;
  map = gsn_csm_contour_map(wks,pwv,res)

;
; Draw markers on the plot in the lat/lon locations.
;
  mkres               = True
  mkres@gsMarkerIndex = 17     ; Filled circle
  mkres@gsMarkerSizeF = 0.03

  gsn_polymarker(wks,map,lon,lat,mkres)

  frame(wks)   ; Now advance the frame.

end