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! 



Comentarios

  1. Excelente, aunque este ejemplo es algo mas simplificado, el video esta todo muy claro.
    Muchas gracias de verdad me ha sido de muchísima ayuda

    ResponderEliminar
  2. Una pregunta:
    Si deseo guardar los raster cortados como .tif o pasarlos .asc solo cambiando en el código:
    writeRaster(corte,paste("E:/Blogger/Post_7/_data/_prec_valle/",names(rasters[[i]]),".tif",sep=""))

    .tif x .asc. Puede funcionar bien??

    ResponderEliminar
    Respuestas
    1. Solamente cambias el .tif por el formato que quieres, en el caso de ascii sería .asc

      Saludos desde Cali - Colombia

      Eliminar
  3. Hola! Gracias por compartir. Seguí el script y me resulto genial en el corte, pero tuve un problema, que las capas no encajan una con otra. ¿Cómo puedo asegurar el snap? esto para que los pixeles concuerden entre todas las capas?

    ResponderEliminar
    Respuestas
    1. un parámetro de la función crop es snap, y allí pones la capa con la cual se alinee

      Eliminar

Publicar un comentario

Entradas populares de este blog

Convertir una tabla a shape en R

¿Cómo descargar información de GBIF usando R?