Zufällig geformte Zellklumpen in einem Raster aus Samen von 1 Zelle / Pixel erstellen?

Wie mein Titel schon sagt, möchte ich Zellklumpen aus Samen in einem Raster "züchten". Mein Basisraster ist voll mit Einsen und Nullen, Einsen zeigen Land und Nullen See / NA-Gebiete an. Aus den Einsen möchte ich 60 zufällige Pixel / Zellen als meine Samen auswählen und dann zufällig einen verbundenen Klumpen einer vordefinierten Nr. Wachsen lassen. Anzahl der Pixel / Zellen wird von diesem Startwert begrenzt. Ich habe gehört, dass die Technik als "Spread Dye" bezeichnet werden kann, hatte aber kein Glück, viel darauf zu finden. Die Keimzelle würde auf einen Wert von 2 gedreht und dann würde die nächste aus den umgebenden Einsen ausgewählte Zelle ebenfalls in 2 umgewandelt. 2 sind dann nicht verfügbar, um in Zukunft konvertiert zu werden.

Dieser Thread hilft ein wenig, da ich dies auch in R tun würde, da ich mit dem Lesen und Bearbeiten von GIS-Daten in R vertraut bin. Ich benötige jedoch eine Reihe von Regeln, um Pixel, die einen vorhandenen Klumpen umgeben, zufällig auszuwählen .

Wenn jemand diese grundlegendere Form von zellularen Automaten in einer GIS-Umgebung ausgeführt hat, würde ich mich über Ratschläge / Anleitungen freuen.

Beispiel:

Ich habe ein Ziel von 250 Zellen. Ich wähle zufällig eine Zelle mit einem Wert von 1 aus. Dies wird auf einen Wert von 2 gesetzt. Dann wird einer der Nachbarn der Startzelle, die = 1 ist, auf 2 gesetzt. Dann einer der Nachbarn einer der Zellen mit einem Wert von 2 wird ausgewählt und auf 2 gesetzt. Dies würde fortgesetzt, bis eine kontinuierliche Form mit einer Nummerierung von 250 Zellen erreicht ist.

Bearbeiten: Weitere Fragen

Aufgrund der großartigen Antwort von whuber habe ich einige Fragen zum Code:

1. Wie kann ich die Werte der gewachsenen Zellen nur einer '2' zuordnen und nicht variablen Werten, die die Reihenfolge darstellen, in der sie erstellt wurden?
2. Ich muss 60 Zellklumpen in meinem Bereich von '1' erstellen. Ich habe Wege gefunden, um zufällige Startpositionen auszuwählen, aber ich habe Mühe, alles innerhalb einer Schleife mit der von expandIhnen freundlicherweise geschriebenen Funktion zum Laufen zu bringen. Können Sie einen Weg vorschlagen, um 60 Klumpen zu erstellen, die nicht miteinander kollidieren und in derselben endgültigen Matrix enthalten sind?

Bearbeiten: Weitere Erklärung des Problems

Jeder Zellklumpen repräsentiert einen geschützten Bereich einer definierten Größe, z. B. 250 Zellen. Jedes Gebiet muss beginnen und zu Zellen mit einem Wert von 1 wachsen, da dies Land darstellt, und Zellen mit einem Wert von 0 vermeiden, da dies Meer darstellt. Ich muss dies 1000 Mal mit 60 geschützten Bereichen in jeder Iteration wiederholen, um ein Nullmodell zu erstellen, das zeigt, welche Verteilungen dieser Bereiche zufällig sein werden. Aus diesem Grund muss die Gesamtzahl der Zellen in allen 60 Bereichen in jeder der 1000 Iterationen identisch sein, damit sie vergleichbar sind. Daher ist es in Ordnung, wenn sich die Bereiche berühren, aber wenn es zu einer Kollision kommt, wächst der Klumpen idealerweise in eine andere verfügbare Richtung, bis das Ziel von 250 erreicht ist.

Sobald jedes dieser 1000 Schutzgebietsnetzwerke erstellt ist, werden sie als Maske gegen andere Rasterdaten wie Biodiversitätsmessungen verwendet, um festzustellen, (a) ob sie bestimmte Artenbereiche überschneiden und (b) wie viel Prozent bestimmter Artenbereiche diese zufälligen Netzwerke umfassen von Schutzgebieten abdecken.

Vielen Dank an @whuber für Ihre bisherige Hilfe. Ich erwarte nicht, dass Sie mehr Zeit damit verbringen, mir zu helfen, aber ich dachte, ich würde versuchen, meine Situation zu klären, wie Sie es gewünscht haben.

Ich werde eine RLösung anbieten , die leicht unkodiert ist R, um zu veranschaulichen, wie sie auf anderen Plattformen angegangen werden könnte.

Das Problem bei R(wie auch bei einigen anderen Plattformen, insbesondere solchen, die einen funktionalen Programmierstil bevorzugen) ist, dass die ständige Aktualisierung eines großen Arrays sehr teuer sein kann. Stattdessen behält dieser Algorithmus seine eigene private Datenstruktur bei, in der (a) alle bisher gefüllten Zellen aufgelistet sind und (b) alle zur Auswahl stehenden Zellen (um den Umfang der gefüllten Zellen) sind aufgelistet. Obwohl das Manipulieren dieser Datenstruktur weniger effizient ist als das direkte Indizieren in ein Array, wird es wahrscheinlich viel weniger Rechenzeit in Anspruch nehmen, wenn die geänderten Daten klein gehalten werden. (Es wurden auch keine Anstrengungen unternommen, um es zu optimieren R. Die Vorbelegung der Zustandsvektoren sollte einige Ausführungszeit sparen, wenn Sie lieber weiterarbeiten möchten R.)

Der Code ist kommentiert und sollte einfach zu lesen sein. Um den Algorithmus so vollständig wie möglich zu gestalten, werden keine Add-Ons verwendet, außer am Ende, um das Ergebnis zu zeichnen. Der einzige schwierige Teil ist, dass es aus Gründen der Effizienz und Einfachheit bevorzugt wird, mithilfe von 1D-Indizes in die 2D-Gitter zu indizieren. Eine Konvertierung erfolgt in der neighborsFunktion, die die 2D-Indizierung benötigt, um herauszufinden, welche Nachbarn einer Zelle zugänglich sein könnten, und konvertiert sie dann in den 1D-Index. Diese Konvertierung ist Standard, daher werde ich sie nicht weiter kommentieren, außer um darauf hinzuweisen, dass Sie in anderen GIS-Plattformen möglicherweise die Rollen von Spalten- und Zeilenindizes umkehren möchten. (In Rändern sich die Zeilenindizes vor den Spaltenindizes.)

Zur Veranschaulichung verwendet dieser Code ein Raster, xdas Land und ein flussähnliches Merkmal unzugänglicher Punkte darstellt, beginnt an einer bestimmten Stelle (5, 21) in diesem Raster (nahe der unteren Flussbiegung) und erweitert es zufällig auf 250 Punkte . Das Gesamt-Timing beträgt 0,03 Sekunden. (Wenn die Größe des Arrays um den Faktor 10.000 bis 3000 Zeilen um 5000 Spalten erhöht wird, steigt das Timing nur auf 0,09 Sekunden - ein Faktor von nur etwa 3 - und zeigt die Skalierbarkeit dieses Algorithmus.) Statt Wenn Sie nur ein Raster mit 0, 1 und 2 ausgeben, wird die Reihenfolge ausgegeben, mit der die neuen Zellen zugewiesen wurden. Auf der Figur sind die frühesten Zellen grün und gehen durch Gold in Lachsfarben über.

Es sollte offensichtlich sein, dass eine Acht-Punkte-Nachbarschaft jeder Zelle verwendet wird. Ändern Sie für andere Nachbarschaften einfach den nbrhoodWert am Anfang von expand: Es handelt sich um eine Liste von Indexversätzen relativ zu einer bestimmten Zelle. Zum Beispiel könnte eine "D4" -Nachbarschaft als angegeben werden matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

Es ist auch offensichtlich, dass diese Ausbreitungsmethode ihre Probleme hat: Sie hinterlässt Löcher. Wenn dies nicht beabsichtigt war, gibt es verschiedene Möglichkeiten, um dieses Problem zu beheben. Halten Sie beispielsweise die verfügbaren Zellen in einer Warteschlange, damit die frühesten gefundenen Zellen auch die frühesten ausgefüllten sind. Eine gewisse Randomisierung kann weiterhin angewendet werden, aber die verfügbaren Zellen werden nicht mehr mit einheitlichen (gleichen) Wahrscheinlichkeiten ausgewählt. Eine andere, kompliziertere Möglichkeit wäre die Auswahl verfügbarer Zellen mit Wahrscheinlichkeiten, die davon abhängen, wie viele gefüllte Nachbarn sie haben. Sobald eine Zelle umgeben ist, können Sie die Auswahlchance so hoch einstellen, dass nur wenige Löcher ungefüllt bleiben.

Abschließend möchte ich darauf hinweisen, dass dies kein zellularer Automat (CA) ist, der nicht zellenweise vorgeht, sondern ganze Zellmengen in jeder Generation aktualisiert. Der Unterschied ist subtil: Mit der CA wären die Auswahlwahrscheinlichkeiten für Zellen nicht einheitlich.

#
# 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)

Mit geringfügigen Änderungen können wir eine Schleife ausführen expand, um mehrere Cluster zu erstellen. Es ist ratsam, die Cluster durch eine Kennung zu unterscheiden, die hier 2, 3, ... usw. ausführt.

Ändern Sie expandzunächst (a) NAin der ersten Zeile, wenn ein Fehler vorliegt, und (b) die Werte indicesanstelle einer Matrix y. (Verschwenden Sie keine Zeit damit, ybei jedem Aufruf eine neue Matrix zu erstellen .) Mit dieser Änderung ist die Schleife einfach: Wählen Sie einen zufälligen Start, versuchen Sie, ihn zu erweitern, akkumulieren Sie die Cluster-Indizes, indiceswenn dies erfolgreich ist, und wiederholen Sie den Vorgang , bis Sie fertig sind. Ein wesentlicher Teil der Schleife besteht darin, die Anzahl der Iterationen zu begrenzen, falls viele zusammenhängende Cluster nicht gefunden werden können count.max.

Hier ist ein Beispiel, bei dem 60 Clusterzentren gleichmäßig zufällig ausgewählt werden.

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

Hier ist das Ergebnis, wenn es auf ein Raster von 310 x 500 angewendet wird (ausreichend klein und grob gemacht, damit die Cluster sichtbar werden). Die Ausführung dauert zwei Sekunden. Bei einem Raster von 3100 x 5000 (100-mal größer) dauert es länger (24 Sekunden), aber das Timing ist recht gut skalierbar. (Auf anderen Plattformen wie C ++ sollte das Timing kaum von der Rastergröße abhängen.)

Ohne jemals Ihre Operation durchgeführt zu haben und keine Zeit zum Spielen zu haben, kann ich nur diese beiden Links zu Ihrer Liste hinzufügen:

Finden Sie den nächstgelegenen Rasterzellenwert basierend auf einem Vektorpunkt (Die erste Antwort (mit 4 Stimmen) hat mich fasziniert).

Außerdem: Würde Hawths Gridspread helfen?