Adding a scale to an image plot

[NOTE: new version of the image.scale function can be found here: http://menugget.blogspot.de/2013/12/new-version-of-imagescale-function.html.]

Here's a function that allows you to add a color scale legend to an image plot (or probably any plot needing a z-level scale). I found myself having to program this over and over again, and just decided to make a plotting function for future use. While I really like the look of levelplot(), the modular aspect of image() makes it much more handy to combine with other plotting commands or overlays.
For example, as far as I can tell, the simple addition of the triangle symbol to mark the highest point in the above map of Maunga Whau volcano is not possible with levelplot.
After adding this symbol, the function below - image.scale() - was used to add the accompanying color scale to another area of the device.



The function...

#This function creates a color scale for use with e.g. the image()
#function. Input parameters should be consistent with those
#used in the corresponding image plot. The "horiz" argument
#defines whether the scale is horizonal(=TRUE) or vertical(=FALSE).
#Depending on the orientation, x- or y-limits may be defined that
#are different from the z-limits and will reduce the range of
#colors displayed.
 
image.scale <- function(z, zlim, col = heat.colors(12),
breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
 if(!missing(breaks)){
  if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
 }
 if(missing(breaks) & !missing(zlim)){
  breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1)) 
 }
 if(missing(breaks) & missing(zlim)){
  zlim <- range(z, na.rm=TRUE)
  zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
  zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
  breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
 }
 poly <- vector(mode="list", length(col))
 for(i in seq(poly)){
  poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
 }
 xaxt <- ifelse(horiz, "s", "n")
 yaxt <- ifelse(horiz, "n", "s")
 if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
 if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
 if(missing(xlim)) xlim=XLIM
 if(missing(ylim)) ylim=YLIM
 plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i", ...)  
 for(i in seq(poly)){
  if(horiz){
   polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
  }
  if(!horiz){
   polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
  }
 }
}

Created by Pretty R at inside-R.org




The script used to create the image above...

source("image.scale.R")
 
pal.1=colorRampPalette(c("black", "red", "yellow"), space="rgb")
pal.2=colorRampPalette(c("black", "blue", "cyan"), space="rgb")
 
png("volcano_w_scale.png", width=7, height=4, units="in", res=200)
layout(matrix(c(1,2,3,0,4,0), nrow=2, ncol=3), widths=c(4,4,1), heights=c(4,1))
layout.show(4)
 
#1st image
breaks <- seq(min(volcano), max(volcano),length.out=100)
par(mar=c(1,1,1,1))
image(seq(dim(volcano)[1]), seq(dim(volcano)[2]), volcano, 
col=pal.1(length(breaks)-1), breaks=breaks, xaxt="n", yaxt="n", ylab="", xlab="")
#Add additional graphics
highest <- which.max(volcano)
points(highest %% dim(volcano)[1], highest %/% dim(volcano)[1], 
pch=2, lwd=2, cex=2,col="blue")
#Add scale
par(mar=c(3,1,1,1))
image.scale(volcano, col=pal.1(length(breaks)-1), breaks=breaks, horiz=TRUE)
box()
 
 
#2nd image
breaks <- c(0,100, 150, 170, 190, 200)
par(mar=c(1,1,1,1))
image(seq(dim(volcano)[1]), seq(dim(volcano)[2]), volcano, 
col=pal.2(length(breaks)-1), breaks=breaks, xaxt="n", yaxt="n", ylab="", xlab="")
#Add additional graphics
highest <- which.max(volcano)
points(highest %% dim(volcano)[1], highest %/% dim(volcano)[1], 
pch=2, lwd=2, cex=2,col="red")
#Add scale
par(mar=c(1,1,1,3))
image.scale(volcano, col=pal.2(length(breaks)-1), breaks=breaks, horiz=FALSE, yaxt="n")
axis(4,at=breaks, las=2)
box()
 
dev.off()

Created by Pretty R at inside-R.org