Adding stippling to image/contour plot

Question

have some data that I would like to add "stippling" to show where it is "important", as they do in the IPCC plots

At the moment I am really struggling with trying to do this in R.

If I make up some test data and plot it:

data <- array(runif(12*6), dim=c(12,6) )
over <- ifelse(data > 0.5, 1, 0 )
image(1:12, 1:6, data)

What I would like to finally do is over-plot some points based on the array "over" on top of the current image.

Any suggestions!??

Answer 1

This should help - I had do do a similar thing before and wrote a function that I posted here.

#required function from www.menugget.blogspot.com
matrix.poly <- function(x, y, z=mat, n=NULL){
 if(missing(z)) stop("Must define matrix 'z'")
 if(missing(n)) stop("Must define at least 1 grid location 'n'")
 if(missing(x)) x <- seq(0,1,,dim(z)[1])
 if(missing(y)) y <- seq(0,1,,dim(z)[2])
 poly <- vector(mode="list", length(n))
 for(i in seq(length(n))){
  ROW <- ((n[i]-1) %% dim(z)[1]) +1
  COL <- ((n[i]-1) %/% dim(z)[1]) +1
 
  dist.left <- (x[ROW]-x[ROW-1])/2
  dist.right <- (x[ROW+1]-x[ROW])/2
  if(ROW==1) dist.left <- dist.right
  if(ROW==dim(z)[1]) dist.right <- dist.left
 
  dist.down <- (y[COL]-y[COL-1])/2
  dist.up <- (y[COL+1]-y[COL])/2
  if(COL==1) dist.down <- dist.up
  if(COL==dim(z)[2]) dist.up <- dist.down
 
  xs <- c(x[ROW]-dist.left, x[ROW]-dist.left, x[ROW]+dist.right, x[ROW]+dist.right)
  ys <- c(y[COL]-dist.down, y[COL]+dist.up, y[COL]+dist.up, y[COL]-dist.down)
  poly[[i]] <- data.frame(x=xs, y=ys)
 }
 return(poly)
}

#make vector of grids for hatching
incl <- which(over==1)

#make polygons for each grid for hatching
polys <- matrix.poly(1:12, 1:6, z=over, n=incl)

    #plot
png("hatched_image.png")
image(1:12, 1:6, data)
for(i in seq(polys)){
    polygon(polys[[i]], density=10, angle=45, border=NA)
    polygon(polys[[i]], density=10, angle=-45, border=NA)
}
box()
dev.off()

Or, and alternative with "stipples":

png("hatched_image2.png")
image(1:12, 1:6, data)
for(i in seq(polys)){
    xran <- range(polys[[i]]$x)
    yran <- range(polys[[i]]$y)
    xs <- seq(xran[1], xran[2],,5)
    ys <- seq(yran[1], yran[2],,5)
    grd <- expand.grid(xs,ys)
    points(grd, pch=19, cex=0.5)
}
box()
dev.off()

Update:

In (very late) response to Paul Hiemstra's comment, here are two more examples with a matrix of higher resolution. The hatching maintains a nice regular pattern, but it is not nice to look at when broken up. The stippled example is much nicer:

n <- 100
x <- 1:n
y <- 1:n
M <- list(x=x, y=y, z=outer(x, y, FUN = function(x,y){x^2 * y * rlnorm(n^2,0,0.2)}))
image(M)
range(M$z)
incl <- which(M$z>5e5)

polys <- matrix.poly(M$x, M$y, z=M$z, n=incl)

png("hatched_image.png", height=5, width=5, units="in", res=400)
op <- par(mar=c(3,3,1,1))
image(M)
for(i in seq(polys)){
  polygon(polys[[i]], density=10, angle=45, border=NA, lwd=0.5)
  polygon(polys[[i]], density=10, angle=-45, border=NA, lwd=0.5)
}
box()
par(op)
dev.off()

png("stippled_image.png", height=5, width=5, units="in", res=400)
op <- par(mar=c(3,3,1,1))
image(M)
grd <- expand.grid(x=x, y=y)
points(grd$x[incl], grd$y[incl], pch=".", cex=1.5)
box()
par(op)
dev.off()

Answer 2

Do it using the coordinate positioning mechanism of ?image [1].

data(volcano)
m <- volcano
dimx <- nrow(m)
dimy <- ncol(m)

d1 <- list(x = seq(0, 1, length = dimx), y = seq(0, 1, length = dimy), z = m)

With your 'image' constructed that way you keep the structure with the object, and its coordinates intact. You can collect multiple matrices into a 3D array or as multiple elements, but you need to augment image() in order to handle that, so I keep them separate here.

Make a copy of the data to specify an interesting area.

d2 <- d1
d2$z <- d2$z > 155

Use the coordinates to specify which cells are interesting. This is expensive if you have a very big raster, but it's super easy to do.

pts <- expand.grid(x = d2$x, y = d2$y)
pts$over <- as.vector(d2$z)

Set up the plot.

op <- par(mfcol = c(2, 1))
image(d1)

image(d1)
points(pts$x[pts$over], pts$y[pts$over], cex = 0.7)

par(op)

Don't forget to modify the plotting of points to get different effects, in particular a very dense grid with lots of points will take ages to draw all those little circles. pch = "." is a good choice.

Now, do you have some real data to plot on that nice projection? See examples here for some of the options: http://spatial-analyst.net/wiki/index.php?title=Global_datasets

[1] R has classes for more sophisticated handling of raster data, see package sp and raster for two different approaches.