Tworzenie losowo ukształtowanych skupisk komórek w rastrze z nasion 1 komórki / piksel?

Jak mówi mój tytuł, chcę „wyhodować” skupiska komórek z nasion w rastrze. Mój podstawowy raster jest pełen 1 i 0, 1 oznacza ląd i 0 morze / obszary NA. Z 1-tych chciałbym wybrać 60 losowych pikseli / komórek jako swoje nasiona, a następnie losowo wyhodować połączoną kępę o wstępnie zdefiniowanym nr. limitu pikseli / komórek z tego materiału siewnego. Słyszałem, że tę technikę można nazwać „barwnikiem do smarowania”, ale nie miałem szczęścia znaleźć jej dużo. Komórka zaszczepiająca zostanie zamieniona na wartość 2, a następnie następna komórka wybrana z otaczających jedynek również zostanie przekonwertowana na 2. 2 są następnie niedostępne do konwersji w przyszłości.

Ten wątek trochę pomaga, ponieważ chętnie bym to zrobił w R, ponieważ znam się na czytaniu i manipulowaniu danymi GIS w R. Jednak wymagam zestawu reguł do losowego wybierania pikseli otaczających istniejącą grupę .

Jeśli ktoś zrobiłby tę bardziej podstawową formę automatów komórkowych w środowisku GIS, byłbym wdzięczny za wszelkie porady / wskazówki.

Przykład:

Mam cel 250 komórek. Losowo wybieram jedną komórkę, która ma wartość 1. To jest zamieniane na wartość 2. Następnie jeden z sąsiadów komórki zaszczepiającej, który = 1 jest zamieniany na 2. Następnie jeden z sąsiadów jednej z komórek z wartością 2 jest wybierana i zamieniana na 2. Trwało to do momentu osiągnięcia ciągłego kształtu liczącego 250 komórek.

Edycja: dalsze pytania

Na podstawie doskonałej odpowiedzi Whubera mam kilka pytań dotyczących kodu:

1. Jak mogę przydzielić wartości hodowanych komórek do wartości „2” zamiast wartości zmiennych reprezentujących kolejność, w której zostały utworzone?
2. Muszę utworzyć 60 skupisk komórek w moim obszarze „1”. Wymyśliłem sposoby wybierania losowych pozycji początkowych, ale staram się, aby wszystko działało w pętli za pomocą expandfunkcji, którą uprzejmie napisałeś. Czy możesz zasugerować sposób utworzenia 60 grup, które nie kolidują ze sobą i są zawarte w tej samej ostatecznej matrycy?

Edycja: Dalsze wyjaśnienie problemu

Każde skupisko komórek będzie reprezentowało chroniony obszar o określonym rozmiarze, np. 250 komórek. Każdy obszar musi zacząć i rosnąć w komórki o wartości 1, ponieważ reprezentuje to ląd, i unikać komórek o wartości 0, ponieważ reprezentuje morze. Muszę iterować to 1000 razy z 60 chronionymi obszarami w każdej iteracji, aby stworzyć model zerowy, pokazujący, jakie rozkłady tych obszarów będą przypadkowe. Z tego powodu łączna liczba komórek we wszystkich 60 obszarach będzie musiała być identyczna w każdej z 1000 iteracji, aby były porównywalne. Dlatego jest ok, jeśli obszary się stykają, ale jeśli dojdzie do kolizji, idealnie kępa będzie rosła w innym dostępnym kierunku, aż do osiągnięcia celu 250.

Po utworzeniu każdej z tych 1000 sieci obszarów chronionych będą one używane jako maska ​​w stosunku do innych danych rastrowych, takich jak miary różnorodności biologicznej, aby zobaczyć (a) czy przecinają one poszczególne zakresy gatunków oraz (b) jaki procent poszczególnych gatunków obejmuje te losowe sieci ochrony obszarów chronionych.

Dzięki @whuber za pomoc do tej pory, nie oczekuję, że poświęcisz mi więcej czasu, ale pomyślałem, że postaram się wyjaśnić moją sytuację, tak jak prosiłeś.

Oferuję Rrozwiązanie, które jest nieco zakodowane w Rcelu zilustrowania, jak można podejść do niego na innych platformach.

Problemem R(podobnie jak niektórych innych platform, szczególnie tych, które preferują funkcjonalny styl programowania) jest to, że ciągła aktualizacja dużej tablicy może być bardzo droga. Zamiast tego algorytm zachowuje własną strukturę prywatnych danych, w której (a) wymienione są wszystkie komórki, które zostały dotychczas wypełnione, i (b) wszystkie komórki, które można wybrać (na obwodzie wypełnionych komórek) są wymienione. Mimo że manipulowanie tą strukturą danych jest mniej wydajne niż bezpośrednie indeksowanie do tablicy, utrzymując zmodyfikowane dane w małym rozmiarze, prawdopodobnie zajmie to znacznie mniej czasu obliczeniowego. (Nie podjęto również żadnych wysiłków, aby go zoptymalizować R. Wstępna alokacja wektorów stanu powinna zaoszczędzić trochę czasu wykonania, jeśli wolisz nadal pracować wewnątrz R).

Kod jest komentowany i powinien być łatwy do odczytania. Aby algorytm był jak najbardziej kompletny, nie wykorzystuje żadnych dodatków poza końcem do wykreślenia wyniku. Jedyną trudną częścią jest to, że dla wydajności i prostoty woli indeksować do siatek 2D za pomocą indeksów 1D. Konwersja zachodzi w neighborsfunkcji, która potrzebuje indeksowania 2D, aby dowiedzieć się, czym mogą być dostępni sąsiedzi komórki, a następnie przekształca je w indeks 1D. Ta konwersja jest standardowa, więc nie będę komentować jej dalej, poza zaznaczeniem, że na innych platformach GIS możesz chcieć odwrócić role indeksów kolumn i wierszy. (W Rindeksy wierszy zmieniają się przed indeksami kolumn).

Aby to zilustrować, ten kod pobiera siatkę xprzedstawiającą ląd i rzekę cechę niedostępnych punktów, rozpoczyna się w określonej lokalizacji (5, 21) na tej siatce (w pobliżu dolnego zakola rzeki) i rozwija ją losowo, aby objąć 250 punktów . Całkowity czas wynosi 0,03 sekundy. (Gdy rozmiar tablicy zostanie zwiększony o współczynnik od 10 000 do 3000 wierszy o 5000 kolumn, czas wzrośnie tylko do 0,09 sekundy - współczynnik tylko około 3 - co pokazuje skalowalność tego algorytmu.) Zamiast po prostu wyprowadzając siatkę zer, 1 i 2, wyświetla sekwencję, z którą przydzielono nowe komórki. Na rycinie najwcześniejsze komórki są zielone, przechodząc przez złoto do kolorów łososiowych.

Powinno być oczywiste, że stosuje się ośmiopunktowe sąsiedztwo każdej komórki. W przypadku innych dzielnic wystarczy zmodyfikować nbrhoodwartość na początku expand: jest to lista przesunięć indeksu względem dowolnej komórki. Na przykład sąsiedztwo „D4” można określić jako matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

Oczywiste jest również, że ta metoda rozprzestrzeniania ma swoje problemy: pozostawia dziury w tyle. Jeśli nie było to zamierzone, istnieją różne sposoby rozwiązania tego problemu. Na przykład trzymaj dostępne komórki w kolejce, aby najwcześniej znalezione komórki były również najwcześniej wypełnione. Nadal można zastosować pewną randomizację, ale dostępne komórki nie będą już wybierane z jednakowymi (równymi) prawdopodobieństwami. Innym, bardziej skomplikowanym sposobem byłby wybór dostępnych komórek z prawdopodobieństwem, które zależą od liczby ich wypełnionych sąsiadów. Gdy komórka zostanie otoczona, jej szansa na wybranie jest tak wysoka, że ​​niewiele otworów pozostanie niewypełnionych.

Zakończę komentując, że nie jest to całkiem automat komórkowy (CA), który nie przechodziłby komórka po komórce, ale zamiast tego aktualizowałby całe warstwy komórek w każdym pokoleniu. Różnica jest subtelna: w przypadku CA prawdopodobieństwo wyboru komórek nie byłoby jednolite.

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

Przy niewielkich modyfikacjach możemy zapętlić się, expandaby utworzyć wiele klastrów. Wskazane jest rozróżnienie klastrów według identyfikatora, który tutaj będzie działał 2, 3, ... itd.

Najpierw zmień, expandaby zwracać (a) NAw pierwszym wierszu, jeśli wystąpił błąd, i (b) wartości indiceszamiast w macierzy y. (Nie marnuj czasu na tworzenie nowej matrycy yprzy każdym wywołaniu.) Po wprowadzeniu tej zmiany zapętlenie jest łatwe: wybierz losowy początek, spróbuj rozszerzyć się wokół niego, akumuluj indeksy klastrów, indicesjeśli się powiedzie, i powtarzaj, aż się zakończy. Kluczową częścią pętli jest ograniczenie liczby iteracji w przypadku, gdy nie można znaleźć wielu sąsiadujących klastrów: odbywa się to za pomocą count.max.

Oto przykład, w którym 60 centrów skupień jest wybieranych losowo równomiernie.

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

Oto wynik po nałożeniu na siatkę 310 na 500 (wykonany wystarczająco mały i gruby, aby klastry były widoczne). Wykonanie zajmuje dwie sekundy; na siatce 3100 na 5000 (100 razy większa) zajmuje więcej czasu (24 sekundy), ale czas jest odpowiednio skalowany. (Na innych platformach, takich jak C ++, taktowanie nie powinno zależeć od rozmiaru siatki.)

Nie wykonując żadnej operacji i nie tracąc czasu na grę, mogę dodać tylko te dwa linki do Twojej listy:

Znajdź najbliższą wartość komórki rastrowej na podstawie punktu wektorowego (zaintrygowała mnie pierwsza odpowiedź (4 głosy)).

Ponadto: czy Hawth's Gridspread pomoże?