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.

  • Have you looked at the raster package yet? It's got a lot of tools to work with raster (noo, rly?) data. Oct 24, 2011 at 17:27
  • Thanks, Roman. Yes, this should give me the tools to read and manipulate my base map.
    – Matt Geary
    Oct 25, 2011 at 10:38

1 Answer 1


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:

Land cover grid

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:

Table of results

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


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)

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
  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.
       function(i) image(y.list[[i]], 
                         col=c("White", "Black"),
  • +1 Very interesting. If you had historic landcover data for a particular area, would it be possible to derive q and/or k? Oct 24, 2011 at 19:15
  • 2
    @Kirk Yes, but I wouldn't recommend it: the model I used for illustration is too simplistic. But you can derive something better: by looking at the empirical frequencies of transitions out of each neighborhood configuration that has occurred, you can create models of future evolution whose transitions statistically emulate the past evolution. If transition frequencies are spatially homogeneous and if the future continues to act like the past, running a few of these simulations can give a clear picture of what the future is likely to hold.
    – whuber
    Oct 24, 2011 at 21:06
  • Thanks, this seems to do exactly what I need. Would it be possible to set a limit on the proportion of the area which changes?
    – Matt Geary
    Oct 25, 2011 at 10:37
  • @Matt Yes, at least in a probabilistic sense. The theory describes how Markov chains reach an asymptotically stable mixture of proportions of each state. This is a dynamic equilibrium: at each generation a lot of cells may be changing, but the net result is to keep their proportions within the grid the same (up to small chance deviations).
    – whuber
    Oct 25, 2011 at 13:56
  • 1
    I'm a terrible R programmer. I can share the Mathematica code I used; with R's apply functions, it ought to port well. You need a kernel, a transition rule, and a procedure to apply them to a 2D 0/1 array. Thus: kernel = ConstantArray[1/3^2, {3,3}] for the kernel; transitionRule [k_] := With[{q = 0.1}, Boole[RandomReal[{0, 1}] > q^k]] for the rule; and next[a_, kernel_, f_] := Map[f, ListConvolve[kernel, a, {1, 1}, 0], {2}] to apply them to an array a. E.g., to plot four generations from start, use ArrayPlot /@ NestList[next[#, kernel, transitionRule] &, start, 3].
    – whuber
    Nov 15, 2011 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.