Realizar extracción por máscara en R haciendo uso de paralelización para optimizar el tiempo en el procesamiento

En este tutorial realizaremos una extracción por máscara de 48 archivos raster, que representan la precipitación, temperatura máxima, media y mínima, para 3 distintos departamentos de Colombia; la idea final es obtener los 48 raster para cada departamento, es decir que se obtendrían al final 144 raster. Para ello haremos uso de paralelización, lo que permite ejecutar el mismo proceso para los tres distintos departamentos en núcleos distintos, es decir, a la misma vez que se cortan los datos para Antioquia se están cortando los datos para Cundinamarca y Valle del Cauca, pues se hará uso de tres núcleos (si su computador tiene menos de cuatro núcleos no podrá replicar este ejercicio tal cual).

Los tres departamentos a los cuales les haremos la extracción por máscara son Antioquia, Valle del Cauca y Cundinamarca, para ello tendremos el shapefile de división político administrativa de Colombia (Fuente: SIGOT); además de ello, los raster climáticos para toda Colombia. El ejercicio incluye la descarga de los datos para que replique el ejercicio (URL).

1. Cargamos las librerías necesarias, en caso dado que usted no tenga la librería esta se descargara automáticamente.

if(!require(raster)){install.packages('raster'); library(raster)} else {library(raster)}
if(!require(rgdal)){install.packages('rgdal'); library(rgdal)} else {library(rgdal)}
if(!require(spdplyr)){install.packages('spdplyr'); library(spdplyr)} else {library(spdplyr)}

if(!require(tidyverse)){install.packages('tidyverse'); library(tidyverse)} else {library(tidyverse)}
if(!require(foreach)){install.packages('foreach'); library(foreach)} else {library(foreach)}
if(!require(parallel)){install.packages('parallel'); library(parallel)} else {library(parallel)}
if(!require(doSNOW)){install.packages('doSNOW'); library(doSNOW)} else {library(doSNOW)}
if(!require(gridExtra)){install.packages('gridExtra'); library(gridExtra)} else {library(gridExtra)}
if(!require(rasterVis)){install.packages('rasterVis'); library(rasterVis)} else {library(rasterVis)}

2. Descargamos y declaramos los archivos a utilizar en este tutorial, tanto los archivos raster como el shapefile de departamentos de Colombia

url    <- 'https://zdl9fa.bn1303.livefilestore.com/y4mV948RDpIKmA015QmtTAON46y0O8vW0kF7B3JQ6xeP006WJDddpGmwQ1wpRS9G3f36LGwVNy-EbwxqA9HPNUX6msaosaUNe3RJf8WynIjTbMohC490Mer7l5lrpEMTdVPR1-XWmPHv9CUaPxLFnqqu1qLIo5a36EjDeR4Ry0QmCsvt4G9TXyoLVoMkWt54VOUpAH4VIsLSuS5yk44wVI8Iw/_data.zip'

mydir  <- 'D:/parallel'
if(!file.exists(mydir)){dir.create(mydir)} 
temp   <- tempfile(tmpdir = mydir, fileext = '.zip')
download.file(url, temp)
unzip(temp, exdir = mydir)
unlink(temp) # Eliminar el archivo .zip

adm   <- shapefile(paste0(mydir, '/_data/_shp/LimiteDptal_geo_ok.shp'))
lyrs  <- list.files(paste0(mydir, '/_data/_raster/_co'), full.names = T, pattern = '.asc$') %>%
            lapply(FUN = raster) %>%
            stack()

3.  Del shape de Colombia extraemos tres objetos distintos que sean cada departamento, y luego los guardamos en un objeto tipo lista. 

ant   <- filter(adm, NOMBRE_DPT == 'ANTIOQUIA')
val   <- filter(adm, NOMBRE_DPT == 'VALLE DEL CAUCA')
cun   <- filter(adm, NOMBRE_DPT == 'CUNDINAMARCA')
shps  <- list(ant, val, cun)

4. Creamos las carpetas donde se guardaran los archivos cortados para cada departamento, en mi caso creare ‘_ant’ para Antioquía, ‘_val’ para Valle del Cauca, ‘_cun’ para Cundinamarca.

dirsOut <- paste0(path, '/', c('_ant', '_val', '_cun'))
if(!file.exists(dirsOut)){lapply(dirsOut, FUN = dir.create)}

5.  Indicamos la cantidad de núcleos a utilizar. Se recomienda utilizar un núcleo menos que el total de núcleos que tenga su computador.

nCores   <- detectCores(all.tests = FALSE, logical = TRUE)
cl       <- makeCluster(nCores - 1) #Número de nucleos a utilizar
registerDoSNOW(cl)

6. Configuramos la barra de progreso para así conocer que tanto hace falta para que se termine el proceso de extracción por máscara.

length_run <- length(shps)
pb         <- txtProgressBar(max = length_run, style = 3)
progress   <- function(n) setTxtProgressBar(pb, n) 
opts       <- list(progress=progress)

7. Creación de la función de extracción por máscara. Primero usamos la función, crop y mask para cortar por departamento y luego realizamos un unstack para tener los archivos raster como lista, y poderlos así escribir mediante un for sencillo.

extractMask <- function(lyrs, shps, dirsOut){
  
  lyrs_cut <- crop(lyrs, shps) %>%
                mask(shps) %>%
                unstack()

  for(j in 1:length(lyrs_cut)){
    
    writeRaster(lyrs_cut[[j]], paste0(dirsOut, '/', names(lyrs_cut[[j]]), '.asc'))
    
  }
  
  return(lyrs_cut)

}

8. Ejecutamos el foreach para correr el proceso paralelamente en los tres núcleos destinados para ello, se recomienda no tener abierto muchos programas en el computador para que así no se ponga lento el PC.

lyrs_cut_ok <- foreach(i = 1:length(shps), .packages = c('raster', 'rgdal', 'dplyr'), .options.snow=opts) %dopar% {
  
   extractMask(lyrs, shps[[i]], dirsOut[i])
  
} 

9. Creamos una función sencilla para plotear los archivos resultantes

plotGrid <- function(lyr){
  
  gg <- gplot(lyr) + 
            geom_tile(aes(fill = value)) + 
            coord_equal() +
            labs(fill = "Precipitación \n Enero (mm)") + 
            xlab('Lon') + 
            ylab('Lat')
  
}

plotAnt <- plotGrid(lyr = lyrs_cut_ok[[1]][[1]])
plotCun <- plotGrid(lyr = lyrs_cut_ok[[2]][[1]])
plotVal <- plotGrid(lyr = lyrs_cut_ok[[3]][[1]])

plotAll <- grid.arrange(plotAnt, plotCun, plotVal, ncol=3)

Aquí podemos observar la precipitación para el mes de enero en cada departamento:

Muchas gracias a Jeison Mesa por la ayuda brindada en la construcción de este tutorial.