[R] Using a background photo in lattice--title and axes missing

[Waichler, Scott R]

I am trying to use a background photo in a lattice plot.  I am using the
rimage and TeachingDemos packages to plot the photo and translate from
the photo coordinates in pixels to geographic coordinates, which is what
I want to use for plotting contours, lines, etc.  The (unrunable) code
below does give me a plot showing the photo, color contours, contour
lines, and colorkey, but not the plot title or axes.  How can I get the
title and axes to appear?

  # Define a panel function that fills the true contoured regions with
color.
  # This requires the gridBase package.  panel.contourplot() gives a
stair-step, rectangular
  # pattern caused by the underlying grid.  This custom function from
Deepayan Sarkar
  # fills the contour color regions right up to the contour lines
themselves.
  # This version uses pre-defined colors.
  panel.filledContour2 <- function(x, y, z, subscripts, at, col =
my.colors, ...) {
    stopifnot(require("gridBase"))
    stopifnot(require("rimage"))  # use for plotting a background photo
    stopifnot(require("TeachingDemos")) # use for changing from photo
coordinates (pixels) to geographic coordinates
    z <- matrix(z[subscripts], nrow = length(unique(x[subscripts])),
ncol = length(unique(y[subscripts])))
    if (!is.double(z)) storage.mode(z) <- "double"
    opar <- par(no.readonly = TRUE)
    on.exit(par(opar))
    if (panel.number() > 1) par(new = TRUE)
    par(fig = gridFIG(), omi = c(0, 0, 0, 0), mai = c(0, 0, 0, 0))
    cpl <- current.panel.limits() # set the plot window boundaries to
the current panel limits
    plot.window(xlim = cpl$xlim, ylim = cpl$ylim, log = "", xaxs = "i",
yaxs = "i")

    # Plot the background photo
    plot(photo)
    
    # Change from photo to geographic coordinate systems.
    # Reference two points in the image coordinate system to the new
coordinate system.
    # x1, y1 are in the photo coordinates; x2, y2 are the map/plotting
coordinates you want to work in.
    updateusr(x1 = photo.resize.factor * x1r, y1 = photo.resize.factor *
y1r, x2 = x2r, y2 = y2r)

    # Plot the color contours
    .Internal(filledcontour(as.double(do.breaks(cpl$xlim, nrow(z) - 1)),
                            as.double(do.breaks(cpl$ylim, ncol(z) - 1)),
                            z, levels = as.double(at), col = col))
    
    # Add contour lines
    contour(as.double(do.breaks(cpl$xlim, nrow(z) - 1)),
                            as.double(do.breaks(cpl$ylim, ncol(z) - 1)),
                            z, levels = as.double(at), col = "black",
drawlabels=T, add=T, labcex=0.8)
  }  # end of panel.filledContour2()

  plotfile <- plotfile2b
  png(file = plotfile, width=11, height=8.5, units="in", res=75,
pointsize=12)  # landscape letter
  
  my.cuts <- seq(0, 90, by=10)
  my.colors <- rev(BluetoOrange.10[1:9]) # orange to blue (warm to cold
with increasing time)
  #my.colors <- rgb(col2rgb(my.colors)/255, alpha=0.2)  # make
translucent; only works for pdf, quartz devices
  
  x.limits <- c(593600, 593950)
  x.ticks <- pretty(x.limits, n=5)
  y.limits <- c(113300, 113530)
  y.ticks <- pretty(y.limits, n=5)
  
  p$time[ind.neg] <- 0
  ind <- ind.all.times.injected
  di <- interp(x=p$x[ind], y=p$y[ind], z=p$time[ind], duplicate="mean",
linear=T, 
               xo=seq(x.limits.photo[1], x.limits.photo[2], length =
300), yo=seq(y.limits.photo[1], y.limits.photo[2], length = 300))
  ia <- which(!is.na(di$z))
  grid <- expand.grid(x=di$x, y=di$y)
  
  plot.new()
  
  # contour plot doesn't print a key by default
  print(contourplot(di$z ~ x * y, grid,
    # Plot head as filled color contours
    panel = panel.filledContour2, subscripts = 1:length(grid$x),

    at = my.cuts,
    col.regions = my.colors,
    colorkey=T,
    contour = TRUE, 
    aspect="iso", 
    as.table=T,
    scales = list(x = list(relation="same", alternating = T,
                           at =x.ticks.photo, labels = x.ticks.photo,
                           limits=x.limits.photo
                          ), 
                  y = list(relation="same", alternating = T,
                           at = y.ticks.photo, labels = y.ticks.photo,
                           rot=0,
                           limits=y.limits.photo
                          )
                  ),
    xlab=list(label="X (m)", cex=1.1),
    ylab=list(label="Y (m)", cex=1.1),
    layout=c(1,1), # ncols, nrows
    main=list(label="Particle Travel Times", cex=1.2),
    strip = F,   # set to FALSE if you don't want strip drawn in plot (a
subtitle)
    plot.args = list(newpage = FALSE)        
  )) # end print(xyplot())
  dev.off()


Thanks,
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler at pnl.gov