Créer des amas de cellules de forme aléatoire dans un raster à partir de graines de 1 cellule / pixel?

Comme mon titre l'indique, je souhaite "faire pousser" des amas de cellules à partir de graines dans un raster. Mon raster de base est plein de 1 et de 0, 1 indique la terre et 0 la mer / NA. À partir des 1, je souhaite sélectionner 60 pixels / cellules aléatoires comme graines, puis développer aléatoirement un groupe connecté d'un non prédéfini. de pixels / cellules limite à partir de cette graine. J'ai entendu dire que la technique peut être qualifiée de «colorant à tartiner», mais je n'ai pas eu de chance d'en trouver beaucoup. La cellule de graine serait transformée à une valeur de 2, puis la cellule suivante choisie parmi les 1 environnants serait également convertie en 2. Les 2 ne sont alors plus disponibles pour être convertis à l'avenir.

Ce fil aide un peu, car je serais également disposé à le faire dans R car je suis familier avec la lecture et la manipulation des données SIG dans R. Cependant, ce dont j'ai besoin est un ensemble de règles pour sélectionner au hasard les pixels entourant un groupe existant .

Si quelqu'un a fait cette forme plus basique d'automates cellulaires dans un environnement SIG, j'apprécierais tout conseil / orientation.

Exemple:

J'ai un objectif de 250 cellules. Je sélectionne au hasard une cellule qui a une valeur de 1. Ceci est transformé en une valeur de 2. Ensuite, l'un des voisins de la cellule de semence qui = 1 est transformé en un 2. Ensuite, l'un des voisins de l'une ou l'autre des cellules avec une valeur 2 est sélectionnée et transformée en 2. Cela continuerait jusqu'à ce qu'une forme continue numérotée de 250 cellules soit atteinte.

Edit: Autres questions

Basé sur la grande réponse de whuber, j'ai quelques questions sur le code:

1. Comment puis-je attribuer les valeurs des cellules cultivées à un «2» plutôt qu'à des valeurs variables qui représentent l'ordre dans lequel elles ont été créées?
2. J'ai besoin de créer 60 amas de cellules dans ma zone de 1. J'ai imaginé des façons de choisir des positions de départ aléatoires, mais j'ai du mal à le faire fonctionner dans une boucle en utilisant la expandfonction que vous avez aimablement écrite. Pouvez-vous suggérer un moyen de créer 60 blocs qui ne se heurtent pas les uns aux autres et sont contenus dans la même matrice finale?

Edit: Explication supplémentaire du problème

Chaque groupe de cellules représentera une zone protégée d'une taille définie, par exemple 250 cellules. Chaque zone doit commencer et se développer en cellules avec une valeur de 1 car cela représente la terre et éviter les cellules avec une valeur de 0, car cela représente la mer. J'ai besoin d'itérer cette 1000 fois avec 60 zones protégées à chaque itération pour créer un modèle nul, montrant quelles distributions de ces zones seront par hasard. Pour cette raison, le nombre total de cellules dans les 60 zones devra être identique dans chacune des 1000 itérations afin qu'elles soient comparables. Par conséquent, c'est ok si les zones se touchent, mais s'il y a une collision, alors idéalement, la touffe se développerait dans une autre direction disponible jusqu'à ce que l'objectif de 250 soit atteint.

Une fois que chacun de ces 1000 réseaux d'aires protégées sera créé, ils seront utilisés comme masque contre d'autres données matricielles telles que les mesures de biodiversité pour voir (a) si elles croisent des gammes d'espèces particulières et (b) quel% d'espèces particulières range ces réseaux aléatoires des aires protégées couvrent.

Merci à @whuber pour votre aide jusqu'à présent, je ne m'attends pas à ce que vous passiez plus de temps à m'aider, mais j'ai pensé que j'essaierais de clarifier ma situation comme vous l'avez demandé.

Je proposerai une Rsolution qui est codée de manière légèrement non Rillustrée pour illustrer comment elle pourrait être abordée sur d'autres plateformes.

La préoccupation de R(ainsi que de certaines autres plates-formes, en particulier celles qui favorisent un style de programmation fonctionnel) est que la mise à jour constante d'un grand tableau peut être très coûteuse. Au lieu de cela, cet algorithme conserve sa propre structure de données privées dans laquelle (a) toutes les cellules qui ont été remplies jusqu'à présent sont répertoriées et (b) toutes les cellules qui peuvent être choisies (autour du périmètre des cellules remplies) sont répertoriés. Bien que la manipulation de cette structure de données soit moins efficace que l'indexation directe dans un tableau, en conservant les données modifiées à une petite taille, cela prendra probablement beaucoup moins de temps de calcul. (Aucun effort n'a été fait pour l'optimiser Rnon plus. La pré-allocation des vecteurs d'état devrait permettre de gagner du temps d'exécution, si vous préférez continuer à travailler à l'intérieur R.)

Le code est commenté et doit être simple à lire. Pour rendre l'algorithme aussi complet que possible, il n'utilise aucun module complémentaire, sauf à la fin pour tracer le résultat. La seule partie délicate est que pour l'efficacité et la simplicité, il préfère indexer dans les grilles 2D en utilisant des index 1D. Une conversion se produit dans la neighborsfonction, qui a besoin de l'indexation 2D afin de comprendre ce que pourraient être les voisins accessibles d'une cellule, puis les convertit en index 1D. Cette conversion est standard, donc je ne commenterai pas plus loin, sauf pour souligner que dans d'autres plates-formes SIG, vous voudrez peut-être inverser les rôles des index de colonnes et de lignes. (Dans R, les index de ligne changent avant les index de colonne.)

Pour illustrer, ce code prend une grille xreprésentant la terre et une caractéristique semblable à une rivière de points inaccessibles, commence à un emplacement spécifique (5, 21) dans cette grille (près du coude inférieur de la rivière) et l'étend au hasard pour couvrir 250 points . Le chronométrage total est de 0,03 seconde. (Lorsque la taille du tableau est augmentée d'un facteur de 10 000 à 3 000 lignes par 5 000 colonnes, le délai ne monte que jusqu'à 0,09 seconde - un facteur de 3 environ - ce qui démontre l'évolutivité de cet algorithme.) Au lieu de produisant juste une grille de 0, 1 et 2, il sort la séquence avec laquelle les nouvelles cellules ont été allouées. Sur la figure, les premières cellules sont vertes, passant de l'or à la couleur saumon.

Il doit être évident qu'un voisinage en huit points de chaque cellule est utilisé. Pour les autres quartiers, modifiez simplement la nbrhoodvaleur vers le début de expand: c'est une liste de décalages d'index par rapport à une cellule donnée. Par exemple, un voisinage "D4" pourrait être spécifié comme matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

Il est également évident que cette méthode d'épandage a ses problèmes: elle laisse des trous. Si ce n'est pas ce qui était prévu, il existe plusieurs façons de résoudre ce problème. Par exemple, conservez les cellules disponibles dans une file d'attente afin que les premières cellules trouvées soient également les plus anciennes remplies. Une certaine randomisation peut toujours être appliquée, mais les cellules disponibles ne seront plus choisies avec des probabilités uniformes (égales). Une autre façon, plus compliquée, serait de sélectionner les cellules disponibles avec des probabilités qui dépendent du nombre de voisins remplis qu'elles ont. Une fois qu'une cellule est entourée, vous pouvez augmenter ses chances de sélection de manière à ce qu'il reste peu de trous non remplis.

Je terminerai en commentant que ce n'est pas tout à fait un automate cellulaire (CA), qui ne procéderait pas cellule par cellule, mais mettrait à jour des pans entiers de cellules à chaque génération. La différence est subtile: avec l'AC, les probabilités de sélection des cellules ne seraient pas 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)

Avec de légères modifications, nous pouvons effectuer une boucle expandpour créer plusieurs clusters. Il convient de différencier les clusters par un identifiant, qui s'exécutera ici 2, 3, ..., etc.

Tout d'abord, changez expandpour retourner (a) NAà la première ligne s'il y a une erreur et (b) les valeurs dans indicesplutôt qu'une matrice y. (Ne perdez pas de temps à créer une nouvelle matrice yà chaque appel.) Une fois cette modification effectuée, le bouclage est facile: choisissez un démarrage aléatoire, essayez de vous développer autour de lui, accumulez les index de cluster en indicescas de succès et répétez jusqu'à ce que cela soit fait. Un élément clé de la boucle est de limiter le nombre d'itérations dans le cas où de nombreux clusters contigus ne peuvent pas être trouvés: cela se fait avec count.max.

Voici un exemple où 60 centres de grappes sont choisis uniformément au hasard.

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

Voici le résultat lorsqu'il est appliqué à une grille de 310 par 500 (rendue suffisamment petite et grossière pour que les grappes soient apparentes). L'exécution prend deux secondes; sur une grille de 3100 par 5000 (100 fois plus grande), cela prend plus de temps (24 secondes) mais le timing évolue assez bien. (Sur d'autres plateformes, telles que C ++, le timing ne devrait guère dépendre de la taille de la grille.)

Sans jamais avoir effectué votre opération, et sans avoir de temps à perdre pour jouer, je ne peux qu'ajouter ces deux liens à votre liste:

Trouvez la valeur de cellule raster la plus proche basée sur un point vectoriel (La première réponse (avec 4 votes) m'a intrigué).

Aussi: Hawth's Gridspread pourrait-il aider?