Terra R - 使用自定义函数加速栅格数据的聚合（）

Posted 2023-03-23

技术标签:

【中文标题】Terra R - 使用自定义函数加速栅格数据的聚合（）【英文标题】：Terra R - Speed up aggregate() of raster data with custom function 【发布时间】：2021-10-22 03:15:11 【问题描述】：

我想使用 terra R 包中的 aggregate 函数以分位数方法作为聚合函数聚合栅格。下面，我使用来自R base 的quantile 函数使用本地包目录中的栅格计算第50 个百分位数（即中位数）。我选择了第 50 个百分位数与中位数进行比较，但我的目标确实是计算其他分位数...

library(terra)
# load elevation coming with the terra pakage
r <- rast( system.file("ex/elev.tif", package="terra") )

plot(r)

# number of iteration
n_it <- 20

# with a custom function
start_time <- Sys.time()
for (i in 1:n_it)
ra <- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T))

end_time <- Sys.time()

我的电脑大约花了。 6秒做20次。

print(end_time-start_time)

时差 6.052727 秒

当我使用中值内置函数运行相同的 aggregate 运行时，它需要大约。执行相同的 20 次迭代的时间减少了 40 倍！

# with a built-in function
start_time <- Sys.time()
for (i in 1:n_it)
  ra <- aggregate(r, 2 , fun = median)

end_time <- Sys.time()
print(end_time-start_time)

0.1456101秒的时间差

由于我想计算第 50 位以外的其他百分位数，有人可以提供一些建议以加快 aggregate 在使用自定义函数时的速度吗？

【问题讨论】：

谷歌"r fast quantile function". 【参考方案1】：

在使用自定义函数时，aggregate() 本身并不慢。相反，使用quantile() 而不是median() 来获得中位数会更昂贵。这可能是由于计算本身的成本（terra 使用比arbitrary quantile 更快的C++ implementation to compute the median），还因为quantile() 执行更多检查，因此在此过程中调用更多附加函数。当aggregate 多次执行该操作时，这种较高的计算成本就会增加。

如果您有一个更大的栅格，则使用 cores 参数在多个内核上分布计算可能是有益的，请参阅 ?terra::aggregate。但是，我认为这不是 elev 数据的选项，因为开销太大。

如果您想为许多不同的probs 调用aggregate，您可以并行化循环，例如使用foreach package。

【讨论】：

谢谢。我认为 foreach 循环并不是 terra 的真正选择。或者至少，使用 raster 包执行此操作要简单得多。在此处查看更多信息：***.com/questions/67445883/…【参考方案2】：

根据回复，我测试了两个选项：使用tdigest package 和terra 包中的内置并行化例程（cores 参数）。 Dunning et al., (2019) 的 t-Digest 构造算法使用一维 k-means 聚类的变体来生成非常紧凑的数据结构，可以准确估计分位数。我建议使用tquantile 函数，它可以将测试数据集的处理时间减少三分之一。

对于那些正在考虑foreach 并行化的人来说，没有简单的解决方案可以使用 terra 对象运行 foreach 循环。对于此类任务，我仍在使用良好的旧光栅包。 It is a planned update but not on short term - see here。更多详情如下。

玩具数据集

library(terra)
# load elevation coming with the terra pakage
r <- rast( system.file("ex/elev.tif", package="terra") )

plot(r)

# number of iteration
n_it <- 20

# With `stats::quantile()` function
start_time <- Sys.time()
for (i in 1:n_it)
  ra <- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T))

end_time <- Sys.time()
print(end_time-start_time)

时差 6.013551 秒

与`tdigest::tquantile()`

library(tdigest)

start_time_tdigest <- Sys.time()
for (i in 1:n_it)
  ra_tdigest <- aggregate(r, 2 , fun = function(x) tquantile(tdigest(na.omit(x)), probs = .5))
  
end_time_tdigest <- Sys.time()
print(end_time_tdigest-start_time_tdigest)

时差1.922526秒

正如 Martin 所怀疑的那样，在 terra:aggregate 函数中使用 cores 参数并没有提高处理时间：

`stats::quantile()` + 并行化

start_time_parallel <- Sys.time() 
for (i in 1:n_it)
ra_tdigest_parallel <- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T), cores = 2) 
 
end_time_parallel <- Sys.time() 
print(end_time_parallel-start_time_parallel)

时差 8.537751 秒

`tdigest::tquantile()` + 并行化

tdigest_quantil_terra <- function(x)    
require(tdigest)   
tquantile(tdigest(na.omit(x)), probs = .5) 
        
start_time_tdigest_parallel <- Sys.time() for (i in 1:n_it)
ra_tdigest_parallel <- aggregate(r, 2 , 
fun = function(x, ff) ff(x), cores = 2 , 
ff = tdigest_quantil_terra)     
 

end_time_tdigest_parallel <- Sys.time() 
print(end_time_tdigest_parallel-start_time_tdigest_parallel)