r - Randomly altering raster map of habitat types?

I have a raster of habitat types for a specific area in Scotland. I need to create future habitat scenarios with changes in habitat to assess the population viability of a bird species.

For example, in the future there might be 10% more forestry in the area. I would like to alter the current map by randomly adding the forestry in blocks of a certain size. I am, so far, thinking along the lines of selecting random points from a raster which identifies areas where the forestry could occur and growing the correct sized blocks using some sort of cellular automata.

Does this seem like the best way of going about this? Is there a better method?

If this is the best way available, how could I do this in, preferably, R? (I am currently looking at the rpoints function in "spatstat" along with the CellularAutomata package)

I also have access to GRASS, QGis and ArcMap 10 if there are simpler ways in any of them.

Answer

Have you thought of using a Markov chain? This is effectively a "probabilistic cellular automaton," thereby supplying the desired randomness. Instead of prescribing the new generation in terms of local neighbors of the existing generation, it specifies a probability distribution for the new generation. That distribution can be estimated from, say, time sequences of images of the same or similar areas.

Intuitively, this model says that a cell will not necessarily make a transition from forested to non-forested (or vice versa), but the chances that it will make the transition depend on the land cover immediately around it. It can handle multiple classes of cover, complex configurations of neighborhoods, and even be generalized to "remember" the recent history of land cover evolution.

The transitions can be implemented using Map Algebra statements, which makes this method practicable in any raster-based GIS, even those without direct or quick access to cell-level data. Using R makes it even easier.

For example, consider this starting configuration with just two classes, white and black:

To illustrate what can happen, I created a parameterized model (not based on any data) in which the transition to black occurs with probability 1 - q^k where k is the average number of black cells within the 3 by 3 neighborhood (k = 0, 1/9, 2/9, ..., 1). When either q is small or most of the neighborhood is already black, the new cell will be black. Here are four independent simulations of the tenth generation for five values of q ranging from 0.25 down to 0.05:

Evidently this model has many of the characteristics of a CA but it also includes a random effect useful for exploring alternative outcomes.

---

Code

The following implements the simulation in R.

#
# Make a transition from state `x` using a kernel having `k.ft` as
# its Fourier transform.
#

transition <- function(x, k.ft, q=0.1) {
  k <- zapsmall(Re(fft(k.ft * fft(x), inverse=TRUE))) / length(x)
  matrix(runif(length(k)) > q^k, nrow=nrow(k))
}
#
# Create the zeroth generation and the fft of a transition kernel.
#
n.row <- 2^7 # FFT is best with powers of 2
n.col <- 2^7
kernel <- matrix(0, nrow=n.row, ncol=n.col)

kernel[1:3, 1:3] <- 1/9
kernel.f <- fft(kernel)

set.seed(17)
x <- matrix(sample(c(0,1), n.row*n.col, replace=TRUE, prob=c(599, 1)), n.row)
#
# Prepare to run multiple simulations.
#
y.list <- list()
parameters <- c(.25, .2, .15, .1, .05)

#
# Perform and benchmark the simulations.
#
i <- 0
system.time({
  for (q in parameters) {
    y <- x
    for (generation in 1:10) {
      y <- transition(y, kernel.f, q)
    }

    y.list[[i <- i+1]] <- y
  }
})
#
# Display the results.
#    
par(mfrow=c(1,length(parameters)))
invisible(sapply(1:length(parameters), 
       function(i) image(y.list[[i]], 
                         col=c("White", "Black"),

                         main=parameters[i])))