Terra R - 使用自定义函数加速栅格数据的聚合（）
Terra R - Speed up aggregate() of raster data with custom function
我想使用 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)
Time difference of 6.052727 secs
当我运行同样aggregate运行用median内置函数时，用了大约。执行完全相同的 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)
Time difference of 0.1456101 secs
因为我想计算除第 50 个以外的其他百分位数，有人可以提供一些建议以在使用自定义函数时加快 aggregate 速度吗？
并不是说 aggregate() 使用自定义函数本身就慢。相反，使用 quantile() 而不是 median() 来获取中位数的成本更高。这大概是由于计算本身的成本（terra 使用 C++ implementation to compute the median that is faster than for an arbitrary quantile），还因为 quantile() 执行了更多的检查，因此在过程中调用了更多的附加函数。当 aggregate.
多次执行该操作时，这种更高的计算成本会增加
如果您有更大的栅格，则使用 cores 参数在多个内核之间分配计算可能会有所帮助，请参阅 ?terra::aggregate。但是，我认为这不是 elev 数据的选项，因为开销太大。
如果你想为许多不同的 probs 调用 aggregate 你可以并行化循环，例如使用 foreach package.
根据回复，我测试了两个选项：使用 tdigest package and the built-in parallelization routine from terra package (cores parameter). The t-Digest construction algorithm, by Dunning et al., (2019)，使用一维 k 均值聚类的变体来生成非常紧凑的数据结构，可以准确估计分位数。我建议使用 tquantile 函数，它可以将测试数据集的处理时间减少三分之一。
对于那些考虑 foreach 并行化的人来说，运行 foreach 循环与 terra 对象没有简单的解决方案。对于此类任务，我仍在使用良好的旧光栅包。 。下面有更多详细信息。
玩具数据集
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)
Time difference of 6.013551 secs
和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)
Time difference of 1.922526 secs
正如 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)
Time difference of 8.537751 secs
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)
Time difference of 7.520231 secs
简而言之：
1 tdigest 1.922526 secs
2 base_quantile 6.013551 secs
3 tdigest_parallel 7.520231 secs
4 base_quantile_parallel 8.537751 secs
我想使用 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)
Time difference of 6.052727 secs
当我运行同样aggregate运行用median内置函数时，用了大约。执行完全相同的 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)
Time difference of 0.1456101 secs
因为我想计算除第 50 个以外的其他百分位数，有人可以提供一些建议以在使用自定义函数时加快 aggregate 速度吗？
并不是说 aggregate() 使用自定义函数本身就慢。相反，使用 quantile() 而不是 median() 来获取中位数的成本更高。这大概是由于计算本身的成本（terra 使用 C++ implementation to compute the median that is faster than for an arbitrary quantile），还因为 quantile() 执行了更多的检查，因此在过程中调用了更多的附加函数。当 aggregate.
如果您有更大的栅格，则使用 cores 参数在多个内核之间分配计算可能会有所帮助，请参阅 ?terra::aggregate。但是，我认为这不是 elev 数据的选项，因为开销太大。
如果你想为许多不同的 probs 调用 aggregate 你可以并行化循环，例如使用 foreach package.
根据回复，我测试了两个选项：使用 tdigest package and the built-in parallelization routine from terra package (cores parameter). The t-Digest construction algorithm, by Dunning et al., (2019)，使用一维 k 均值聚类的变体来生成非常紧凑的数据结构，可以准确估计分位数。我建议使用 tquantile 函数，它可以将测试数据集的处理时间减少三分之一。
对于那些考虑 foreach 并行化的人来说，运行 foreach 循环与 terra 对象没有简单的解决方案。对于此类任务，我仍在使用良好的旧光栅包。
玩具数据集
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)
Time difference of 6.013551 secs
和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)
Time difference of 1.922526 secs
正如 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)
Time difference of 8.537751 secs
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)
Time difference of 7.520231 secs
简而言之：
1 tdigest 1.922526 secs
2 base_quantile 6.013551 secs
3 tdigest_parallel 7.520231 secs
4 base_quantile_parallel 8.537751 secs