Here’s a (slight hack of an) example of displaying two time series (5-min mean wind speed and wind direction), each with their own dependent axis, on a set of two independent axes (local time and universal time) using
.
For data, I’m using a netCDF archive file from the NCAR Mesa Lab (ML) weather station. Information about the weather station, including a link to the FTP site to download other archived data files, can be found here:
.
You can get the data I used in this example, as well as the program on which the example is based,
After the break, I’ll explain in detail how I made this plot.
I chose a
file from the archive on the date of a
chinook, or windstorm, in Boulder. I think it’s neat to see how the weather station at the ML, which is situated close to the Flatirons, captures the winds in the storm
[1].
Start by reading the “time_offset”, “wspd” and “wdir” variables from the archive file:
ml_file = file_which('mlab.20100615.cdf', /include_current_dir)
ml_id = ncdf_open(ml_file)
ncdf_varget, ml_id, ncdf_varid(ml_id, 'time_offset'), time
ncdf_varget, ml_id, ncdf_varid(ml_id, 'wspd'), wspd
ncdf_varget, ml_id, ncdf_varid(ml_id, 'wdir'), wdir
ncdf_close, ml_id
TIME is an offset, in seconds, from 00 UTC (18 local time). Convert it to hours:
time /= 60.0^2
Plot the wind speeds first, in red, overriding the default plot margins (expressed in normalized coordinates) to leave room for the extra axes.
plot_margin = [0.15, 0.25, 0.15, 0.15]
plot_xrange = [0,24]
p_wspd = plot(time, wspd, 'r', $
axis_style=1, $ ; make only x & y axes, not box axes
margin=plot_margin, $
xrange=plot_xrange, $
xmajor=9, $ ; 3-hr intervals
dimensions=[700,600], $ ; embiggen window to fit extra axes
name='Speed', $
xtitle='Time (UTC)', $
ytitle='Wind Speed ($m s^{-1}$)', $
title='NCAR Mesa Lab Weather Station Winds!C2010-06-15')
Here, the short code “r” makes a red plot line. Setting AXIS_STYLE to 1 shows only left and bottom plot axes. NG allows TeX-like format codes in text annotations, like the superscript in the y-axis title. Looking ahead, the XRANGE and MARGIN values will be needed in the plot of wind direction and the NAME property will be used in the plot legend.
Next, calculate local time from the UTC values in the wind speed plot P_WSPD. Use AXIS to display these values with a second x-axis, positioned just below the first:
local_time = strtrim(round(p_wspd.xtickvalues + 18) mod 24, 2)
a_time = axis('x', $
tickname=local_time, $
location=[0,min(p_wspd.yrange)-2,0], $ ; data coordinates
title='Time (LST)')
I used the minimum y-axis value from P_WSPD plus an offset to position this axis slightly below the first. A better technique would be to calculate the offset as a fraction of the entire y-axis range of P_WSPD. (I’ll show this in a subsequent post.)
Now plot the wind direction, in blue, in the same window as P_WSPD:
p_wdir = plot(time, wdir, 'b', $
/current, $
axis_style=0, $ ; display no axes
margin=plot_margin, $ ; need to use the same margin as above
xrange=plot_xrange, $ ; and the same xrange
yrange=[0,360], $
name='Direction')
The CURRENT keyword (in lieu of OVERPLOT) is the key here: it puts this plot in the same window as the first, but it doesn’t use the same data coordinate system; that’s why I have to use the same MARGIN and XRANGE as in P_WSPD to align the plot correctly in the window.
Because I chose not to display axes in P_WDIR, use AXIS to make a y-axis on the right side of the plot:
a_wdir = axis('y', $
target=p_wdir, $
major=5, $ ; [0, 90, 180, 270, 360]
minor=2, $
location=[max(p_wdir.xrange),0,0], $ ; right axis, data coordinates
textpos=1, $ ; text faces outward
tickdir=1, $ ; ticks face inward
title='Wind Direction (deg from N)')
The TARGET keyword picks up data coordinates from P_WDIR.
Note how the wind speeds peak when the wind direction is from 270 degrees, or west – a sure sign of a windstorm in Boulder!
To finish, display a legend for the two plots in the default location:
!null = legend(target=[p_wspd, p_wdir])
Again, TARGET picks up properties of the plots P_WSPD and P_WDIR for display. I chose to throw away the reference returned from LEGEND by directing it to !NULL.
Save this visualization to an encapsulated PostScript file with:
p.save, 'ng_multiple_axes.eps'
So, I’ve shown a slightly hacky (although not necessarily worse than some of the DG plots I’ve made) way of producing a NG plot with multiple axes. Please use this as a guide for making your own NG plots.