[R] Creating smooth color regions with panel.contourplot()

[Deepayan Sarkar]

On Thu, Sep 18, 2008 at 10:23 AM, Waichler, Scott R
<Scott.Waichler at pnl.gov> wrote:
> Thank you very much, Deepayan.  There is just one more feature I'd like
> to get, the ability to add the contour lines.  My revision to your code
> below prints too many lines.  What needs to be changed?
>
> --Thanks,
> Scott Waichler
> scott.waichler at pnl.gov
>
> library(gridBase)
> library(lattice)
> data(volcano)
>
> panel.filledcontour <- function(x, y, z, subscripts, at, col.regions =
> cm.colors,
>                                col = col.regions(length(at) - 1), ...)
> {
>  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
>  .Internal(filledcontour(as.double(do.breaks(cpl$xlim, nrow(z) - 1)),
>                          as.double(do.breaks(cpl$ylim, ncol(z) - 1)),
>                          z, as.double(at), col = col))
>
>  # add the contour lines--this prints too many of them.  I really want
> just
>  # the lines dividing the color regions.
>  contour(as.double(do.breaks(cpl$xlim, nrow(z) - 1)),
>                            as.double(do.breaks(cpl$ylim, ncol(z) - 1)),
>                            z, as.double(at), add=T,

You need to name arguments here. as.double(at) is being matched to
'nlevels', but you want 'levels'.

Another option is to use panel.levelplot() for the contours. E.g.,
(also with more accurate 'x' and 'y'):

panel.filledcontour <-
   function(x, y, z, subscripts,
            at,
            col.regions = cm.colors,
            col = col.regions(length(at) - 1),
            ...)
{
   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")
   .Internal(filledcontour(as.double(sort(unique(x[subscripts]))),
                           as.double(sort(unique(y[subscripts]))),
                           z, as.double(at), col = col))
   panel.contourplot(x, y, z, subscripts = subscripts, at = at,
                   region = FALSE, contour = TRUE, labels = FALSE)
}

-Deepayan

>          col = "gray", # color of the lines
>          drawlabels=F  # add labels or not
>         )
> }
>
> pdf("volcano.pdf")
> plot.new()
>
> print(levelplot(volcano, panel = panel.filledcontour,
>          col.regions = terrain.colors,
>          cuts = 10,
>          plot.args = list(newpage = FALSE)))
> dev.off()
>
>