| This post was kindly contributed by SAS and R - go there to comment and to read the full post. |
Before explaining how to make and interpret the plot above, Nick and I want to make a plea for questions–it’s hard to come up with useful questions to explore each week!
As shown in Example 8.12, data from the Cyclemeter app can be used to make interesting plots of a bike ride. But in the previous application, questions remain. Why did my speed vary so much, for example? And why do some sections of the route appear to be very straight while others are relatively bumpy? To investigate these questions, we’ll add some more data into the plot, using color. First, I’ll shade the line according to my elevation– my speed is probably faster when I’m going down hill. Then I’ll add the points where the GPS connected; this may explain the straight sections. Doing all of this will require more complicated coding tasks.
(Data were read in previously.)
SAS
My initial thought was that I would use proc univariate to generate a histogram (section 5.1.4) and generate the color for the line from the histogram categories. But it turns out that you can’t force the histogram to have your desired number of categories. (This is true for the R hist() function as well.) So I’m going to find the minimum and maximum elevation using proc summary (section 2.1), saving the results in a data set. Then I’ll add those min and max values to each line of the data, using the tip found here. Finally, I’ll make my own categories for elevation by looping through the category boundaries until I find one bigger than my observation.
proc summary data=ride;
var elevation__feet_;
output out=minmax max=maxelev min=minelev;;
run;
data aride;
set ride;
setme = 1;
set minmax point=setme;
colorindex = 0;
range = maxelev-minelev;
do i = 0 to 5;
if elevation__feet_ ge (minelev + (i * range/6) )
then colorindex = i + 1;
end;
run;
To make the plot, I’m going to use the annotate facility, as before. However, the annotate %line macro won’t work, for reasons I won’t go into. So I need to make an annotate data set directly. Briefly, the way annotate data sets work is that they have a pre-specified variable list in which only certain values are allowed. Some of these are discussed in section 5.1, but in addition to these, below we use the function variable; when function="MOVE" the conceptual plotting pen moves without drawing. When function="DRAW" a line is plotted from the last location to the current one. Locations are determined by x and y variables. The line and size variables affect the line style and thickness.
The other interesting bit of the code below is the creation and use of a temporary array (section 1.11.5) for the colors. I choose the color from this array by indexing on the colorindex variable created above.
For help on the annotate facility, see Contents; SAS Products; SAS/GRAPH; The Annotate Facility. Colors are discussed in section 5.3.11, but we use them in a slightly different way, here. For help with color specification, see Contents; SAS Products; SAS/GRAPH; Concepts; Colors.
data annoride2;
set aride;
if elevation__feet_ ne .;
retain xsys '2' ysys '2' hsys '6';
array mycolor [6] $ _temporary_
("CX252525" "CX636363" "CX969696" "CXBDBDBD" "CXD9D9D9" "CXF7F7F7");
x = longitude;
y = latitude;
/* If this is the first observation, we need to move the pen to
the first point. Otherwise we draw a line */
if _n_ ne 1 then function = "DRAW";
else function = "MOVE";
color = mycolor[colorindex];
line = 1;
size = speed__miles_h_ * 2;
output;
run;
Finally, we can plot the figure. I use the symbol statement (section 5.2.2) to plot nice dots where the GPS connected, and the axis statement (section 5.3.7) to remove the values of latitude and longitude, which don’t interest me much. Analysis of the graph I leave for the R version.
axis1 major=none minor=none value=none;
symbol1 i=none v=dot c=black h=.2;
proc gplot data = ride;
plot latitude * longitude /
vaxis=axis1 haxis=axis1 annotate=annoride2;
run;
quit;
The resulting plot is shown above.
R
In R, I’m going to modify my vectorized version of the lines() function to additionally assign colors to each line. As in the SAS example, I will use categories of elevation to do this. Assigning the elevations to categories is a little trickier if I want to avoid for() loops. To do it, I will first find the category boundaries. I’ll then use the which() function (as in section 1.13.2) to figure out the category boundaries smaller than the observation, and the max() function (section 1.8.1) to select the largest of these. Unfortunately, I also need the sapply() function (section B.5.3) so that I can repeat this process for each elevation and return a vector. I’ve annotated below to explain how this works. Finally, I use the RColorBrewer package to generate a vector of colors. (I also used it to find the colors I put in the SAS code.)
veclines2 = function(x, y, z, c) {
require(RColorBrewer)
breaks = min(c) + (0:5 * (max(c) - min(c))/6)
# The sapply() function applies the function named in the
# second argument to each element of the vector named in the
# first argument. Since I don't need this function again,
# I write it out in within the sapply() function call and don't
# even have to name it.
ccat = sapply(c, function (r) {max(which(r >= breaks))} )
mypalette = brewer.pal(6,"Greys")
for (i in 1:(length(x)-1)) {
lines(x[i:(i+1)], y[i:(i+1)], lwd=z[i], col=mypalette[7 - ccat[i]])
}
}
Reader JanOosting pointed out that the segments() function will do exactly what my for-looped lines() function does. The segments() function requires “to” and “from” x and y locations. Here’s a version of the above function that uses it.
bikeseg = function(x,y,z,c) {
l=length(x)
require(RColorBrewer)
mypalette <-brewer.pal(6,"Greys")
breaks = min(c) + (0:5 * (max(c) - min(c))/6)
ccat = sapply(c,function (r) {max(which(r >= breaks))})
segments(x0 = x[1:(l-1)], y0 = y[1:(l-1)],
x1 = x[2:l], y1 = y[2:l],lwd=z,col=mypalette[7-ccat])
}
Then I call the function, after creating an empty plot and making the RColorBrewer library available. Finally, I add the points and some street names (see sections 5.2.3 and 5.2.11), to aid in interpretation. The way R draws the names live, as you enter commands, makes it much easier to place and rotate the names than in SAS, where you would have to re-make the annotate data set and run it to see the results.
plot(Longitude, Latitude, type="n")
library(RColorBrewer)
bikeseg(Longitude, Latitude, Speed..miles.h./2, Elevation..feet.)
points(Longitude, Latitude, cex=2, pch='.')
text(-72.495, 42.34, "Station Road")
text(-72.5035, 42.33, "Middle Street", srt=90)
text(-72.5035, 42.36, "Bike Path", srt=-39)
text(-72.54, 42.36, "Bike Path", srt=15)
text(-72.55, 42.347, "Maple", srt=107)
text(-72.54, 42.338, "Pomeroy")
text(-72.515, 42.331, "Potwine")
In the final image, thicker lines show greater speed, darker colors indicate lower elevations, and the dots are where the GPS connects. My best speed is achieved along Station Road, which is a long, straight drop. There aren’t enough categories of color to explain most of the other variations in speed, although you can see me get slower as I near Middle Street along Potwine. The dots are fairly evenly spaced except along the bike path, where there is a lot of tree cover in some spots. This causes the long straight pieces near the top of the plot. Actually, since the phone is also sitting in my pocket as I ride, so I’m more pleased than disappointed with the perfomance of the iPhone and Cyclemeter.
| This post was kindly contributed by SAS and R - go there to comment and to read the full post. |