R - Accélérer le remplacement des valeurs dans les grandes piles raster

Question

Récemment, j'avais des difficultés à accélérer le code que j'utilisais pour appliquer un masque de nuage pour certaines données satellite.

Des conseils pour accélérer le remplacement des valeurs dans les grandes matrices?

Answer 1

Le code suivant définit deux grandes piles raster et en utilise une pour masquer les valeurs d'une autre. Veuillez noter que ces tests ont été effectués sur un serveur 32 core. Les améliorations de l'utilisation de foreach dépendent fortement de l'utilisation d'ordinateurs haute performance.

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 

système utilisateur

écoulé

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 

système utilisateur écoulé

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 

système utilisateur

écoulé

57,654 1,692 59,378