0

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.

joran
  • 169,992
  • 32
  • 429
  • 468
KKL234
  • 367
  • 1
  • 5
  • 23
  • BTW, you can have [two y-axes](http://rpubs.com/kohske/dual_axis_in_ggplot2) labels in GGplot – Jthorpe Oct 12 '15 at 22:27
  • 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 Oct 13 '15 at 02:37
  • 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 Oct 13 '15 at 05:08
  • 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 Oct 13 '15 at 20:08

1 Answers1

1

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()
}
Pierre Lapointe
  • 16,017
  • 2
  • 43
  • 56
  • Thanks! This worked great; the output graph is just what I was looking for. The only issue I had was the actual loop: when I ran it in R, only the first plot from the data.frame was graphed. I couldn't locate any of the other graphs. I applied it both to the subsetted data.frame I used above, and also to my own full data.frame, and in both cases I could only get the last graph to appear. Any ideas? – KKL234 Oct 15 '15 at 15:21
  • I should also add that I'm getting this error message: "Error in savePlot(filename = paste0("c:/temp/", i, ".png"), type = "png") : can only copy from 'windows' devices" – KKL234 Oct 15 '15 at 15:25
  • Yes, I wrote that I assumed that you were on Windows in my original answer. I added an update, which will work on non-windows machines. Make sure to change the path to where you want the charts saved (should not be c:/temp/) – Pierre Lapointe Oct 15 '15 at 15:56
  • By the way, you could only see the first chart because of the savePlot error. The loop stops when it encounters an error. – Pierre Lapointe Oct 15 '15 at 15:58
  • The fix worked great! I had to change ".png" to type="png" and then it ran exactly as needed! – KKL234 Oct 15 '15 at 16:56