[R-sig-Geo] contour plots with lattice

[Martin Ivanov]

Dear R users,

I need to overlay the plot of one spatial variable with the contours of 
another spatial variable. And I want to do this with lattice.
Here is an example with the meuse.grid data set that comes with the R 
package sp.
First of all, in the classical function "contour" labels are drawn so 
that each contour line is broken at the place of the label,
so that the label does not overlap the line; besides care is taken that 
the labels look good as labels are plotted only for
those closed contours that are large enough for the label - if the 
closed contour is too small, no label is drawn:

library(sp);
library(lattice);
library(latticeExtra);
library(gridBase);
data(meuse.grid);
coordinates(meuse.grid) <- c("x", "y");
pixDF <- as(meuse.grid, "SpatialPixelsDataFrame")["dist"];
zlevs.conts <- zlevs.fill <- seq(0, 1, 0.1);
Z <- matrix(NA, nrow = pixDF at grid@cells.dim[1], ncol = 
pixDF at grid@cells.dim[2]); # convert the SpatialPixelsDataFrame to a matrix
Z[pixDF at grid.index] <- pixDF at data$dist;
Z <- Z[, ncol(Z):1];
contour(Z, levels=zlevs.conts);

In the trellis contourplot the labels are placed beside or over the 
line, controlled by the label.style parameter:
levelplot(Z, contour=TRUE, labels=TRUE, at=zlevs.conts, region=FALSE);
Labels are drawn at every single contour, even at the closed ones that 
are too small and the map looks nasty.

I do not like this and I could not find an option that would place the 
labels just as in the classical contour: breaking the contour line to 
make place for
the label, and with a smart consideration where to place a label and 
where not. I made several posting in the R and r-sig-geo forum on this 
issue,
and I was lucky to receive generous help from Dr. Tim Appelhans from the 
r-sig-geo forum. He wrote a function for me that combines the 
traditional R contour with grid graphics. With his permission am posting 
his function here:

panel.filledcontour <- function(x, y, z, subscripts, at, fill.cont = T,
                                 col.regions,
                                 contours = T,
                                 col = col.regions(length(zlevs.fill)),
                                 ...)
{
   stopifnot(require("gridBase"))
   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()
   plot.window(xlim = cpl$xlim, ylim = cpl$ylim,
               log = "", xaxs = "i", yaxs = "i")
   # paint the color contour regions
   if (isTRUE(fill.cont))
     .Internal(filledcontour(as.double(do.breaks(cpl$xlim,
                                                 nrow(z) - 1)),
                             as.double(do.breaks(cpl$ylim,
                                                 ncol(z) - 1)),
                             z, levels = as.double(zlevs.fill),
                             col = col))
   else NULL

   #add contour lines
   if (isTRUE(contours))
     contour(as.double(do.breaks(cpl$xlim, nrow(z) - 1)),
             as.double(do.breaks(cpl$ylim, ncol(z) - 1)),
             z, levels = as.double(zlevs.conts),
             add=T, col = "grey10", # color of the lines
             drawlabels = T, ...
              # add labels or not
     )
   else NULL

}

Practically this function solves my problem:

plot1 <- spplot(pixDF, at=zlevs.conts);
plot2 <- levelplot(Z, panel = function(fill.cont, contours, ...) {
           panel.filledcontour(fill.cont = F, contours = T, ...) }, 
col.regions = colorRampPalette(c("snow", "lightblue")),
           plot.args = list(newpage = FALSE))

plot.new();print( plot1 + as.layer(plot2));

My data are all spatial, so they are available as spatial objects. 
However, spplot does not seem to be able to overlay
a plot with a contourplot of another spatial variable, it can only 
overlay spatial points, spatial lines, spatial polygons and text (at least
according to the spplot documentation: see the sp.layout argument). That 
is why one has to use levelplot for the second field. But levelplot
does not recognise spatial objects, so one has to convert the spatial 
object to a matrix. So far so good. But if spplot is given xlim and ylim
arguments:

bbx <- bbox(pixDF);
xlim <- c(.999*bbx[1, 1], 1.001*bbx[1, 2]);
ylim <- c(.999*bbx[2, 1], 1.001*bbx[2, 2]);
plot11 <- spplot(pixDF, at=zlevs.conts, xlim=xlim, ylim=ylim);

plot2 <- levelplot(Z, panel = function(fill.cont, contours, ...) {
           panel.filledcontour(fill.cont = F, contours = T, ...) }, 
col.regions = colorRampPalette(c("snow", "lightblue")),
           plot.args = list(newpage = FALSE))

plot.new();print(plot11 + as.layer(plot2))


things go wrong: the second field is plotted in the old scales of 
course, while spplot plots the first field according to the
newly given scales. This can only be solved by converting the matrix of 
the second field appropriately so that it matches the
new scales of the first field. So my problem is in principle over, but 
nevertheless I have 3 questions:

1. Is there a simpler way to achieve the final result I need with 
trellis/grid graphics? Or possibly some not necessarily simpler way, but 
a way
not relying on traditional R graphics, as so long as I know combinig 
traditional and grid graphics does have drawbacks.
2. Is it possible to get the final result with one or two calls to 
spplot, that is, working
only with spatial objects, without having to transform one of the 
objects to a matrix?
3. How can I see the code of spplot? I tried spplot, sp::spplot, 
sp:::spplot, but neither of the commands returned the code.

I am looking forward to Your comments. Any suggestions will be appreciated.
I would again like to thank Dr. Appelhans for his generous help.


Best regards,

-- 
Dr. Martin Ivanov
Eberhard-Karls-Universität Tübingen
Mathematisch-Naturwissenschaftliche Fakultät
Fachbereich Geowissenschaften
Water & Earth System Science (WESS)
Keplerstraße 17, 72074 Tübingen, Deutschland
Tel. +4970712975023