Criando grupos de células de formato aleatório em uma varredura a partir de sementes de 1 célula / pixel?

Como meu título diz, desejo "cultivar" grupos de células a partir de sementes dentro de uma varredura. Minha base raster está cheia de zeros e zeros, zeros indicam áreas terrestres e zeros para o mar / NA. A partir dos anos 1, eu desejo selecionar 60 pixels / células aleatórios como minhas sementes e, em seguida, cultivar aleatoriamente um grupo de um não predefinido. limite de pixels / células a partir dessa semente. Ouvi dizer que a técnica pode ser chamada de 'espalhar corante', mas não tive muita sorte em encontrar muita coisa sobre ela. A célula-semente seria transformada em um valor de 2 e, em seguida, a próxima célula escolhida entre os 1s circundantes também seria convertida em 2. 2 estão indisponíveis para serem convertidos no futuro.

Esse segmento ajuda um pouco, pois eu também estaria disposto a fazer isso em R, pois estou familiarizado com a leitura e manipulação de dados GIS em R. No entanto, o que eu preciso é de um conjunto de regras para selecionar aleatoriamente pixels ao redor de um grupo existente .

Se alguém tiver feito essa forma mais básica de autômato celular em um ambiente GIS, eu agradeceria qualquer conselho / orientação.

Exemplo:

Eu tenho um alvo de 250 células. Seleciono aleatoriamente uma célula que tem um valor de 1. Isso é transformado em um valor de 2. Em seguida, um dos vizinhos da célula-semente que = 1 é transformado em um 2. Em seguida, um dos vizinhos das células com um valor 2 é selecionado e girado para um 2. Isso continuará até que uma forma contínua numerada em 250 células seja atingida.

Editar: Outras perguntas

Com base na ótima resposta do whuber, tenho algumas perguntas sobre o código:

1. Como faço para alocar os valores das células cultivadas para apenas um '2', em vez de valores variáveis ​​que representam a ordem em que foram criadas?
2. Preciso criar 60 aglomerados de células na minha área de '1's. Eu inventei maneiras de escolher posições iniciais aleatórias, mas luto para fazer tudo funcionar em um loop usando a expandfunção que você gentilmente escreveu. Você pode sugerir uma maneira de criar 60 grupos que não se chocam entre si e estão contidos na mesma matriz final?

Editar: Explicação adicional do problema

Cada grupo de células representará uma área protegida de tamanho definido, por exemplo, 250 células. Cada área precisa iniciar e crescer em células com o valor 1, pois isso representa terra e evitar células com o valor 0, pois isso representa o mar. Preciso iterar isso 1000 vezes com 60 áreas protegidas em cada iteração para criar um modelo nulo, mostrando quais distribuições dessas áreas serão por acaso. Por esse motivo, a contagem total de células em todas as 60 áreas precisará ser idêntica em cada uma das 1000 iterações para que sejam comparáveis. Portanto, tudo bem se as áreas tocarem, mas se houver uma colisão, o ideal é que o grupo cresça em outra direção disponível até que a meta de 250 seja atingida.

Depois que cada uma dessas 1000 redes de áreas protegidas for criada, elas serão usadas como uma máscara contra outros dados rasterizados, como medidas de biodiversidade, para verificar (a) se eles cruzam faixas de espécies específicas e (b) qual% de espécies específicas varia essas redes aleatórias de áreas protegidas.

Graças a @whuber por sua ajuda até agora, não espero que você gaste mais tempo me ajudando, mas pensei em tentar esclarecer minha situação conforme você solicitou.

Oferecerei uma Rsolução que é codificada de maneira um pouco não Rilustrativa para ilustrar como ela pode ser abordada em outras plataformas.

A preocupação R(assim como em algumas outras plataformas, especialmente aquelas que favorecem um estilo de programação funcional) é que atualizar constantemente uma grande variedade pode ser muito caro. Em vez disso, esse algoritmo mantém sua própria estrutura de dados privada, na qual (a) todas as células que foram preenchidas até agora são listadas e (b) todas as células que estão disponíveis para serem escolhidas (em torno do perímetro das células preenchidas) estão listadas. Embora manipular essa estrutura de dados seja menos eficiente do que indexar diretamente em uma matriz, mantendo os dados modificados em um tamanho pequeno, provavelmente levará muito menos tempo de computação. (Também não foi feito nenhum esforço para otimizá-lo R. A pré-alocação dos vetores de estado deve economizar algum tempo de execução, se você preferir continuar trabalhando dentro R.)

O código é comentado e deve ser fácil de ler. Para tornar o algoritmo o mais completo possível, ele não usa nenhum complemento, exceto no final, para plotar o resultado. A única parte complicada é que, por eficiência e simplicidade, ele prefere indexar nas grades 2D usando índices 1D. Uma conversão acontece na neighborsfunção, que precisa da indexação 2D para descobrir quais podem ser os vizinhos acessíveis de uma célula e depois os converte no índice 1D. Essa conversão é padrão, portanto, não vou comentar mais, exceto para salientar que em outras plataformas GIS você pode querer reverter as funções dos índices de coluna e linha. (Em R, os índices de linha mudam antes dos índices da coluna.)

Para ilustrar, esse código pega uma grade que xrepresenta a terra e um recurso semelhante a um rio de pontos inacessíveis, começa em um local específico (5, 21) nessa grade (perto da curva inferior do rio) e a expande aleatoriamente para cobrir 250 pontos . O tempo total é de 0,03 segundos. (Quando o tamanho da matriz é aumentado em um fator de 10.000 a 3000 linhas por 5.000 colunas, o tempo sobe apenas para 0,09 segundos - um fator de apenas 3 ou mais - demonstrando a escalabilidade desse algoritmo.) Em vez de apenas produzindo uma grade de 0, 1 e 2, gera a sequência com a qual as novas células foram alocadas. Na figura, as células mais antigas são verdes, passando por douradas até as cores salmão.

Deve ser evidente que uma vizinhança de oito pontos de cada célula está sendo usada. Para outras vizinhanças, basta modificar o nbrhoodvalor próximo ao início de expand: é uma lista de compensações de índice em relação a qualquer célula. Por exemplo, um bairro "D4" pode ser especificado como matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

Também é evidente que esse método de propagação tem seus problemas: deixa buracos para trás. Se não é esse o objetivo, existem várias maneiras de corrigir esse problema. Por exemplo, mantenha as células disponíveis em uma fila para que as células mais antigas encontradas também sejam as primeiras preenchidas. Alguma randomização ainda pode ser aplicada, mas as células disponíveis não serão mais escolhidas com probabilidades uniformes (iguais). Outra maneira mais complicada seria selecionar as células disponíveis com probabilidades que dependem de quantos vizinhos cheios eles têm. Quando uma célula fica cercada, você pode ter uma chance de seleção tão alta que poucos furos seriam deixados sem preenchimento.

Termino comentando que esse não é um autômato celular (CA), que não procede célula por célula, mas atualiza faixas inteiras de células em cada geração. A diferença é sutil: com a CA, as probabilidades de seleção para células não seriam 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)

Com pequenas modificações, podemos fazer um loop expandpara criar vários clusters. É aconselhável diferenciar os clusters por um identificador, que aqui executará 2, 3, ..., etc.

Primeiro, mude expandpara retornar (a) NAna primeira linha se houver um erro e (b) os valores em indicesvez de em uma matriz y. (Não perca tempo criando uma nova matriz ya cada chamada.) Com essa alteração, o loop é fácil: escolha um início aleatório, tente expandi-lo, acumule os índices de cluster indicesse for bem-sucedido e repita até concluir. Uma parte importante do loop é limitar o número de iterações, caso muitos clusters contíguos não possam ser encontrados: isso é feito com count.max.

Aqui está um exemplo em que 60 centros de cluster são escolhidos uniformemente aleatoriamente.

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=" "))

Aqui está o resultado quando aplicado a uma grade de 310 por 500 (feita suficientemente pequena e grossa para que os clusters sejam aparentes). Demora dois segundos para executar; em uma grade de 3100 x 5000 (100 vezes maior), leva mais tempo (24 segundos), mas o tempo está escalando razoavelmente bem. (Em outras plataformas, como C ++, o tempo dificilmente deve depender do tamanho da grade.)

Sem nunca ter realizado sua operação e sem tempo de sobra para jogar, só posso adicionar esses dois links à sua lista:

Encontre o valor mais próximo da célula de varredura com base em um ponto vetorial (a primeira resposta (com 4 votos) é o que me intrigou).

Além disso: a Gridspread de Hawth ajudaria?