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

This is a solution in the spirit of @mdsummer's comment using ggplot2. I first draw the grid, and then draw +'es at the locations where a certain value has been exceeded. Note that ggplot2 works with data.frame's, not with multi-dimensional arrays or matrices. You can use melt from the reshape package to convert from an array / marix to a data.frame flat structure.

Here is a concrete example using the example data from the geom_tile documentation:

pp <- function (n,r=4) { 
 x <- seq(-r*pi, r*pi, len=n) 
 df <- expand.grid(x=x, y=x) 
 df$r <- sqrt(df$x^2 + df$y^2) 
 df$z <- cos(df$r^2)*exp(-df$r/6) 
 df 
} 

require(ggplot2)
dat = pp(200)
over = dat[,c("x","y")]
over$value = with(dat, ifelse(z > 0.5, 1, 0))
ggplot(aes(x = x, y = y), data = dat) + 
   geom_raster(aes(fill = z)) + 
   scale_fill_gradient2() +
   geom_point(data = subset(over, value == 1), shape = "+", size = 1)

Answer 3

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.

Answer 4

This is probably coming too late, but I'd like to post my answer as a reference too.

One nice option for spatial data is to use the rasterVis package. Once you have a "base" raster object, and the "mask" object, which you will use to draw the stippling, you can do something like:

require(raster)
require(rasterVis)

# Scratch raster objects
data(volcano)
r1 <- raster(volcano)

# Here we are selecting only values from 160 to 180.
# This will be our "mask" layer.
over <- ifelse(volcano >=160 & volcano <=180, 1, NA) 
r2 <- raster(over)

# And this is the key step:
# Converting the "mask" raster to spatial points
r.mask <- rasterToPoints(r2, spatial=TRUE)

# Plot
levelplot(r1, margin=F) +
layer(sp.points(r.mask, pch=20, cex=0.3, alpha=0.8))

which resembles the map that the OP was looking for. Parameters of the points such as color, size and type can be fine tuned. ?sp.points provides all the arguments that can be used to do that.