¿Crear grupos de celdas de forma aleatoria en una trama a partir de semillas de 1 celda / píxel?

Como dice mi título, deseo "hacer crecer" grupos de células a partir de semillas dentro de una trama. Mi ráster base está lleno de 1 y 0, 1 indica áreas terrestres y 0 mar / NA. De los 1 deseo seleccionar 60 píxeles / celdas aleatorias como mis semillas, y luego crecer aleatoriamente un grupo conectado de un no predefinido. de píxeles / celdas límite de esa semilla. He oído que la técnica se puede denominar 'tinte extendido', pero no he tenido suerte de encontrarla. La celda de semilla se convertiría en un valor de 2 y luego la siguiente celda elegida de los alrededores de 1 se convertiría también en 2. Los 2 no están disponibles para convertirse en el futuro.

Este hilo ayuda un poco, ya que también estaría dispuesto a hacerlo 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 píxeles que rodean un grupo existente .

Si alguien ha hecho esta forma más básica de autómatas celulares en un entorno SIG, agradecería cualquier consejo / orientación.

Ejemplo:

Tengo un objetivo de 250 celdas. Selecciono aleatoriamente una celda que tiene un valor de 1. Esto se convierte en un valor de 2. Luego, uno de los vecinos de la celda semilla que = 1 se convierte en un 2. Luego, uno de los vecinos de cualquiera de las celdas con un valor 2 seleccionado y convertido en 2. Esto continuaría hasta que se alcance una forma continua que numere 250 celdas.

Editar: más preguntas

Según la gran respuesta de whuber, tengo algunas preguntas sobre el código:

1. ¿Cómo puedo asignar los valores de las celdas crecidas a solo un '2' en lugar de valores variables que representan el orden en que se crearon?
2. Necesito crear 60 grupos de celdas dentro de mi área de '1'. He ideado formas de elegir posiciones iniciales al azar, pero me cuesta hacer que todo funcione en un bucle utilizando la expandfunción que amablemente escribió. ¿Puede sugerir una forma de crear 60 grupos que no choquen entre sí y que estén contenidos en la misma matriz final?

Editar: Explicación adicional del problema

Cada grupo de celdas representará un área protegida de un tamaño definido, por ejemplo, 250 celdas. Cada área necesita comenzar y crecer en celdas con un valor de 1, ya que esto representa la tierra y evitar 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, que muestre qué distribuciones de estas áreas serán por casualidad. Por esta razón, el recuento total de celdas en las 60 áreas deberá ser idéntico en cada una de las 1000 iteraciones para que sean comparables. Por lo tanto, está bien si las áreas se tocan, pero si hay una colisión, lo ideal sería que el grupo creciera en otra dirección disponible hasta alcanzar el objetivo de 250.

Una vez que se crean cada una de estas 1000 redes de áreas protegidas, se utilizarán como una máscara contra otros datos ráster, como las medidas de biodiversidad para ver (a) si intersecan rangos de especies particulares y (b) qué% de especies particulares rangos de estas redes aleatorias de áreas protegidas cubiertas.

Gracias a @whuber por tu ayuda hasta ahora, no espero que pases más tiempo ayudándome, pero pensé que intentaría aclarar mi situación como me lo solicitaste.

Ofreceré una Rsolución que está codificada de una manera ligeramente diferente Rpara ilustrar cómo podría abordarse en otras plataformas.

La preocupación en R(así como en algunas otras plataformas, especialmente aquellas que favorecen un estilo de programación funcional) es que actualizar constantemente una gran matriz puede ser muy costoso. En cambio, este algoritmo mantiene su propia estructura de datos privados en la que (a) se enumeran todas las celdas que se han llenado hasta ahora y (b) todas las celdas que están disponibles para ser elegidas (alrededor del perímetro de las celdas llenas) están listados. Si bien la manipulación de esta estructura de datos es menos eficiente que la indexación directa en una matriz, al mantener los datos modificados en un tamaño pequeño, es probable que tome mucho menos tiempo de cálculo. (Tampoco se ha hecho ningún esfuerzo para optimizarlo R. La preasignación de los vectores de estado debería ahorrar algo de tiempo de ejecución, si prefiere seguir trabajando dentro R).

El código está comentado y debe ser fácil de leer. Para que el algoritmo sea lo más completo posible, no hace uso de ningún complemento, excepto al final para trazar el resultado. La única parte difícil es que, por eficiencia y simplicidad, prefiere indexar en las cuadrículas 2D mediante el uso de índices 1D. Se produce una conversión en la neighborsfunción, que necesita la indexación 2D para determinar cuáles podrían ser los vecinos accesibles de una celda y luego los convierte en el índice 1D. Esta conversión es estándar, por lo que no la comentaré más, excepto para señalar que en otras plataformas SIG es posible que desee revertir los roles de los índices de columna y fila. (En R, los índices de fila cambian antes que los índices de columna).

Para ilustrar, este código toma una cuadrícula que xrepresenta la tierra y una característica similar a un río de puntos inaccesibles, comienza en una ubicación específica (5, 21) en esa cuadrícula (cerca de la curva inferior del río) y la expande al azar para cubrir 250 puntos . El tiempo total es de 0.03 segundos. (Cuando el tamaño de la matriz se incrementa en un factor de 10,000 a 3000 filas en 5000 columnas, el tiempo aumenta solo a 0.09 segundos, un factor de solo 3 más o menos, lo que demuestra la escalabilidad de este algoritmo). simplemente generando una cuadrícula de 0, 1 y 2, genera la secuencia con la que se asignaron las nuevas celdas. En la figura, las primeras celdas son verdes, graduándose a través de dorados en colores salmón.

Debería ser evidente que se está utilizando una vecindad de ocho puntos de cada celda. Para otros vecindarios, simplemente modifique el nbrhoodvalor cerca del comienzo de expand: es una lista de compensaciones de índice en relación con cualquier celda dada. Por ejemplo, un vecindario "D4" podría especificarse 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: deja agujeros atrás. Si eso no es lo que se pretendía, hay varias formas de solucionar este problema. Por ejemplo, mantenga las celdas disponibles en una cola para que las primeras celdas encontradas también se llenen de las primeras. Todavía se puede aplicar algo de aleatorización, pero las celdas disponibles ya no se elegirán con probabilidades uniformes (iguales). Otra forma más complicada sería seleccionar celdas disponibles con probabilidades que dependan de cuántos vecinos llenos tengan. Una vez que una celda se rodea, puede hacer que su posibilidad de selección sea tan alta que quedarían pocos agujeros sin llenar.

Terminaré comentando que este no es un autómata celular (CA), que no procedería celda por celda, sino que actualizaría franjas enteras de células en cada generación. La diferencia es sutil: con la CA, las probabilidades de selección para las celdas no serían uniformes.

#
# 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 hacer un bucle expandpara crear múltiples grupos. Es aconsejable diferenciar los grupos por un identificador, que aquí ejecutará 2, 3, ..., etc.

Primero, cambie expandpara devolver (a) NAen la primera línea si hay un error y (b) los valores en indiceslugar de una matriz y. (No pierda el tiempo creando una nueva matriz ycon cada llamada). Con este cambio realizado, el bucle es fácil: elija un inicio aleatorio, intente expandirse, acumule los índices de clúster indicessi tiene éxito y repita hasta que termine. Una parte clave del ciclo es limitar el número de iteraciones en caso de que no se puedan encontrar muchos grupos contiguos: esto se hace con count.max.

Aquí hay un ejemplo en el que se seleccionan 60 centros agrupados 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 una cuadrícula de 310 por 500 (hecha lo suficientemente pequeña y gruesa para que los grupos sean evidentes). Tarda dos segundos en ejecutarse; en una cuadrícula de 3100 por 5000 (100 veces más grande) lleva más tiempo (24 segundos), pero el tiempo está escalando razonablemente bien. (En otras plataformas, como C ++, el tiempo no debería depender del tamaño de la cuadrícula).

Sin haber realizado su operación y sin tiempo de sobra para jugar, solo puedo agregar estos dos enlaces a su lista:

Encuentre el valor de celda ráster más cercano basado en un punto vectorial (La primera respuesta (con 4 votos) es lo que me intrigó).

Además: ¿ayudaría Hawth's Gridspread ?