Вычисление растровой статистики по последовательному списку растровых слоев

Question

Предположим, у меня есть rasterstack, состоящий из 5 слоев:

library(raster)
r <- raster(nrows=10,ncols=10)
r[] <- rnorm(100)
stk <- stack(r,r*2,r+10,r-7, r*6)

Я бы хотел вычислить статистику по последовательному набору слоев на основе другого списка чисел.

my.seq<-as.numeric(c("2", "1", "2"))

Как я могу вычислить статистику (например, среднее) для:

Слои 1: 2 в s, т.е. my.seq [1]

Уровень 3 в s, т.е. my.seq [2]

Слои 4: 5 в s, т.е. my.seq [3]

Я чувствую, что лучше всего это сделать с помощью цикла, но я не могу понять, как это будет работать.

Любая помощь будет принята с благодарностью.

---

EDIT

Это то, что я пытаюсь сделать IRL. Я подумал, что включу его сюда, чтобы помочь придать вопросу более широкий контекст.

Я пытаюсь вычислить среднемесячные значения на основе данных о суточной температуре. Растровый стек, с которым я работаю, имеет сотни слоев, по одному на каждый день месяца. Таким образом, слои 1:31, взятые вместе, представляют собой январь 1970 года. В приведенном выше примере my.seq представляет собой список, содержащий количество дней в каждом месяце. Поэтому я пытаюсь вычислить среднее значение для дней 1:31 (январь), затем 32:60 (февраль) и т. Д.

Answer 1

Эта функция использует my.seq для извлечения интересующих слоев (включая потенциально повторяющиеся слои, как в примере).

Вы можете поместить их в новый stack или вычислить / объединить статистику на уровне слоя в той части кода, которая приведена ниже.

new_stack <- function(x){
  tmp <- stack()
  for(i in x){
    new <- stk[[i]]
    tmp <- stack(tmp, new)  
  }
  return(tmp)
}

tmp2 <- new_stack(my.seq) 

Так вот, я не был уверен на 100%, как именно вы хотели вычислить среднее значение, поэтому я сделал это несколькими способами ниже:

stack_mean <- function(x){
  tmp         <- stack()
  layer_means <- numeric()

  for(i in x){
    new <- stk[[i]]

    layer_mean <- mean(stk[[i]]@data@values) 
    print(paste("In layer ", i, "the layer mean is ", layer_mean))

    layer_means <- c(layer_means, layer_mean)

    tmp <- stack(tmp, new)  

  }
  print(paste("Mean of means:", mean(layer_means)))
  return(tmp)
}

stack_mean(my.seq)

Это вычисляет несколько средств. Во-первых, средства слоев (которые я знаю не , что вы хотите - я просто показываю, как для потомков):

[1] "In layer  2 the layer mean is  -0.0981706786020096"
[1] "In layer  1 the layer mean is  -0.0490853393010048"
[1] "In layer  2 the layer mean is  -0.0981706786020096"

Тогда среднее из 3 слоев вместе:

[1] "Mean of means: -0.081808898835008"

Наконец, у вас есть среднее значение объекта stack, что, по-видимому, вещь, хотя я не очень много знаю об этом и понятия не имею, что все это на самом деле означает:

mean(tmp2) # equivalently: mean(stack_means(my.seq))

class       : RasterLayer 
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer 
values      : -5.742679, 4.206359  (min, max)

Answer 2

Вы можете использовать stackApply для этого

Пример данных

library(raster)
r <- raster(nrows=10, ncols=10, vals=rnorm(100))
stk <- stack(r,r*2,r+10,r-7, r*6)
my.seq <- c(2, 1, 2)

Получить индексы и использовать stackApply

i <- rep(1:length(my.seq), my.seq)
x <- stackApply(stk, i, fun='mean')

Answer 3

Я не понимаю ваш пример my.seq в отношении того, что вы объясняете в желании различных подмножеств растров. Вот обработанный пример, где я создаю объект списка с вашими определенными подмножествами. Все, что вам нужно сделать, это передать желаемый числовой индекс подмножества растров в двойной скобочный индекс объекта стека. Объект списка должен быть также подмножеством с использованием двойной скобки, поэтому он выглядит беспорядочно s[[idx[[i]]]]. Однако, если вы разобьете его, для первого набора растров это просто: raster::calc(s[[1:2]], mean)

library(raster)
file <- system.file("external/test.grd", package="raster")
s <- stack(file, file, file)
s <- addLayer(s, raster(file)/2, raster(file)*2)

idx <- list( c(1:2), c(3), c(4:5) )
s.mean <- stack()
  for(i in 1:length(idx)) {
    s.mean <- addLayer(s.mean, calc(s[[idx[[i]]]], mean) )
  }

s.mean
plot(s.mean)

Что касается вашего расширенного вопроса, вы можете использовать класс даты в R для создания индексов, а еще удобнее, пакет rts позволяет использовать эти типы сводок временных рядов.

Здесь давайте создадим стек с 365 слоями, представляющими дни года.

library(raster)
f <- raster(nrows=50, ncols=50, xmn=0, xmx=10)
s <- stack()
  for(i in 1:365) {
    x <- f
    x[] <- runif(ncell(x),0,10)
    s <- addLayer(s, x)
}

Теперь мы можем создать соответствующий вектор дат.

( d <- seq(as.Date("1970/1/1"), as.Date("1970/12/31"), "days") ) 

Вектор дат можно запросить для получения индекса растра. Скажем, мы хотим получить среднее значение за декабрь, мы можем использовать which для получения индекса.

( dec.idx <- which( months(d) == "December") )
( dec.mean <- calc(s[[dec.idx]], mean) ) 

Вы можете легко создать список, содержащий индексы растрового стека для каждого месяца, который будет эквивалентен тому, что вы описываете в своем объекте my.seq.

months.idx <- list()   
  for(m in unique( months(d) ) ) {
    months.idx[[m]] <- which( months(d) == m)  
  }  
months.idx

Тем не менее, это встроенная функциональность для выполнения такого типа временных суммирований в пакете rts, поэтому мы можем сократить любой цикл for, приведя стек к объекту rts, а затем используя одну из функций apply.

library(rts)
s.date <- rts::rts(s, d)
( s.month <- apply.monthly(s.date, mean) )