Membuat gumpalan sel berbentuk acak dalam raster dari biji 1 sel / piksel?

Seperti judul saya katakan, saya ingin "menumbuhkan" gumpalan sel dari biji dalam raster. Basis raster saya penuh dengan 1 dan 0, 1 menunjukkan tanah dan 0 laut / daerah NA. Dari angka 1, saya ingin memilih 60 piksel / sel acak sebagai benih saya, dan kemudian secara acak menumbuhkan rumpun terhubung dari no yang telah ditentukan sebelumnya. batas piksel / sel dari benih itu. Saya pernah mendengar teknik ini dapat disebut sebagai 'menyebar pewarna' tetapi belum beruntung menemukan banyak di atasnya. Sel benih akan berubah menjadi nilai 2 dan kemudian sel berikutnya yang dipilih dari sekitar 1 akan dikonversi menjadi 2 juga. 2's tidak tersedia untuk dikonversi di masa depan.

Utas ini memang sedikit membantu, karena saya juga bersedia melakukan ini di R karena saya terbiasa membaca dan memanipulasi data GIS di R. Namun, yang saya butuhkan adalah seperangkat aturan untuk secara acak memilih piksel di sekitar rumpun yang ada .

Jika ada yang telah melakukan bentuk automata seluler yang lebih mendasar ini dalam pengaturan GIS, maka saya sangat menghargai saran / bimbingan.

Contoh:

Saya memiliki target 250 sel. Saya secara acak memilih satu sel yang memiliki nilai 1. Ini berubah menjadi nilai 2. Kemudian, salah satu tetangga dari sel benih yang = 1 berubah menjadi 2. Kemudian, salah satu tetangga dari salah satu sel dengan nilai 2 dipilih dan diubah menjadi 2. Ini akan berlanjut sampai bentuk kontinu penomoran 250 sel tercapai.

Sunting: Pertanyaan Lebih Lanjut

Berdasarkan jawaban hebat whuber, saya punya beberapa pertanyaan tentang kode:

1. Bagaimana cara saya mengalokasikan nilai sel yang ditumbuhkan ke hanya '2' daripada nilai variabel yang mewakili urutan pembuatannya?
2. Saya perlu membuat 60 rumpun sel di dalam area '1's. Saya telah menemukan cara untuk memilih posisi awal yang acak tetapi berjuang untuk membuatnya semuanya bekerja dalam satu lingkaran menggunakan expandfungsi yang Anda tulis dengan baik. Bisakah Anda menyarankan cara untuk membuat 60 rumpun yang tidak saling berbenturan dan terkandung dalam matriks akhir yang sama?

Sunting: Penjelasan Masalah Lebih Lanjut

Setiap rumpun sel akan mewakili area terlindungi dengan ukuran yang ditentukan misalnya 250 sel. Setiap area harus mulai & tumbuh menjadi sel dengan nilai 1 karena ini mewakili daratan dan menghindari sel dengan nilai 0, karena ini mewakili laut. Saya perlu mengulang ini 1000 kali dengan 60 area yang dilindungi di setiap iterasi untuk membuat model nol, menunjukkan distribusi apa dari area-area ini secara kebetulan. Untuk alasan ini, jumlah total sel di seluruh 60 area harus identik di masing-masing dari 1000 iterasi sehingga mereka dapat dibandingkan. Oleh karena itu, tidak masalah jika area bersentuhan, tetapi jika ada tabrakan, maka idealnya rumpun akan tumbuh ke arah lain yang tersedia sampai target 250 tercapai.

Setelah masing-masing dari 1000 jaringan kawasan lindung ini dibuat, mereka akan digunakan sebagai topeng terhadap data raster lain seperti langkah-langkah keanekaragaman hayati untuk melihat (a) apakah mereka memotong rentang spesies tertentu dan (b) berapa% dari spesies tertentu yang rentang jaringan acak ini tutupan kawasan lindung.

Terima kasih kepada @whuber atas bantuan Anda sejauh ini, saya tidak berharap Anda menghabiskan lebih banyak waktu untuk membantu saya, tetapi saya pikir saya akan mencoba dan menjelaskan situasi saya seperti yang Anda minta.

Saya akan menawarkan Rsolusi yang dikodekan dengan sedikit Rcara untuk mengilustrasikan bagaimana hal itu dapat didekati di platform lain.

Kekhawatiran dalam R(serta beberapa platform lain, terutama yang mendukung gaya pemrograman fungsional) adalah bahwa memperbarui array besar secara terus-menerus bisa sangat mahal. Alih-alih, algoritma ini mempertahankan struktur data pribadinya sendiri di mana (a) semua sel yang telah diisi sejauh ini terdaftar dan (b) semua sel yang tersedia untuk dipilih (sekitar perimeter sel yang diisi) terdaftar. Meskipun memanipulasi struktur data ini kurang efisien daripada langsung mengindeks ke dalam array, dengan menjaga data yang dimodifikasi ke ukuran kecil, itu akan memakan waktu komputasi yang jauh lebih sedikit. (Tidak ada upaya telah dilakukan untuk mengoptimalkannya R, baik. Pra-alokasi vektor negara harus menghemat waktu eksekusi, jika Anda lebih memilih untuk tetap bekerja di dalam R.)

Kode ini dikomentari dan harus langsung dibaca. Untuk membuat algoritma selengkap mungkin, itu tidak menggunakan add-on apa pun kecuali di akhir untuk plot hasilnya. Satu-satunya bagian yang sulit adalah untuk efisiensi dan kesederhanaan itu lebih suka untuk mengindeks ke dalam kisi 2D dengan menggunakan indeks 1D. Konversi terjadi pada neighborsfungsi, yang membutuhkan pengindeksan 2D untuk mengetahui seperti apa tetangga yang dapat diakses dari sebuah sel dan kemudian mengubahnya menjadi indeks 1D. Konversi ini standar, jadi saya tidak akan berkomentar lebih lanjut kecuali untuk menunjukkan bahwa di platform GIS lain Anda mungkin ingin membalikkan peran kolom dan indeks baris. (Dalam R, indeks baris berubah sebelum indeks kolom dilakukan.)

Untuk mengilustrasikan, kode ini mengambil grid yang xmewakili tanah dan fitur seperti sungai dari titik-titik yang tidak dapat diakses, dimulai pada lokasi tertentu (5, 21) dalam grid itu (dekat tikungan sungai yang lebih rendah) dan meluaskannya secara acak untuk mencakup 250 poin . Total waktu adalah 0,03 detik. (Ketika ukuran array ditingkatkan dengan faktor 10.000 hingga 3.000 baris dengan 5.000 kolom, waktunya naik hanya menjadi 0,09 detik - faktor hanya 3 atau lebih - menunjukkan skalabilitas algoritma ini.) Alih-alih hanya mengeluarkan kisi-kisi 0's, 1's, dan 2's, ia mengeluarkan urutan dengan mana sel-sel baru dialokasikan. Pada gambar sel-sel awal berwarna hijau, lulus melalui emas menjadi warna salmon.

Harus jelas bahwa lingkungan delapan poin dari setiap sel sedang digunakan. Untuk lingkungan lain, cukup modifikasi nbrhoodnilai di dekat awal expand: ini adalah daftar offset indeks relatif terhadap sel yang diberikan. Misalnya, lingkungan "D4" dapat ditentukan sebagai matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

Juga jelas bahwa metode penyebaran ini memiliki masalah: ia meninggalkan lubang. Jika bukan itu yang dimaksudkan, ada berbagai cara untuk memperbaiki masalah ini. Misalnya, simpan sel yang tersedia dalam antrian agar sel yang paling awal ditemukan juga yang paling awal terisi. Beberapa pengacakan masih bisa diterapkan, tetapi sel yang tersedia tidak akan lagi dipilih dengan probabilitas yang seragam (sama). Cara lain yang lebih rumit adalah memilih sel yang tersedia dengan probabilitas yang bergantung pada berapa banyak tetangga yang diisi. Setelah sel dikepung, Anda bisa membuat peluang seleksi yang sangat tinggi sehingga beberapa lubang tidak terisi.

Saya akan selesai dengan berkomentar bahwa ini bukan automaton seluler (CA), yang tidak akan memproses sel demi sel, tetapi akan memperbarui seluruh petak sel dalam setiap generasi. Perbedaannya halus: dengan CA, probabilitas pemilihan untuk sel tidak akan seragam.

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

Dengan sedikit modifikasi, kami dapat mengulang expanduntuk membuat beberapa kluster. Dianjurkan untuk membedakan cluster dengan pengidentifikasi, yang di sini akan menjalankan 2, 3, ..., dll.

Pertama, ubah expanduntuk mengembalikan (a) NApada baris pertama jika ada kesalahan dan (b) nilai dalam indicesbukan matriks y. (Jangan buang waktu membuat matriks baru ydengan setiap panggilan.) Dengan perubahan ini dibuat, perulangan mudah: pilih mulai acak, cobalah untuk memperluas di sekitarnya, mengumpulkan indeks cluster indicesjika berhasil, dan ulangi sampai selesai. Bagian penting dari loop adalah untuk membatasi jumlah iterasi jika banyak cluster yang berdekatan tidak dapat ditemukan: ini dilakukan dengan count.max.

Berikut adalah contoh di mana 60 pusat cluster dipilih secara seragam secara acak.

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

Inilah hasilnya ketika diterapkan pada kisi 310 x 500 (dibuat cukup kecil dan kasar untuk kluster agar terlihat jelas). Diperlukan dua detik untuk mengeksekusi; pada 3100 x 5000 grid (100 kali lebih besar) dibutuhkan lebih lama (24 detik) tetapi waktunya scaling cukup baik. (Pada platform lain, seperti C ++, waktunya hampir tidak tergantung pada ukuran grid.)

Tanpa pernah melakukan operasi Anda, dan tidak ada waktu luang untuk bermain, saya hanya dapat menambahkan dua tautan ini ke daftar Anda:

Temukan Nilai Sel Raster Terdekat Berdasarkan Titik Vektor (Jawaban pertama (dengan 4 suara) adalah yang membuat saya penasaran).

Juga: apakah Hawth's Gridspread akan membantu?