Wednesday, August 31, 2011

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

8 comments:

  1. I'm just curious, couldn't you do it with RasterVis[1] package ?

    [1] http://rastervis.r-forge.r-project.org/

    ReplyDelete
  2. I'm not familiar with RasterVis, but it looks like a nice package with some interesting plotting functions - I will have to look into it. Also, I just came acros the function image.plot() from the "fields" packages, which seems to do a similar thing as I showed here.

    ReplyDelete
  3. Thanks ! very useful

    ReplyDelete
  4. Thank you for sharing this good information with us all. I would like to first thank you for helping me on the stackoverflow to create a contour for an estuary. I have a quick question about that problem. I added the shapefile to the plot but now I want to clip the contour by the shapefile and only display the part which is bounded by the polygon. Is that possible ? The label of legend is also not right in my case. It is showing in the order of 100,000. Any suggestions ? Thank you.


    Jdbaba

    ReplyDelete
  5. Hi Bishwamitra - I have added a comment to the stackoverflow Q.

    ReplyDelete
  6. for the first script, how do you change the colours from red to for example dark blue?

    ReplyDelete
    Replies
    1. Hi Guillaume - The first plot uses pal.1 - so you would need to change the colors used in that palette; e.g. pal.1=colorRampPalette(c("black", "blue4", "yellow"), space="rgb")

      Delete
  7. Great Post. I have a minor comment:

    I noticed that your approach for identifying the highest point on the volcano is not entirely accurate:

    > volcano[highest %% dim(volcano)[1], highest %/% dim(volcano)[1]]
    [1] 194

    I have used a similar command in the past that I find more intuitive and that is also more accurate:

    > volcano[which(volcano == max(volcano), arr.ind = TRUE)[1], which(volcano == max(volcano), arr.ind = TRUE)[2]]
    [1] 195

    Zeno

    ReplyDelete