[R] overlaid filled contour plots: solved

[Michael Friendly]

MerryXmas() to R-help!

No one answered this query, but I found an easy solution that others may 
find useful, using a wrapper for contour() that calls contourLines()
and then replots the contours using polygon().  It could use a little 
tuning (easier way to specify fill.alpha) to make it a general extension
of contour() or perhaps incorporate the additional fill.* arguments 
directly into contour, where they probably belong.

# Add facility to fill contours in contour() and to return
#  the contourLines()

#  All arguments of contour.default commented out can be
#    passed to contour() as ...
#  fill.col is a vector of fill colors corresponding to levels.
#    if not supplied, calculate a set of transparent colors
#    based on col and fill.alpha
#  TODO:  add logic for fill.col=NULL to nothing extra, other
#    than return contourLines
#  TODO:  make fill.alpha accept a more standard 0-1 numeric

contourf <- function(
		x = seq(0, 1, length.out = nrow(z)),
     y = seq(0, 1, length.out = ncol(z)),
     z,
     nlevels = 10,
     levels = pretty(zlim, nlevels),
#    labels = NULL,
#    xlim = range(x, finite = TRUE),
#    ylim = range(y, finite = TRUE),
     zlim = range(z, finite = TRUE),
#    labcex = 0.6, drawlabels = TRUE, method = "flattest",
#    vfont, axes = TRUE, frame.plot = axes,
     col = par("fg"),
     fill.col,
     fill.alpha = "80",   # hex value for alpha transparency
#    lty = par("lty"), lwd = par("lwd"),
     add = FALSE, ...) {

	contour(x,y,z, nlevels=nlevels, levels=levels, zlim=zlim, col=col, 
add=add, ...)
	line.list <- contourLines(x, y, z, nlevels=nlevels, levels=levels)
	# contourLines returns a list of lists, each with components
	# 'level', 'x', 'y'

	if(missing(fill.col)) {
		colramp = colorRampPalette(c("white", col))
		fill.col <- paste(colramp(nlevels+1), fill.alpha, sep="")
	}
	for (i in seq_along(line.list)) {
		polygon(line.list[[i]][2:3], col=fill.col[i], border=NA)
	}
	invisible(line.list)
}


best,
-Michael



On 12/22/2011 9:54 AM, Michael Friendly wrote:
> I'm trying to make a set of contour plots of bivariate kernel density
> estimates, showing three such plots overlaid,
> similar to this plot
>
> http://euclid.psych.yorku.ca/SCS/Private/Test/ridge-boot2.pdf
>
> except that I would like to have the contours *filled* (using
> transparent colors). To make this reproducible, I've
> saved the results of KernSmooth::bkde2D() in the following file:
>
> http://euclid.psych.yorku.ca/SCS/Private/Test/bkde2D-boot.RData
>
>  From this, the essence of my plot can be produced by the following:
>
> # show calls to bkde2D()
> #library(KernSmooth)
> #dest0 <- bkde2D(lbt[,wh], bandwidth=.5, range.x=list(c(-16, 8), c(-5,5)))
> #dest1 <- bkde2D(lbt[,1+wh], bandwidth=.5, range.x=list(c(-16, 8),
> c(-5,5)))
> #dest2 <- bkde2D(lbt[,2+wh], bandwidth=.5, range.x=list(c(-16, 8),
> c(-5,5)))
> #save(dest0, dest1, dest2, cm0, cm1, cm2, file="bkde2D-boot.RData")
>
> con <-
> url("http://euclid.psych.yorku.ca/SCS/Private/Test/bkde2D-boot.RData")
> ## print the value to see what objects were saved.
> print(load(con))
> close(con)
>
> col <- c("black", rainbow(5, start=.6, end=.1))
> clr <- col[c(1,3,6)]
>
> contour(dest0$x1, dest0$x2, dest0$fhat, nlevels=4,
> xlim=c(-8,4), ylim=c(-3,0), cex.lab=1.25,
> xlab="GNP", ylab="Unemployed" )
> contour(dest1$x1, dest1$x2, dest1$fhat, nlevels=4, add=TRUE, col=clr[2])
> contour(dest2$x1, dest2$x2, dest2$fhat, nlevels=4, add=TRUE, col=clr[3])
> lines(rbind(cm0, cm1, cm2))
> points(rbind(cm0, cm1, cm2), col=clr, pch=16, cex=1.5)
>
> I've tried filled.contour(), as below, except that (a) it doesn't allow
> add=TRUE, and (b) has no way to suppress the legend.
>
> # try filled contours, using transparent colors
> colramp0 = colorRampPalette(c("white", clr[1]))
> col <- paste(colramp0(5), "80", sep="") # make transparent
> filled.contour(dest0$x1, dest0$x2, dest0$fhat, nlevels=4,
> xlim=c(-8,4), ylim=c(-3,0), cex.lab=1.25, col=col,
> xlab="GNP", ylab="Unemployed" )
>
> colramp1 = colorRampPalette(c("white", clr[2]))
> col <- paste(colramp1(5), "80", sep="") # make transparent
> filled.contour(dest1$x1, dest1$x2, dest1$fhat, nlevels=4,
> xlim=c(-8,4), ylim=c(-3,0), cex.lab=1.25, col=col,
> xlab="GNP", ylab="Unemployed" )
>
> colramp2 = colorRampPalette(c("white", clr[3]))
> col <- paste(colramp2(5), "80", sep="") # make transparent
> filled.contour(dest2$x1, dest2$x2, dest2$fhat, nlevels=4,
> xlim=c(-8,4), ylim=c(-3,0), cex.lab=1.25, col=col,
> xlab="GNP", ylab="Unemployed" )
>
>
> Note that the structure of the result of bkde2D() is a list of the
> following form, so I think I'm stuck with
> plotting the 3 density estimates separately as I do above.
>
>  > str(dest0)
> List of 3
> $ x1 : num [1:51] -16 -15.5 -15 -14.6 -14.1 ...
> $ x2 : num [1:51] -5 -4.8 -4.6 -4.4 -4.2 -4 -3.8 -3.6 -3.4 -3.2 ...
> $ fhat: num [1:51, 1:51] 2.08e-04 1.60e-04 5.88e-05 1.08e-05 1.20e-06 ...
>  >
>
> Any ideas?
>
> -Michael
>