Создание беспорядочно сформированных скоплений ячеек в растре из семян размером 1 ячейка / пиксель?

Как говорится в моем названии, я хочу «вырастить» скопления клеток из семян в растре. Мой базовый растр полон 1 и 0, 1 указывает на сушу, а 0 - на море / NA. Из числа 1 я хочу выбрать 60 случайных пикселей / ячеек в качестве своих семян, а затем случайным образом вырастить связную группу заранее определенного числа нет. пикселей / ячеек предел от этого семени. Я слышал, что эту технику можно назвать «спред-красителем», но мне не очень-то повезло. Начальная ячейка будет установлена ​​в значение 2, а затем следующая ячейка, выбранная из окружающих 1, будет также преобразована в 2. 2 тогда недоступны для преобразования в будущем.

Этот поток немного помогает, поскольку я также хотел бы сделать это в R, так как я знаком с чтением и манипулированием данными ГИС в R. Однако мне требуется набор правил для случайного выбора пикселей, окружающих существующий скопление ,

Если кто-то сделал эту более базовую форму клеточных автоматов в условиях ГИС, я был бы признателен за любые советы / рекомендации.

Пример:

У меня цель 250 клеток. Я случайным образом выбираю одну ячейку, которая имеет значение 1. Это превращается в значение 2. Затем один из соседей исходной ячейки = 1 превращается в 2. Затем один из соседей любой из ячеек со значением 2 выбирается и превращается в 2. Это будет продолжаться до тех пор, пока не будет достигнута непрерывная форма, насчитывающая 250 ячеек.

Изменить: дополнительные вопросы

Основываясь на отличном ответе Уубер, у меня есть несколько вопросов о коде:

1. Как я могу распределить значения выращенных клеток только на «2», а не на значения переменных, которые представляют порядок, в котором они были созданы?
2. Мне нужно создать 60 скоплений клеток в моей области 1. Я разработал способы выбора случайных стартовых позиций, но изо всех сил стараюсь заставить все это работать в цикле, используя expandфункцию, которую вы любезно написали. Можете ли вы предложить способ создания 60 скоплений, которые не конфликтуют друг с другом и содержатся в одной и той же конечной матрице?

Изменить: дальнейшее объяснение проблемы

Каждый скопление клеток будет представлять собой защищенную область определенного размера, например, 250 клеток. Каждая область должна начинаться и расти в ячейки со значением 1, поскольку это представляет землю и избегать ячеек со значением 0, поскольку это представляет море. Мне нужно повторить эту 1000 раз с 60 защищенными областями в каждой итерации, чтобы создать нулевую модель, показывающую, какие распределения этих областей будут случайными. По этой причине общее количество ячеек во всех 60 областях должно быть одинаковым в каждой из 1000 итераций, чтобы они были сопоставимы. Следовательно, это нормально, если области касаются, но если есть столкновение, то в идеале сгусток будет расти в другом доступном направлении, пока не будет достигнута цель 250.

Как только каждая из этих 1000 сетей охраняемых районов будет создана, они будут использоваться в качестве маски против других растровых данных, таких как меры по биоразнообразию, чтобы увидеть (а) пересекают ли они ареалы определенных видов и (б) какой процент конкретных видов охватывает эти случайные сети покрытия охраняемых территорий.

Спасибо @whuber за вашу помощь, я не ожидаю, что вы потратите больше времени, помогая мне, но подумал, что постараюсь прояснить ситуацию, как вы просили.

Я предложу Rрешение, которое написано немного не Rтак, чтобы проиллюстрировать, как оно может подходить на других платформах.

Проблема R(как и некоторых других платформ, особенно тех, которые предпочитают стиль функционального программирования) заключается в том, что постоянное обновление большого массива может быть очень дорогим. Вместо этого этот алгоритм поддерживает свою собственную структуру частных данных, в которой (а) перечислены все ячейки, которые были заполнены до сих пор, и (б) все ячейки, которые доступны для выбора (по периметру заполненных ячеек) перечислены. Хотя манипулирование этой структурой данных менее эффективно, чем прямое индексирование в массив, при сохранении измененных данных небольшого размера это, вероятно, займет гораздо меньше времени вычислений. (Также не было предпринято никаких усилий для его оптимизации R. Предварительное распределение векторов состояния должно сэкономить некоторое время выполнения, если вы предпочитаете продолжать работать внутри R.)

Код комментируется и должен быть простым для чтения. Чтобы сделать алгоритм максимально полным, он не использует никаких дополнений, кроме как в конце, чтобы отобразить результат. Единственная сложность заключается в том, что для эффективности и простоты он предпочитает индексировать в двумерные сетки, используя одномерные индексы. Преобразование происходит в neighborsфункции, которая нуждается в двумерном индексировании, чтобы выяснить, какими могут быть доступные соседи ячейки, а затем преобразует их в одномерный индекс. Это преобразование является стандартным, поэтому я не буду комментировать его далее, за исключением того, что хочу указать, что на других платформах ГИС вы можете изменить роли индексов столбцов и строк. (In R, индексы строк изменяются раньше, чем индексы столбцов.)

Чтобы проиллюстрировать, этот код берет сетку, xпредставляющую сушу и речную особенность недоступных точек, начинается в определенном месте (5, 21) в этой сетке (около нижнего изгиба реки) и расширяет ее случайным образом, чтобы покрыть 250 точек , Общее время составляет 0,03 секунды. (Когда размер массива увеличивается в 10 000–3 000 строк на 5000 столбцов, время увеличивается только до 0,09 секунды - в 3 раза или более - демонстрируя масштабируемость этого алгоритма.) Вместо просто выводя сетку из 0, 1 и 2, она выводит последовательность, с которой были выделены новые ячейки. На рисунке самые ранние ячейки зеленого цвета, переходящие из золота в цвета лосося.

Должно быть очевидно, что используется восьмиточечная окрестность каждой ячейки. Для других окрестностей просто измените nbrhoodзначение в начале expand: это список смещений индекса относительно любой данной ячейки. Например, район "D4" может быть указан как matrix(c(-1,0, 1,0, 0,-1, 0,1), nrow=2).

Также очевидно, что у этого метода распространения есть свои проблемы: он оставляет дыры позади. Если это не то, что было задумано, существуют различные способы решения этой проблемы. Например, сохраняйте доступные ячейки в очереди, чтобы самые ранние найденные ячейки были также заполнены самыми ранними. Некоторая рандомизация все еще может быть применена, но доступные ячейки больше не будут выбираться с одинаковыми (равными) вероятностями. Другим, более сложным способом будет выбор доступных ячеек с вероятностями, которые зависят от того, сколько у них заполненных соседей. Как только клетка окажется в окружении, вы можете сделать так, чтобы ее выбор был настолько велик, что несколько отверстий останутся незаполненными.

В заключение я прокомментирую, что это не совсем клеточный автомат (CA), который не будет обрабатывать ячейку за ячейкой, а вместо этого будет обновлять целые полосы ячеек в каждом поколении. Разница тонкая: с CA вероятности выбора для ячеек не будут одинаковыми.

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

С небольшими изменениями мы можем зацикливаться, expandчтобы создать несколько кластеров. Желательно дифференцировать кластеры по идентификатору, который здесь будет работать 2, 3, ... и т. Д.

Во-первых, измените expandна return (a) NAв первой строке, если есть ошибка, и (b) значения в indicesматрице y. (Не тратьте время на создание новой матрицы yпри каждом вызове.) После внесения этого изменения цикл становится легким: выберите случайный старт, попробуйте расширить его, накапливайте кластерные индексы, indicesесли все прошло успешно, и повторяйте до тех пор, пока не закончите. Ключевой частью цикла является ограничение количества итераций в случае, если невозможно найти много смежных кластеров: это делается с помощью count.max.

Вот пример, где 60 кластерных центров выбираются случайным образом равномерно.

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

Вот результат при применении к сетке 310 на 500 (сделан достаточно малым и грубым, чтобы кластеры были заметны). Выполнение занимает две секунды; на сетке 3100 на 5000 (в 100 раз больше) это занимает больше времени (24 секунды), но время масштабируется достаточно хорошо. (На других платформах, таких как C ++, время вряд ли должно зависеть от размера сетки.)

Даже не выполнив свою операцию и не пожалев времени на игру, я могу добавить только две ссылки в ваш список:

Найти значение ближайшей растровой ячейки на основе векторной точки (первый ответ (с 4 голосами) меня заинтриговал).

Кроме того: поможет ли Gridspread от Hawth ?