0
votes

I'm trying to produce a very specific graphical output, and can't find a package that can do all of the things I'm looking for. Primarily, I need to: (1) Produce multiple plots (facet_wrap in ggplot works great) (2) Have 2 y-axes with different scales (lattice seems the best for this) (3) Produce a graph similar to the output from "plot" using dplR

I have a data frame with multiple unique Sites. I'd like to plot each site's average annual value as a line graph, with a second line graph on the other y-axis representing sample depth. I'd like to produce a separate plot for each site.

Below is some reproducible data(small subset), and what I've come up with so far:

> sites <- structure(list(Year = c("2009", "2009", "2009", "2009", "2009", 
"2009", "2009", "2009", "2009", "2009", "2009", "2009", "2009", 
"2009", "2009", "2009", "2009", "2009", "2009", "2009", "2009", 
"2010", "2010", "2010", "2010", "2010", "2010", "2010", "2010", 
"2010", "2010", "2010", "2010", "2010", "2010", "2010", "2010", 
"2010", "2010", "2010", "2010", "2010", "2010", "2010", "2010", 
"2011", "2011", "2011", "2011", "2011", "2011", "2011", "2011", 
"2011", "2011", "2011", "2011", "2011", "2011", "2011", "2011", 
"2011", "2011", "2011", "2011", "2011", "2011", "2012", "2012", 
"2012"), plot = c("FA6", "FA7", "FK1", "FK2", "FK3", "FO1", "FO2", 
"FO6", "GA2", "GA4", "GR1", "GR2", "HE2", "HE3", "LY1", "LY3", 
"LY8", "NM2", "NM3", "TH3", "TH5", "BR1", "BR8", "FA5", "FA6", 
"FA7", "FK1", "FK2", "FK3", "FO1", "FO2", "FO6", "GA2", "GA4", 
"GR1", "GR2", "HE2", "HE3", "LY1", "LY3", "LY8", "NM2", "NM3", 
"TH3", "TH5", "FA5", "FA6", "FA7", "FK1", "FK2", "FK3", "FO1", 
"FO2", "FO6", "GA2", "GA4", "GR1", "GR2", "HE2", "HE3", "LY1", 
"LY3", "LY8", "NM2", "NM3", "TH3", "TH5", "HE2", "HE3", "TH5"
), AvgRW = c(0.628666666666667, 0.485027777777778, 0.479269230769231, 
0.826875, 0.633269230769231, 1.01830769230769, 1.34580555555556, 
1.13061764705882, 0.422375, 1.377625, 0.535375, 0.366384615384615, 
0.493119047619048, 0.300777777777778, 0.971923076923077, 1.02302941176471, 
1.47245833333333, 1.00654166666667, 0.56425, 1.66342857142857, 
1.28477586206897, 0.860666666666667, 2.10155130769231, 1.74626923076923, 
0.616148148148148, 0.42775, 0.402576923076923, 0.859333333333333, 
0.608961538461538, 1.28303846153846, 1.84344444444444, 1.52214705882353, 
0.425546875, 1.66179166666667, 0.647208333333333, 0.390461538461538, 
0.565892857142857, 0.237388888888889, 1.60419230769231, 1.16611764705882, 
1.95329166666667, 1.18795833333333, 0.655928571428571, 2.009, 
1.36198275862069, 2.61165384615385, 0.873296296296296, 0.596, 
0.485884615384615, 1.13633333333333, 0.684461538461538, 1.30946153846154, 
1.69747222222222, 1.64197058823529, 0.40740625, 1.40716666666667, 
0.641625, 0.428576923076923, 0.729011904761905, 0.376222222222222, 
1.52984615384615, 1.15317647058824, 1.66183333333333, 1.17904166666667, 
0.604857142857143, 1.57425, 1.55772222222222, 0.7315, 0.119, 
1.125875), SampDepth = c(27L, 18L, 13L, 12L, 13L, 13L, 18L, 17L, 
32L, 12L, 12L, 13L, 14L, 18L, 13L, 17L, 12L, 12L, 14L, 14L, 29L, 
21L, 13L, 13L, 27L, 18L, 13L, 12L, 13L, 13L, 18L, 17L, 32L, 12L, 
12L, 13L, 14L, 18L, 13L, 17L, 12L, 12L, 14L, 14L, 29L, 13L, 27L, 
18L, 13L, 12L, 13L, 13L, 18L, 17L, 32L, 12L, 12L, 13L, 14L, 18L, 
13L, 17L, 12L, 12L, 14L, 4L, 9L, 2L, 1L, 8L)), .Names = c("Year", 
"plot", "AvgRW", "SampDepth"), row.names = 2330:2399, class = "data.frame")

I can get a great layout using ggplot and ggplus, but can't seem to add second y-axis:

Plot1 <- ggplot(sites,aes(x=Year, y=AvgRW, group=plot)) + 
    theme(axis.text=element_text(size=4), 
          panel.grid.minor=element_blank(), 
          panel.grid.major=element_line(colour = "gray",size=0.5), 
          strip.text=element_text(size=8)) + 
    geom_line()

facet_multiple(plot = Plot1, facets = 'plot', ncol = 1, nrow = 1, scales="free")

Sample graph

Then I tried doing a doubleYplot with lattice and latticeExtra, but I can't get the lines to appear:

Plot2a <- xyplot(AvgRW ~ Year, data=sites, groups=plot,
              xlab=list("Year", fontsize = 12),
              ylab = list("Avg RW (mm)", fontsize = 12),
              ylab.right = list("Sample Depth", fontsize = 12),
              par.settings = simpleTheme(col = 1),
              type = c("l","g"),
              lty=1,
              scales=list(x=list(rot=90,tick.number=25,
                                 cex=1,axs="r")))

> Plot2b <- xyplot(SampDepth ~ Year,data=sites, groups=plot, type = "o",col="black",
              lty=9)

> doubleYScale(Plot2a, Plot2b)

My ultimate goal is to produce a graph that looks similar to the dplR output, for every site in my data.frame:

Intended output

Any ideas on how to combine my three goals would be greatly appreciated.

1
BTW, you can have two y-axes labels in GGplot - Jthorpe
Thanks. I know ggplot can do 2 axis labels, but I don't know how to graph 2 lines with the same x axis but different y values. - KKL234
The link in my comment illustrates how to have the same scale on the x-axis and separate scales on the left and right y-axes. The use of ncol = 1, nrow = 1, scales="free" produces one figure per slide, which could be accomplished via pdf(blah,blah,blah); for(p in unique(sites$plot) show(Plot1 %+% sites[sites$plot == p]); dev.off(). If you're trying to add a second line on the same y-scale, use Plot1 + geom_smooth(blah,blah,blah) - Jthorpe
Thanks for the suggestion. I looked at the link and applied the code to my own data, however it's producing all plots on 1 graph. I tried adding the for loop in, but loops aren't my specialty: I can't get the loop to run correctly. If you could help me figure out how to apply the loop, I'd appreciate it. - KKL234

1 Answers

1
votes

ggplot2 does not support a secondary y-axis because his creator believe that it can be misleading. However, you can do it with base R. Here's a loop that solves your problems: multiple plots in a loop, secondary y-axis, similar look to the attached plot. I assumed that you are on Windows and included the code to save every chart in your /temp directory using savePlot.

for (i in sites$plot){                          #start loop
par(mar=c(5.1,4.1,4.1,5.1))                     #increase right margin
plot(sites[sites$plot==i,c(1,3)],xaxt="n", type="l", ylab="Avg RW")
title(main=i,line=2.5)                              #add title
axis(1,at=as.integer(sites$Year),labels=sites$Year) #bottom axis
axis(3,at=as.integer(sites$Year),labels=sites$Year) #top axis
par(new = TRUE)                                     #overlay for secondary y
plot(sites[sites$plot==i,c(1,4)],xaxt="n",yaxt="n", ylab="",type="l", col="red")
axis(4)                                             #add secondary y axis
mtext("Sample Depth", side = 4, line=2)             #add secondary y label
savePlot(filename = paste0("c:/temp/",i, ".png"), type = "png") # save
}

enter image description here

UPDATE If you are not on Windows, you cannot use savePlot. Use the normal png function. You will not see the charts on screen but they will be saved. Make sure to change the path to where you want the charts saved. I am using c:/temp/ in the following code, which will not work on non-windows machines.

for (i in sites$plot){                          #start loop
png(filename = paste0("c:/temp/",i, ".png")) # save
par(mar=c(5.1,4.1,4.1,5.1))                     #increase right margin
plot(sites[sites$plot==i,c(1,3)],xaxt="n", type="l", ylab="Avg RW")
title(main=i,line=2.5)                              #add title
axis(1,at=as.integer(sites$Year),labels=sites$Year) #bottom axis
axis(3,at=as.integer(sites$Year),labels=sites$Year) #top axis
par(new = TRUE)                                     #overlay for secondary y
plot(sites[sites$plot==i,c(1,4)],xaxt="n",yaxt="n", ylab="",type="l", col="red")
axis(4)                                             #add secondary y axis
mtext("Sample Depth", side = 4, line=2)             #add secondary y label
dev.off()
}