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

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

R - Speeding up replacement of values in large raster stacks

multicore

r

raster

geospatial