Самый быстрый способ вычисления клеточных процентилей для растрового кирпича [R]

Question

Я хотел бы рассчитать несколько статистических данных для каждой ячейки растрового кирпича (+50 слоев, 70 000 000 ячеек), например, медиана.Этот маленький кусочек кода работает, но невероятно медленно ...:

    r <- raster(ncol=3000, nrow=2000)
    r[] <- rnorm(n=ncell(r))
    s <- stack(x=c(r, r*2, r))
    file <- brick(x=c(r, r*2, r))
    q25 <- file[[1]]
    i=0
    for(i in 1:ncell(file)){
      progress(i, max.value = ncell(file))
      q25[i] <- unname(quantile(file[i],.5,na.rm=TRUE))
    }

q25 - растровый слой (точно такой же, как у 'file'), но только один слой (не кирпич).

Буду рад любым предложениям, как ускорить процесс.

Спасибо!

Answer 1

Избегайте зацикливания на ячейках, поскольку оно слишком медленное.

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

library(raster)
r <- raster(ncol=300, nrow=200)
set.seed(1)
values(r) <- rnorm(n=ncell(r))
s <- stack(r, r*2, r)

Ваш подход

q25 <- raster(s)
for(i in 1:ncell(q25)){
   q25[i] <- quantile(s[i], .5, na.rm=TRUE)
}

лучший подход (~ в 10 раз быстрее, см. Ниже)

q1 <- calc(s, fun=function(x) quantile(x, .5, na.rm=TRUE))

Но для медианы быстрее (~ 37x) использовать

q2 <- calc(s, fun=function(x)median(x, na.rm=TRUE))

Сравнение скорости:

system.time(for(i in 1:ncell(q25)){ q25[i] <- quantile(s[i], .5, na.rm=TRUE)  })
#user  system elapsed 
#69.24   13.90   83.71 
system.time( q1 <- calc(s, fun=function(x)quantile(x, .5, na.rm=TRUE)))
#user  system elapsed 
#7.81    0.00    7.81 
system.time(q2 <- calc(s, fun=function(x)median(x, na.rm=TRUE)))
#user  system elapsed 
#2.23    0.00    2.23 

Answer 2

Вы можете ускорить процесс, запустив функции calc выше, параллельно.Эта ссылка имеет некоторую информацию о том, как.Удачи ... https://www.gis -blog.com / увеличение скорости обработки растра с помощью r-part-33-cluster /