[R] lattice custom axis function -- right side margins

[Timothy W. Hilton]

Dear R help list,

I have a plot with two different vertical scales that I want to display
on either side of the plot.  It's quite similar to the
Fahrenheit-Centigrade example in the examples section of the
documentation for axis.default.  

The right-side axis is clipped off, though, and I haven't been able to
figure out anything with viewport() and clipping or trellis.par.set to
fix that...  Any help greatly appreciated!  Minimal example below.

I would also like to add a label to the right-side vertical axis similar
to the "sill..." label on the left.  Bonus points if anyone can throw
that in...

Many thanks,
Tim

--

Timothy W. Hilton
PhD Candidate, Department of Meteorology
The Pennsylvania State University
503 Walker Building, University Park, PA   16802
hilton at meteo.psu.edu

--------------------------------------------------
code to produce the plot with right-side labels clipped off



example_data <-
structure(list(year = structure(c(4L, 2L, 2L, 7L, 2L, 2L, 4L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 4L, 4L, 2L, 7L, 3L, 2L, 5L), .Label = c("2000", 
"2001", "2002", "2003", "2004", "2005", "2006"), class = "factor"), 
    var = structure(c(2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 2L, 
    2L, 2L, 2L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 1L), .Label = c("NEE.model", 
    "NEE.model.res", "NEE.obs"), class = "factor"), par_set = structure(c(6L, 
    1L, 3L, 6L, 7L, 2L, 7L, 5L, 1L, 6L, 3L, 4L, 1L, 1L, 4L, 1L, 
    7L, 5L, 9L, 9L, 9L), .Label = c("all.all", "all.ann", "all.mon", 
    "pft.all", "pft.ann", "pft.mon", "site.all", "site.ann", 
    "site.mon"), class = "factor"), sigmasq = c(11430.2595455547, 
    12118.5387166954, 12982.4722525337, 16366.3059675243, 16650.2206047512, 
    19730.2121989498, 19958.3416187217, 20491.4117984889, 20647.8829877428, 
    21389.0300281264, 21413.7674128747, 21445.7255788782, 22002.8026436862, 
    22042.9802472953, 22201.0461487030, 22340.9959465200, 24782.8974616218, 
    27207.1283451608, 59450.6758048182, 94725.119215293, 694716.769010273
    )), .Names = c("year", "var", "par_set", "sigmasq"), row.names = c(94L, 
8L, 20L, 43L, 44L, 68L, 100L, 86L, 2L, 92L, 74L, 80L, 62L, 1L, 
82L, 64L, 98L, 37L, 57L, 56L, 59L), class = "data.frame")

sigsq2sig <- function(sigmasq) sqrt(2 * sigmasq)
sig2sigsq <- function(sig) 0.5 * (sig)^2

# axis method to add a std deviation axis to the right side of a sill plot 
axis.sigmasq <- function(side, ...) {
  switch(side,
         left = {
           ylim <- current.panel.limits()$ylim
           pretty_sigmasq <- pretty(ylim)
           panel.axis(side = side, outside = TRUE,
                      at = pretty_sigmasq, labels = pretty_sigmasq)
         },
         right = {
           ylim <- current.panel.limits()$ylim
           pretty_sigmasq <- pretty(ylim)
           pos_sigmasq <- pretty_sigmasq[pretty_sigmasq >= 0]
           pretty_sigma <- pretty(sigsq2sig(pos_sigmasq))
           panel.axis(side = side, outside = TRUE,
                      at = sig2sigsq(pretty_sigma),
                      labels = pretty_sigma)
         },
         axis.default(side = side, ...))
}

my_plot <- function(best.fits, ...) {

  y.main.label <- expression(sill~group("[", group("(", mu*mol~s^-1~m^-2 ,")")^2, "]"))
  #plot the parameter values, one per year
  plt <- xyplot(sigmasq~interaction(var, par_set),
                data=best.fits,
                groups=year,
                axis = axis.sigmasq,
                scales=list(x=list(rot=45)),
                xlab=list(label="NEE measure"),
                ylab=list(label=y.main.label),
                ...)
  return(plt)
}

print(my_plot(example_data))