Как разделить выборку растра на основе вероятностей на три категории и задать значения этих категорий в R?

Question

Прошу прощения за громоздкий вопрос. Я пытаюсь создать растр, в котором случайный 1% пикселей слоя имеет значения, и из этих оцененных пикселей случайные 35% имеют значение 1, 55% имеют значение 2 и 10% имеют значение 3. остальные пиксели должны занимать маркер «нет данных» в R («NA»).

Достаточно просто создать растр с 1% пикселей с однородным значением, используя следующий код:

pixels <- raster(ext = extent(-120, -119, 49, 50), resolution = c(0.001, 0.001), crs = CRS("+proj=longlat +datum=WGS84"), vals = 1)

testing <- sampleRandom(pixels, size = as.integer(0.01*ncell(pixels)), asRaster = TRUE)

Однако я не уверен, как затем разделить оцененные пиксели testing на три категории и установите значения этих категорий в соответствии с приведенным выше.

Возможно ли это или еще один способ достичь того, что я хочу?

Answer 1

Это не слишком громоздко,

Функция f принимает общее количество ячеек

f <- function(N) {
    n <- N/100  # 1% sample
    # create a vector with the values you want
    v <- c(rep(1, 0.35*n), rep(2, 0.55*n), rep(3, 0.1*n))
    # sample these values (that is, put them in random order)
    v <- sample(v)
    # create output vector 
    out <- rep(NA, N)
    # put the values in random places
    out[sample(N, length(v))] <- v
    out
}

library(raster)
# create a RasterLayer
r <- raster(ncol=100, nrow=100)
set.seed(10) # for reproducibility
values(r) <- f(ncell(r))

Показать, что это работает

table(values(r))
# 1  2  3 
#35 55 10 

Другой способ: быть использовать вероятности с выборкой

set.seed(10) 
N <- ncell(r)
v <- sample(3, N/100, prob=c(0.35, 0.55, 0.1), replace=TRUE)
table(v)
# 1  2  3 
#30 67  3 

Но поскольку здесь используются вероятности, пропорции не точные . В данном случае это кажется далеким, но это из-за небольшого размера выборки.

Следуя вашему примеру, вы могли бы также go по этому маршруту

library(raster)
r <- raster(ext=extent(-120, -119, 49, 50), resolution=c(0.001, 0.001), crs="+proj=longlat +datum=WGS84", vals = 1)
r <- sampleRandom(r, size = (0.01*ncell(r)), asRaster = TRUE)

sfun <- function(x) {
   i <- !is.na(x)
   x[i] <- sample(1:3, sum(i), prob=c(0.35, 0.55, 0.1), replace=TRUE)
   x
}

set.seed(101)
x <- calc(r, sfun)

Опять примерно правильно

tab <- table(values(x))
100 * tab / sum(tab)
#   1     2     3 
#35.45 54.62  9.93