¿Alterar aleatoriamente el mapa ráster de los tipos de hábitat?

Tengo una trama de tipos de hábitat para un área específica en Escocia. Necesito crear escenarios de hábitat futuros con cambios en el hábitat para evaluar la viabilidad de la población de una especie de ave.

Por ejemplo, en el futuro podría haber un 10% más de silvicultura en el área. Me gustaría alterar el mapa actual agregando aleatoriamente la silvicultura en bloques de cierto tamaño. Hasta ahora, estoy pensando en seleccionar líneas aleatorias de un ráster que identifica áreas donde podría ocurrir la silvicultura y hacer crecer los bloques del tamaño correcto utilizando algún tipo de autómata celular.

¿Parece esta la mejor manera de hacerlo? hay algun metodo mejor?

Si esta es la mejor manera disponible, ¿cómo podría hacer esto, preferiblemente, en R? (Actualmente estoy viendo la función rpoints en "spatstat" junto con el paquete CellularAutomata)

También tengo acceso a GRASS, QGis y ArcMap 10 si hay formas más simples en cualquiera de ellos.

¿Has pensado en usar una cadena de Markov ? Esto es efectivamente un "autómata celular probabilístico", proporcionando así la aleatoriedad deseada. En lugar de prescribir la nueva generación en términos de vecinos locales de la generación existente, especifica una distribución de probabilidad para la nueva generación. Esa distribución puede estimarse a partir de, por ejemplo, secuencias temporales de imágenes de las mismas áreas o áreas similares.

Intuitivamente, este modelo dice que una célula no necesariamente hará una transición de bosque a bosque (o viceversa ), pero las posibilidades de que haga la transición dependerán de la cobertura del suelo inmediatamente a su alrededor. Puede manejar múltiples clases de cobertura, configuraciones complejas de vecindarios, e incluso generalizarse para "recordar" la historia reciente de la evolución de la cobertura del suelo.

Las transiciones se pueden implementar utilizando declaraciones de Álgebra de mapas, lo que hace que este método sea factible en cualquier SIG basado en ráster, incluso aquellos sin acceso directo o rápido a datos a nivel de celda. Usar R lo hace aún más fácil.

Por ejemplo, considere esta configuración inicial con solo dos clases, blanco y negro:

Para ilustrar lo que puede suceder, creé un modelo parametrizado (no basado en ningún dato) en el que la transición al negro ocurre con probabilidad 1 - q ^ k donde k es el número promedio de celdas negras dentro del vecindario 3 por 3 (k = 0, 1/9, 2/9, ..., 1). Cuando q es pequeño o la mayor parte del vecindario ya es negro, la nueva celda será negra. Aquí hay cuatro simulaciones independientes de la décima generación para cinco valores de q que van desde 0.25 hasta 0.05:

Evidentemente, este modelo tiene muchas de las características de una AC pero también incluye un efecto aleatorio útil para explorar resultados alternativos.

Código

Lo siguiente implementa la simulación en 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])))