Converting values to color levels

     Adding color to a plot is helpful in many situations for visualizing an additional dimension of the data. Again, I wrote the below function "val2col" in R after having coded this manually over and over in the past. It uses similar arguments as the image function in that one defines the colors to be used as well as optional break points or z-limits by which the data is binned into those colors. The typical situation where I use the function is with the mapping of climate data, yet the addition of color to an XY plot can often be easier on the eyes than adding an additional z-axis spatial dimension. In combination with the image.scale function, that I previously posted, the data can be quickly interpretted.
     As an example, gridded sea level pressure is plotted above as projected polygons using the packages maps and mapproj. Values were converted to colors with the val2col function and the image.scale function plotted a corresponding scale. For those interested in using netcdf files, the example also uses the ncdf package for reading the data files into R.

the val2col function ...

#this function converts a vector of values("z") to a vector of color
#levels. One must define the number of colors. The limits of the color
#scale("zlim") or the break points for the color changes("breaks") can 
#also be defined. when breaks and zlim are defined, breaks overrides zlim.
val2col<-function(z, zlim, col = heat.colors(12), breaks){
 if(!missing(breaks)){
  if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
 }
 if(missing(breaks) & !missing(zlim)){
  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)) 
 }
 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))
 }
 CUT <- cut(z, breaks=breaks)
 colorlevels <- col[match(CUT, levels(CUT))] # assign colors to heights for each point
 return(colorlevels)
}

Created by Pretty R at inside-R.org



and the script to produce the above figure...

#Hadley SLP monthly mean dataset
#http://www.esrl.noaa.gov/psd/gcos_wgsp/Gridded/data.hadslp2.html
#ftp://ftp.cdc.noaa.gov/Datasets.other/hadslp/slp.mnmean.nc
 
library(maps)
library(mapproj)
source("val2col.R")
source("image.scale.R")
 
#slp
library(ncdf)
nc <- open.ncdf("slp.mnmean.nc")   # opens nc file 
slp.lon <- get.var.ncdf( nc, "lon")
slp.lat <- get.var.ncdf( nc, "lat")
slp.t <- get.var.ncdf( nc, "time")
slp.raw <- get.var.ncdf(nc, "slp")
close.ncdf(nc)
slp.t <- as.Date(slp.t, origin="1800-01-01")
temp <- which(slp.lon>180)
slp.lon[temp] <- slp.lon[temp]-360
slp.grd <- expand.grid(slp.lon, slp.lat)
colnames(slp.grd) <- c("lon", "lat")
slp <- matrix(c(slp.raw), nrow=length(slp.t), ncol<-length(slp.lon)*length(slp.lat), byrow=TRUE)
row.names(slp) <- as.character(slp.t)
 
#polygons for the grids
spacing=5 # degrees
slp.poly <- vector(mode="list", dim(slp.grd)[1])
for(i in seq(slp.poly)){
 x=c(slp.grd[i,1]-spacing/2, slp.grd[i,1]+spacing/2, slp.grd[i,1]+spacing/2, slp.grd[i,1]-spacing/2)
 y=c(slp.grd[i,2]-spacing/2, slp.grd[i,2]-spacing/2, slp.grd[i,2]+spacing/2, slp.grd[i,2]+spacing/2)
 slp.poly[[i]] <- data.frame(x=x, y=y)
}
 
#figure settings
heights=c(4)
widths=c(4,1)
pal=colorRampPalette(c("cyan", "white", "magenta"))
ncol=100
res=200
YYYYMMDD <- match("1998-01-01", rownames(slp))
colorvalues <- val2col(slp[YYYYMMDD,], col=pal(ncol)) #color levels for the polygons
 
#mapproj settings
project="perspective"
YLIM=c(-90, 90)
XLIM=c(-180, 180)
orientation=c(75, -50, 0)
PAR=3
 
#figure
png(filename="map_slp.png", width = sum(widths), height = sum(heights), units="in", res=res)
par(omi=c(0.1, 0.1, 0.1, 0.1), ps=12)
layout(matrix(c(1,2),nrow=1,ncol=2,byrow=TRUE), widths = widths, heights = heights, respect=TRUE)
layout.show(2)
 
##plot of map
par(mai=c(0.2, 0.2, 0.2, 0.2))
map("world",project=project, orientation=orientation, par=PAR, ylim=YLIM, xlim=XLIM)
for(i in 1:length(slp.poly)){
 polygon(mapproject(x=slp.poly[[i]][,1], y=slp.poly[[i]][,2]), col=colorvalues[i], border=colorvalues[i], lwd=0.3)
}
map("world",project=project, orientation=orientation, par=PAR, fill=FALSE, add=TRUE, col="black", ylim=YLIM, xlim=XLIM)
map.grid(c(-180, 180, -90, 90), nx=36, ny=18, labels=FALSE, col="grey", lwd=1)
 
#add scale
par(mai=c(0.2, 0, 0.2, 0.6))
image.scale(slp[1,], col=pal(ncol), horiz=FALSE, yaxt="n")
axis(4)
mtext("mb", side=4, line=2.5)
box()
 
dev.off()

Created by Pretty R at inside-R.org