R - 加速替换大型栅格堆栈中的值

Question

最近我遇到了一些麻烦，这些代码正在加速我用来为某些卫星数据应用云遮罩的代码。

有关大型栅格数据加速替换值的提示？

Answer 1

以下代码设置两个大型栅格堆栈，并使用一个掩盖另一个栅格的值。请注意，这些测试是在32核心服务器上完成的。使用foreach的改进高度依赖于高性能计算机的使用。

library(foreach) 
library(doParallel) 
library(raster) 

xy = matrix(rnorm(1000^2),1000,1000) 
# Turn the matrix into a raster 
rast = raster(xy) 
# Give it lat/lon coords for 36-37°E, 3-2°S 
extent(rast) = c(36,37,-3,-2) 
# ... and assign a projection 
projection(rast) = CRS("+proj=longlat +datum=WGS84") 

# create first random raster stack 
raster_list = list() 
for(i in 1:300){ 
    set.seed(i) 
    rast[]= matrix(rnorm(1000^2),1000,1000) 
    raster_list = c(raster_list,rast) 
} 

rast_stack = stack(raster_list) 
backup_raster_stack = rast_stack # store a backup to reset 
plot(rast_stack[[3]]) 

# create second random raster stack 
raster_list2 = list() 
for(i in 1:300){ 
    set.seed(i+25) 
    rast[]= matrix(rnorm(1000^2),1000,1000) 
    raster_list2 = c(raster_list2,rast) 
} 

rast_stack2 = stack(raster_list2) 
plot(rast_stack2[[3]]) 


################## 
# do value replacement the normal way 
ptm <- proc.time() 
rast_stack[rast_stack2>1]=NA 
proc.time() - ptm 

用户系统经过

426.872 23.734 462.304

# reset values 
raster_stack = backup_raster_stack 

################# 
# do value replacement in calc 
ptm <- proc.time() 
calc(rast_stack , function(x) { rast_stack[ rast_stack2 > 1 ] = NA; return(x) }) 
proc.time() - ptm 

用户系统经过

491.732 30.242 530.230

# reset values 
raster_stack = backup_raster_stack 

################# 
# do value replacement in foreach loop 
registerDoParallel(32) 
ptm <- proc.time() 
foreach(i=1:dim(rast_stack)[3]) %do% { rast_stack[[i]][rast_stack2[[i]]>1]=NA} 
proc.time() - ptm 

用户系统经过

57.654 1.692 59.378