Here is the text of sonde.ncl. I'll provide a touch additional commentary here.
The basic idea is just as usual. First load the libaries we'll use. Second load the data. Third, parse the data, which in this case is combined in loading, since we have to get rid of the header lines. Then extract the data by defining new variables that mean something to us, including any necessary transformations. Then plot the data.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
filename = "/Users/brian/Documents/ncl_scripts/tutorial_scripts/VBG.txt"
; we know that there are 4 lines of headers,
; which we can read and print out to the screen:
header = readAsciiHead(filename,4)
;print(header)
; now read the rest of the file into an array
data = readAsciiTable(filename,7,"integer",4)
; print the info about the data
;printVarSummary(data)
; pressure is in column 1 (the SECOND column):
pressure = data(:,1)/10.
; temperature is in the fourth column, with a scaling factor of 10:
temperature = data(:,3)/10.
; let's make a simple xy plot of the data
; 1. open a "workstation" -- like figure in Matlab
wks = gsn_open_wks("x11","sonde")
; 2. plot the data
plot = gsn_csm_xy(wks,temperature,pressure,False)
; note the high-level script is from a loaded library
; note that we have no additional "resources" for the plot (False)
; that plot looks pretty bad,
; try marking the correct missing values:
temperature = data(:,3)
temperature@_FillValue = 99999 ;; "@" is for "attribute"
temperature = (/temperature/10./) ;; re-scale, will ignore missing value
; plot again, in same workstation
plot = gsn_csm_xy(wks,temperature,pressure,False)
; now let's flip the Y-axis, so pressure is like height
; we can do this easily by setting a "resource"
res = True ; resource variables are logical
res@trYReverse = True ; then set attributes of resources
plot = gsn_csm_xy(wks,temperature,pressure,res) ;; make sure to change "False" to "res"
; now do it again, but "cheat" by connecting the missing segments
; easiest way to do it is to make new dummy variables that
; don't have the missing values
; 1. find indices where values exist
inds = ind(.not.ismissing(temperature))
; 2. allocate arrays
t = new(dimsizes(inds),typeof(temperature)) ; allocate a new array, size of inds because that's how many non-missing values we have, and data type of temp.
p = t
; 3. fill in arrays
t = temperature(inds)
p = pressure(inds)
; just for fun, at the same time, give
; name and unit attributs to the variables
t@units = "degC"
t@long_name = "Temperature"
p@units = "mb"
p@long_name = "Pressure"
plot = gsn_csm_xy(wks,t,p,res)
; show your cheating ways:
res@xyMarkLineMode = "MarkLines" ; draws markers at data points
res@xyMarkerColor = "red"
res@xyMarker = 16 ; index number tells to draw big dot
plot = gsn_csm_xy(wks,t,p,res)
; now let's do some clean up
; "@tm..." are tickmark resources
res@tmXTBorderOn = False
res@tmXTOn = False
res@tmYRBorderOn = False
res@tmYROn = False
; "@tr..." are transform or axis resources
res@trYMaxF = max(p)
plot = gsn_csm_xy(wks,t,p,res)
; you could try to make things a little different:
td = new(dimsizes(p),"float")
z = td
wd = td
ws = td
td = data(inds,4)/10.
td@long_name = "Dew Point Temperature"
z = data(inds,2)
z@long_name = "Height"
wd = data(inds,5)
ws = data(inds,6)
ws@_FillValue = 99999
ws = (/ws*.1/)
; To delete something, use the delete command. We delete wks here just to make it easy to move on to the next example.
delete(wks)
; Now we'll use the magic of NCL (i.e., the fact that meteorologists are the main users) to make a really fancy looking skew-T plot. This is actually some pretty hard-core NCL, but the code to do it is provided, making our lives much easier. This is taken almost verbatim from: http://www.ncl.ucar.edu/Applications/skewt.shtml
wks = gsn_open_wks("x11","skew")
; --- Create background skew-T and plot sounding----------------
skewtOpts = True
skewtOpts@DrawColAreaFill = True
skewtOpts@DrawFahrenheit = False
skewtOpts@DrawHeightScale = False
skewtOpts@DrawHeightScaleFt = False
skewtOpts@DrawMoistAdiabat = False
dataOpts = False
;dataOpts@PlotWindH = True
;dataOpts@HspdHdir = True
;dataOpts@Height = z
;dataOpts@Hspd = ws
;dataOpts@Hdir = wd
skewt_bkgd = skewT_BackGround(wks,skewtOpts)
skewt_data = skewT_PlotData(wks, skewt_bkgd, p,t,td,z, ws,wd, dataOpts)
draw (skewt_bkgd)
draw (skewt_data)
frame(wks)