Функция суммирования каждой ячейки сетки растрового стека с использованием других растров в качестве индикатора

## input raster
s <- stack(list.files("~/dailyraster", full.names=TRUE)) # daily raster stack
r_start <- raster("~/stackSumSTART.asc") # this raster contain starting Julian day
r_end <- raster("~/stackSumEND.asc") # this raster contain ending Julian day
noNAcells <- which(!is.na(r[])) # cell numbers which contain values

## dummy raster
x <- r
x[] <- NA 

## loop
for (i in noNAcells) {      
  x[i] <- sum(s[[r_start[i]:r_end[i]]][i])
}

Я хотел бы создать функцию типа stackApply(), но я хочу, чтобы она работала на основе ячеек.

Выше представлена ​​for() версия цикла, и она работает хорошо, но занимает слишком много времени.

Дело в том, что каждая ячейка получает диапазон sum() из двух растровых слоев, r_start, r_end в приведенном выше скрипте.

Сейчас я пытаюсь преобразовать этот код с помощью семейства apply().

Есть ли возможность улучшить скорость с помощью for() loop? или дайте мне несколько советов по написанию этого кода на apply()

Любые комментарии мне помогут, спасибо.


r raster
person Kongsol    schedule 17.04.2018    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Ваш подход

x <- s$layer.1

system.time(
for (i in 1:ncell(x)) {
     x[i] <- sum(s[[r_start[i]:r_end[i]]][i], na.rm = T)
   }
)

   user  system elapsed 
  0.708   0.000   0.710

Мое предложение

Вы можете добавить растры, используемые в качестве индексов, в конец вашего стека, а затем использовать calc, чтобы значительно ускорить процесс (~ 30-50x).

s2 <- stack(s, r_start, r_end)
sum_time <- function(x) {sum(x[x[6]:x[7]], na.rm = T)}

system.time(
   output <- calc(s2, fun = sum_time)
)

   user  system elapsed 
  0.016   0.000   0.015

all.equal(x, output)

[1] TRUE

Образец данных

library(raster)

# Generate rasters of random values
r1 <- r2 <- r3 <- r4 <- r5 <- r_start <- r_end <- raster(ncol=10, nrow=10)

r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
r4[] <- rnorm(ncell(r4), 1, 0.2)
r5[] <- rnorm(ncell(r5), 1, 0.2)
s <- stack(r1,r2,r3,r4,r5)

r_start[] <- sample(1:2, ncell(r_start),replace = T)
r_end[] <- sample(3:5, ncell(r_end),replace = T)
person DJack    schedule 17.04.2018
comment
Большое вам спасибо, DJack! - person Kongsol; 17.04.2018
comment
У меня есть следующие вопросы. Я изменил sum_time() как sum_time2 <- function(x) {sum(x[x[ifelse(nlayers(x)==7, 6, 7)]:x[ifelse(nlayers(x)==7, 7, 8)]], na.rm = T)}, но консоль возвращает ошибку cannot use this function. Когда я использую nlayers(s2), кроме nlayers(x) в функции ifelse() (которая немного утомлена), это работает. Я думаю, что есть какое-то правило для определения или использования функции, которую я не знаю. - person Kongsol; 18.04.2018
comment
Внутри функции x - это вектор, соответствующий значениям каждого растра в стеке для одной ячейки. Это означает, что вы должны использовать length(x) вместо nlayers(x). При этом добавление этой функции (немного) уменьшит скорость calc. Лучше сохранить a = nlayers(s2) - 1 и b = nlayers(s2) перед вызовом calc и определить функцию как sum(x[x[a]:x[b]]. - person DJack; 18.04.2018