R plot 2 axes multiple plots

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")

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:

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

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
}

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()
}