Terra R -使用自定义函数加速栅格数据的aggregate()

Question

我想使用terraR包中的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时使用自定义函数?

Answer 1

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

对于那些考虑foreach并行化的人来说，没有简单的解决方案可以使用terra对象运行foreach循环。对于这样的任务，我仍然使用旧的栅格包。这是一个有计划的更新，但不是在短期内-看到这里。详情如下。

<标题>玩具集

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秒

Withtdigest::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()+ parallelization

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()+ parallelization

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)

时差7.520231秒

简而言之:

1 tdigest 1.922526 secs

2 base_quantile 6.013551 secs

3 tdigest_parallel 7.520231 secs

4 base_quantile_parallel 8.537751 secs

Answer 2

这并不是说aggregate()在使用自定义函数时本身很慢。相反，使用quantile()而不是median()来获得中位数的成本更高。这大概是由于计算本身的成本(terra使用c++实现来计算中位数，比任意分位数更快)，也因为quantile()执行更多的检查，因此在此过程中调用了更多的附加函数。当像aggregate那样多次执行操作时，这一较高的计算成本就会增加。

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

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