计算光栅砖的单元百分位数的最快方法 [R]答案

【问题标题】：Fastest way to calculate cell-wise percentiles for raster brick [R]计算光栅砖的单元百分位数的最快方法 [R]
【发布时间】：2019-09-21 03:59:31
【问题描述】：

我想为栅格砖的每个单元（+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 是一个栅格层（与“文件”完全相同），但只有一个层（不是砖）。

如果有任何关于如何加快速度的建议，我会很高兴。

谢谢！

【问题讨论】：

你能展示一个可复制的小例子供其他人测试吗
谢谢，我刚刚添加了一个可重现的光栅砖。请注意，我的砖块有更多的单元格。

标签： r loops statistics raster

【解决方案1】：

您可以通过并行运行上面的计算函数来加快速度。这个链接有一些关于如何的信息。祝你好运... https://www.gis-blog.com/increasing-the-speed-of-raster-processing-with-r-part-33-cluster/

【讨论】：

【解决方案2】：

避免循环遍历单元格，因为它太慢了。

示例数据

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

【讨论】：

感谢您的建议。这个过程仍然需要 6 个多小时，如果可能的话，我想进一步加快速度。