Creación de grupos de células de forma aleatoria en una trama a partir de semillas de 1 célula/píxel

Question

Deseo "cultivar" grupos de células a partir de semillas dentro de una trama. Mi trama base está llena de 1's y 0's, los 1's indican tierra y los 0's mar/áreas NA. De los 1's deseo seleccionar 60 pixels/células al azar como mis semillas, y entonces hacer crecer al azar un grupo conectado de un límite predefinido de pixels/células de esa semilla. He oído que la técnica puede ser referido como 'tinte propagación', pero no he tenido suerte de encontrar mucho en él. La célula semilla se convertiría en un valor de 2 y luego la siguiente célula elegida de los 1 circundantes se convertiría en 2 también. Los 2 no están disponibles para ser convertidos en el futuro.

¿Alteración aleatoria del mapa raster de tipos de hábitat? ayuda un poco, ya que también estaría dispuesto a hacer esto en R, ya que estoy familiarizado con la lectura y manipulación de datos SIG en R. Sin embargo, lo que necesito es un conjunto de reglas para seleccionar al azar los píxeles que rodean un grupo existente.

Tal vez alguien haya realizado esta forma más básica de autómatas celulares en un entorno SIG.

Ejemplo:

Tengo un objetivo de 250 células. Selecciono aleatoriamente una celda que tiene un valor de 1. Esta se convierte en un valor de 2. A continuación, uno de los vecinos de la celda semilla que = 1 se convierte en un 2. A continuación, uno de los vecinos de cualquiera de las celdas con un valor de 2 se selecciona y se convierte en un 2. Esto continuaría hasta que se alcance una forma continua que numere 250 celdas.

Basado en la gran respuesta de @whuber, tengo algunas preguntas sobre el código:

1. ¿Cómo consigo asignar los valores de las celdas crecidas a sólo un '2'? en lugar de valores variables que representan el orden en que fueron creadas?
2. Necesito crear 60 grupos de celdas dentro de mi área de "1". He ideado formas de elegir posiciones de inicio aleatorias, pero ¿cómo puedo hacer que todo funcione dentro de un bucle utilizando la función expand función sobre la que escribió?

¿Puede sugerirnos una forma de crear 60 grupos que no choquen entre sí y estén contenidos en la misma matriz final?

Explicación del problema

Cada grupo de celdas representará un área protegida de un tamaño definido, por ejemplo, 250 celdas. Cada área tiene que empezar y crecer en celdas con un valor de 1, ya que esto representa la tierra, y evitar las celdas con un valor de 0, ya que esto representa el mar. Necesito iterar esto 1000 veces con 60 áreas protegidas en cada iteración para crear un modelo nulo, mostrando cuales serán las distribuciones de estas áreas por azar. Por esta razón, el recuento total de celdas en las 60 áreas tendrá que ser idéntico en cada una de las 1000 iteraciones para que sean comparables. Por lo tanto, no pasa nada si las áreas se tocan, pero si se produce una colisión, lo ideal sería que el grupo creciera en otra dirección disponible hasta alcanzar el objetivo de 250. Futuras modificaciones a este proceso implicarán variar el tamaño de cada una de las 60 áreas y potencialmente la forma aunque entiendo que esto sería difícil.

Una vez creadas cada una de estas 1.000 redes de áreas protegidas, se utilizarán como máscara frente a otros datos rasterizados, como las medidas de biodiversidad, para ver (a) si se cruzan con áreas de distribución de especies concretas y (b) qué % de áreas de distribución de especies concretas cubren estas redes aleatorias de áreas protegidas.

Answer 1

Ofrezco un R solución que se codifica en un ligeramente no-R manera de ilustrar cómo podría ser abordado en otras plataformas.

La preocupación en R (así como de otras plataformas, especialmente aquellos que están a favor de un estilo de programación funcional) es que la actualización constante de una gran cantidad puede ser muy caro. En su lugar, a continuación, este algoritmo mantiene su propia estructura de datos privados en el que (a) todas las células que han sido llenados hasta el momento están en la lista y (b) todas las células que están disponibles para ser elegido (alrededor del perímetro de las células llenas) están en la lista. Aunque la manipulación de esta estructura de datos es menos eficiente que directamente la indización en una matriz, manteniendo los datos modificados en un tamaño pequeño, es probable que tome mucho menos tiempo de cálculo. (No se ha hecho un esfuerzo para optimizarlo para R,. Pre-asignación de los vectores de estado debe ahorrar algo de tiempo de ejecución, si usted prefiere seguir trabajando en R.)

El código está comentado y deben ser fáciles de leer. Para hacer que el algoritmo sea lo más completo posible, que no hace uso de add-ons, excepto en la final a la trama el resultado. La única parte difícil es que para la eficiencia y la simplicidad prefiere índice en las mallas 2D mediante el uso de 1D índices. Una conversión que sucede en la neighbors función de las necesidades de la 2D de indexación con el fin de averiguar lo que el acceso de los vecinos de una célula puede ser y luego los convierte en la 1D índice. Esta conversión es estándar, así que no voy a comentar más, excepto para señalar que en otras plataformas GIS puede que desee invertir los roles de la columna y la fila de los índices. (En R, la fila de los índices de cambio antes de los índices de las columnas.)

Para ilustrar, este código tiene una cuadrícula x que representa a la tierra y a un río-como característica de puntos inaccesibles, se inicia en un lugar específico (5, 21) en esa red (cerca de la parte inferior curva del río) y se expande al azar para cubrir los 250 puntos. Total de temporización es de 0,03 segundos. (Cuando el tamaño de la matriz es mayor por un factor de 10.000 a 3000 filas por 5000 columnas, el momento en que sube sólo a 0.09 segundos-un factor de sólo 3--demostración de la escalabilidad de este algoritmo.) En lugar de simplemente generar una cuadrícula de 0, 1, y 2, salidas de la secuencia con la que las nuevas células fueron asignados. En la figura las primeras células son de color verde, que se graduó en la medallas de oro en el salmón colores.

Debe ser evidente que un niño de ocho punto vecindario de cada celda está siendo utilizado. Para otros barrios, simplemente modificar el nbrhood valor cerca del principio de la expand: es una lista de índice de compensaciones en relación a cualquier célula dada. Por ejemplo, una "D4" barrio podría ser especificado como matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

También es evidente que este método de propagación tiene sus problemas: se deja orificios detrás. Si eso no es lo que era la intención, hay varias formas de solucionar este problema. Por ejemplo, mantener la disposición de las células en una cola, de modo que las primeras células que se encuentran son también los más antiguos se llenó. Algunos aleatorización puede aplicarse, pero la disposición de las células dejará de ser elegido con uniforme (igual) de probabilidades. Otra, más complicada, sería para seleccionar celdas disponibles con probabilidades que dependen de la cantidad de lleno a los vecinos que tienen. Una vez que una célula se convierte rodeado, usted puede hacer su probabilidad de selección tan alto que algunos de los agujeros quedaría sin cubrir.

Voy a terminar comentando que este no es un autómata celular (CA), que no procedería de la célula por célula, sino de la actualización de bloques enteros de células en cada generación. La diferencia es sutil: con la CA, las probabilidades de selección de las células no sería uniforme.

#
# Expand a patch randomly within indicator array `x` (1=unoccupied) by
# `n.size` cells beginning at index `start`.
#
expand <- function(x, n.size, start) {
  if (x[start] != 1) stop("Attempting to begin on an unoccupied cell")
  n.rows <- dim(x)[1]
  n.cols <- dim(x)[2]
  nbrhood <- matrix(c(-1,-1, -1,0, -1,1, 0,-1, 0,1, 1,-1, 1,0, 1,1), nrow=2)
  #
  # Adjoin one more random cell and update `state`, which records
  # (1) the immediately available cells and (2) already occupied cells.
  #
  grow <- function(state) {
    #
    # Find all available neighbors that lie within the extent of `x` and
    # are unoccupied.
    #
    neighbors <- function(i) {
      n <- c((i-1)%%n.rows+1, floor((i-1)/n.rows+1)) + nbrhood
      n <- n[, n[1,] >= 1 & n[2,] >= 1 & n[1,] <= n.rows & n[2,] <= n.cols, 
             drop=FALSE]             # Remain inside the extent of `x`.
      n <- n[1,] + (n[2,]-1)*n.rows  # Convert to *vector* indexes into `x`.
      n <- n[x[n]==1]                # Stick to valid cells in `x`.
      n <- setdiff(n, state$occupied)# Remove any occupied cells.
      return (n)
    }
    #
    # Select one available cell uniformly at random.
    # Return an updated state.
    #
    j <- ceiling(runif(1) * length(state$available))
    i <- state$available[j]
    return(list(index=i,
                available = union(state$available[-j], neighbors(i)),
                occupied = c(state$occupied, i)))
  }
  #
  # Initialize the state.
  # (If `start` is missing, choose a value at random.)
  #
  if(missing(start)) {
    indexes <- 1:(n.rows * n.cols)
    indexes <- indexes[x[indexes]==1]
    start <- sample(indexes, 1)
  }
  if(length(start)==2) start <- start[1] + (start[2]-1)*n.rows
  state <- list(available=start, occupied=c())
  #
  # Grow for as long as possible and as long as needed.
  #
  i <- 1
  indices <- c(NA, n.size)
  while(length(state$available) > 0 && i <= n.size) {
    state <- grow(state)
    indices[i] <- state$index
    i <- i+1
  }
  #
  # Return a grid of generation numbers from 1, 2, ... through n.size.
  #
  indices <- indices[!is.na(indices)]
  y <- matrix(NA, n.rows, n.cols)
  y[indices] <- 1:length(indices)
  return(y)
}
#
# Create an interesting grid `x`.
#
n.rows <- 3000
n.cols <- 5000
x <- matrix(1, n.rows, n.cols)
ij <- sapply(1:n.cols, function(i) 
      c(ceiling(n.rows * 0.5 * (1 + exp(-0.5*i/n.cols) * sin(8*i/n.cols))), i))
x[t(ij)] <- 0; x[t(ij - c(1,0))] <- 0; x[t(ij + c(1,0))] <- 0
#
# Expand around a specified location in a random but reproducible way.
#
set.seed(17)
system.time(y <- expand(x, 250, matrix(c(5, 21), 1)))
#
# Plot `y` over `x`.
#
library(raster)
plot(raster(x[n.rows:1,], xmx=n.cols, ymx=n.rows), col=c("#2020a0", "#f0f0f0"))
plot(raster(y[n.rows:1,] , xmx=n.cols, ymx=n.rows), 
     col=terrain.colors(255), alpha=.8, add=TRUE)

---

Con ligeras modificaciones podemos recorrer expand crear varios grupos. Es conveniente diferenciar los grupos de un identificador, que aquí se va a ejecutar 2, 3, ..., etc.

En primer lugar, el cambio expand a devolver (a) NA en la primera línea si no es un error, y (b) los valores en indices en lugar de una matriz y. (No pierdas el tiempo en la creación de una nueva matriz y con cada llamada.) Con este cambio, el bucle es fácil: elija un arranque aleatorio, intenta ampliar su alrededor, se acumulan los índices clúster en indices en caso de éxito, y repetir hasta que esté hecho. Una parte clave del bucle es limitar el número de iteraciones en el caso de que muchos de los clústeres contiguos no se encuentra: esto se hace con count.max.

Aquí está un ejemplo en 60 centros de cluster son elegidos de manera uniforme al azar.

size.clusters <- 250
n.clusters <- 60
count.max <- 200
set.seed(17)
system.time({
  n <- n.rows * n.cols
  cells.left <- 1:n
  cells.left[x!=1] <- -1 # Indicates occupancy of cells
  i <- 0
  indices <- c()
  ids <- c()
  while(i < n.clusters && length(cells.left) >= size.clusters && count.max > 0) {
    count.max <- count.max-1
    xy <- sample(cells.left[cells.left > 0], 1)
    cluster <- expand(x, size.clusters, xy)
    if (!is.na(cluster[1]) && length(cluster)==size.clusters) {
      i <- i+1
      ids <- c(ids, rep(i, size.clusters))
      indices <- c(indices, cluster)
      cells.left[indices] <- -1
    }                
  }
  y <- matrix(NA, n.rows, n.cols)
  y[indices] <- ids
})
cat(paste(i, "cluster(s) created.", sep=" "))

Aquí está el resultado cuando se aplica a un 310 por 500 rejilla (de hecho lo suficientemente pequeña y gruesa para los grupos a ser evidente). Se tarda dos segundos en ejecutarse; en un 3100 de 5000 cuadrícula (100 veces más grande) lleva más tiempo (24 segundos) pero el momento es de escala razonablemente bien. (En otras plataformas, tales como C++, el momento en que debiera depender del tamaño de la cuadrícula.)

Answer 2

Sin haber realizado la operación, y no tenía tiempo de sobra para jugar, sólo puedo añadir estos dos enlaces a la lista:

Encontrar el más Cercano de Celda Ráster Valor Basado en un Vector de Punto (La primera respuesta (con 4 votos) es lo que me ha intrigado).

También: ¿ Hawth del Gridspread ayuda?