Extracción por mascara en R

Extracción por máscara en R

El objetivo del presente ejercicio es realizar la extracción por mascara (corte) de un grupo de archivos raster a partir de un shapefile, este usado como máscara, mediante un ciclo en R, es un proceso común cuando se tienen, por ejemplo, la precipitación (o alguna otra variable) para los 12 meses del año de todo Colombia y se quiere obtener solamente la precipitación de una unidad administrativa de menor área, para el presente caso es el departamento del Valle del Cauca; evitando hacer el proceso repetitivo, por ejemplo, en ArcGIS usando la herramienta “Extract by Mask”. Se hace uso de los ciclos, básicamente consiste en decirle al sistema que repita en cada una de las capas de un listado el corte a partir de una capa que representa el límite de nuestra área de estudio.

 Fuente: esri

Primero habilitamos las librerías de las cuales haremos uso:

require(raster) #si la librería no la tiene incluida se puede descargar usando install.packages("raster")

require(rgdal)

Luego abrimos el shapefile, que representa la máscara, bajo el comando “shapefile”.

mascara = shapefile("E:/Blogger/Post_7/_data/_mask/valle_cacuca.shp")

Para poder visualizar ambos archivos hacemos uso del comando “plot”.

plot(mascara)

class(mascara) #permite ver el tipo de archivo que es mascara, para este caso es un shapefile tipo polígono.

Luego hacemos el uso del comando “paste” para indicarle a R que lea todos los archivos que contienen la dirección escrita, esto permite unir la dirección, el 1 a 12 significa cada mes (es como finaliza el archivo) y por último la extensión que es .asc.

precipitacion=paste("E:/Blogger/Post_7/_data/_prec_col/prec_",1:12,".asc", sep="")

Seguido de ello hacemos uso del comando “lapply” para indicar que todos los archivos son raster.

rasters=lapply(precipitacion, FUN = raster)

Luego realizamos el ciclo, aspectos a tener en cuenta:

·         El for es el comando que se utiliza

·         La “i” es para indicar que cada vez que esta letra cambie, conforme a la secuencia (1:12) será un archivo distinto al cual se le hace el proceso.

·         Lo que está dentro de las “{}” será lo que repetirá una y otra vez hasta finalizar.

·         El comando “crop” se usa para realizar el corte (equivalente al “extract by mask” de ArcGIS), en la primera posición va el dato de entrada, y en la segunda posición va el dato de máscara crop(rasters, mask),

·         El comando “mask” se usa para indicar que la máscara (extensión del archivo) será equivalente al objeto “mascara”.

·         El comando writeRaster es para guardar los archivos en una dirección de nuestro computador, y el comando “names" se utiliza para que guarde el archivo bajo el mismo nombre que el archivo de entrada.

for(i in 1:19){

  corte = crop(rasters[[i]],mascara)

  corte = mask(corte, mascara)

  writeRaster(corte,paste("E:/Blogger/Post_7/_data/_prec_valle/",names(rasters[[i]]),".tif",sep=""))

}

plot(corte[[1]])

De este link se puede descargar los datos que se hacen uso en este tutorial.
Gracias a @AntonioPantojaC por su colaboración!