drawImage <- function(m, xtimes) {
        # for grey background of NaN values
	grey <- mat.or.vec(nrow(m),ncol(m))
	for (i in 1:nrow(m)) { for (j in 1:ncol(m)) { 
	  if (!is.na(m[i,j])) {
	    if (m[i,j] > 3.0) { m[i,j] = 3.0 }
	    if (m[i,j] < -3.0) { m[i,j] = -3.0 }
	    grey[i,j] = NA
	  } else { 
	    grey[i,j] = 0 
	  }
        
	xr = mat.or.vec(nrow(xtimes),1)
	for(i in 1:nrow(xtimes)) {
          xr[i]=xtimes[i,1]
        }
        cx <- data.matrix(m)
        colnames(cx) = c("Antarctica","South\nAmerica","Australasia","N.America\n(pollen)","N.America\n(tree rings)","Asia","Europe","Arctic")

        # palette for temperature deviations, using 2008 Roger D. Pengs color.palette function (GNU)
        # see 
	tempDevPal = color.palette(c("dark blue", "cyan", "white", "yellow", "dark red"), n.steps.between=c(10,3,3,10))
        side2 <- 0.1
	# row name width, N.America more dense
	cols <- c(0.05,0.1,0.15,0.2,0.225,0.275,0.325,0.375)
        par(las = 2, cex.axis = 0.5, lab=c(20,8,7))
        side2 <- max(side2, max(strwidth(colnames(cx), "inches")) + 0.1)
        par(mai = c(3.5, side2, 0.2, 0.1))
        # two columns, left the actual plot, right the color scale
        layout(matrix(data=c(1,2), nrow=1, ncol=2), widths=c(4,1), heights=c(1,1))

        # draw temperature dev.
        image(unclass(xr), cols, cx, col = tempDevPal(61),
              xlim = c(1,2000), zlim = c(-3.0,3.0), xaxt = "n", yaxt = "n", ylab = "", xlab = "") # range(xr)
        # draw grey background for NaN values
        image(unclass(xr), cols, grey, col = "light grey",
              xaxt = "n", yaxt = "n", ylab = "", xlab = "", add=TRUE)
        axis(2, at = cols, colnames(cx))
        Axis(xr, side = 1)
        box()

	# draw color scale
        ColorLevels <- seq(-3.0, 3.0, length=61)
        par(mai = c(3.5, 0.6, 0.2, 0.8))
	par(mfcol=c(1,4),new=TRUE)
	image(1, ColorLevels, matrix(data=ColorLevels, ncol=length(ColorLevels),nrow=1), col=tempDevPal(61), xlab="", ylab="", xaxt = "n")
}