Skip to content

Creating smooth color regions with panel.contourplot()

10 messages · David Carslaw, Deepayan Sarkar, Waichler, Scott R

#
When I use panel.contourplot() with filled color regions, the coloring
follows the stair-step edge of the underlying grid instead the smooth
contour lines themselves.  How can I get the latter behavior?  I would
guess there is a much simpler way than manually creating polygons with
contourLines(), especially since a contour interval/region can have
holes inside it (different contour intervals).

Thanks,
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler at pnl.gov
1 day later
#
On 9/15/08, Waichler, Scott R <Scott.Waichler at pnl.gov> wrote:
Manually creating polygons with contourLines will not really help
because (1) as you noted, there will be holes, and (2) contours that
cross edges will be open. The only real solution is to color each
rectangle individually, and that would be very inefficient in R code
(grid does not have a C-level interface).

The good news is that filled.contour() already does this efficiently,
and you can use that through the gridBase package:


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(do.breaks(cpl$xlim, nrow(z) - 1)),
                            as.double(do.breaks(cpl$ylim, ncol(z) - 1)),
                            z, as.double(at), col = col))
}

plot.new()

levelplot(volcano, panel = panel.filledcontour,
          col.regions = terrain.colors,
          cuts = 25)

-Deepayan
#
I think this is a very useful function that I imagine has wide appeal -
thanks.  Using the code below produces the plot OK but when I try and
copy/save it (as a metafile) I receive the following error and an hourglass:

Error: invalid graphics state

[using XP, 2.72, lattice 0.17.13]

Regards,
David
Deepayan Sarkar wrote:
-----
Institute for Transport Studies
University of Leeds
#
On Wed, Sep 17, 2008 at 1:12 PM, David Carslaw
<d.c.carslaw at its.leeds.ac.uk> wrote:
Yes, that's the drawback of mixing grid and base graphics (even
resizing the screen device will not work). However, if you explicitly
start a device and plot to it, instead of copying from a screen
device, that should work fine.

-Deepayan
#
On Wed, Sep 17, 2008 at 1:25 PM, Deepayan Sarkar
<deepayan.sarkar at gmail.com> wrote:
That isn't completely accurate. To ensure that the plot doesn't start
a new page, you need something like

pdf()
plot.new()
levelplot(volcano, panel = panel.filledcontour,
          col.regions = terrain.colors,
          cuts = 25,
          plot.args = list(newpage = FALSE))
dev.off()

-Deepayan
#
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,
          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()
#
On Thu, Sep 18, 2008 at 10:23 AM, Waichler, Scott R
<Scott.Waichler at pnl.gov> wrote:
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
#
The contour lines do not exactly line up with some of the color breaks.
I see this in the volcano example as well as the figure I'm trying to
make.  

--Scott Waichler
#
On Thu, Sep 18, 2008 at 1:36 PM, Waichler, Scott R
<Scott.Waichler at pnl.gov> wrote:
Does the contour() solution work? If not, there's not really much that
can be done.

-Deepayan
#
Yes, I figured out how to add the naming of arguments.  This works
great:

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, 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), add=T,
          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()

--Scott Waichler