Alterando aleatoriamente o mapa raster dos tipos de habitat?

Eu tenho uma variedade de tipos de habitat para uma área específica na Escócia. Preciso criar futuros cenários de habitat com mudanças no habitat para avaliar a viabilidade populacional de uma espécie de ave.

Por exemplo, no futuro, pode haver 10% mais florestas na área. Gostaria de alterar o mapa atual adicionando aleatoriamente a silvicultura em blocos de um determinado tamanho. Até agora, estou pensando nas linhas de seleção de pontos aleatórios de uma varredura que identifica áreas onde a silvicultura pode ocorrer e o crescimento dos blocos de tamanho correto usando algum tipo de autômato celular.

Essa parece a melhor maneira de fazer isso? Há um método melhor?

Se esta é a melhor maneira disponível, como eu poderia fazer isso, preferencialmente, R? (Atualmente, estou olhando a função rpoints em "spatstat" junto com o pacote CellularAutomata)

Também tenho acesso ao GRASS, QGis e ArcMap 10, se houver maneiras mais simples em qualquer uma delas.

Você já pensou em usar uma cadeia de Markov ? Este é efetivamente um "autômato celular probabilístico", fornecendo assim a aleatoriedade desejada. Em vez de prescrever a nova geração em termos de vizinhos locais da geração existente, especifica uma distribuição de probabilidade para a nova geração. Essa distribuição pode ser estimada a partir de, digamos, seqüências de tempo de imagens da mesma área ou de áreas semelhantes.

Intuitivamente, este modelo diz que uma célula não fará necessariamente uma transição de floresta para não florestal (ou vice-versa ), mas as chances de que ela faça a transição dependem da cobertura da terra imediatamente ao seu redor. Ele pode lidar com várias classes de cobertura, configurações complexas de bairros e até ser generalizado para "lembrar" a história recente da evolução da cobertura do solo.

As transições podem ser implementadas usando instruções de Álgebra de Mapa, que tornam esse método viável em qualquer GIS baseado em varredura, mesmo naquelas sem acesso direto ou rápido aos dados no nível da célula. Usar R torna ainda mais fácil.

Por exemplo, considere esta configuração inicial com apenas duas classes, branca e preta:

Para ilustrar o que pode acontecer, criei um modelo parametrizado (não com base em dados) no qual a transição para o preto ocorre com probabilidade 1 - q ^ k, em que k é o número médio de células negras na vizinhança 3 por 3 (k = 0, 1/9, 2/9, ..., 1). Quando q é pequeno ou a maior parte da vizinhança já está preta, a nova célula ficará preta. Aqui estão quatro simulações independentes da décima geração para cinco valores de q variando de 0,25 a 0,05:

Evidentemente, esse modelo possui muitas das características de uma AC, mas também inclui um efeito aleatório útil para explorar resultados alternativos.

Código

O seguinte implementa a simulação em 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])))