R/cut_raster.R

#' @title Cut raster based on a shape
#' @name cut.raster
#'
#' @description A function to cut a raster based on an informed shapefile.
#'
#' @param abio an object with the rasters to be cut. Accepts an object of type RasterStack, generated by the function \code{\link[raster]{stack}}
#' @param shape.dir here are two options. The first is to inform an object in the shapefile format. The second option is to inform a directory containing the shapefile.
#' @param extension extent of the cut rasters. Default is .asc, see\code{\link[raster]{write.raster}} for more output format possibilities.
#' @param br logical. If TRUE, use the Brazil shape provided by \code{\link[maptools]{wrld_simpl}}. If FALSE (default), use the informed shape.
#' @param plot logical. Plots one of the cut rasters.
#' @param trim logical. If TRUE, exclude pixels with NA (very slow). If FALSE (default), keep the NAs.
#'
#' @details When the shapeir argument is a folder where it contains the reference shape, this folder must contain only the shape and its auxiliary files.
#'
#' @return Raster files in a folder "Cortados".
#'
#' @author Diogo S. B. Rocha
#'
#' @seealso \code{\link[raster]{crop}}, \code{\link[raster]{mask}}, \code{\link[raster]{stack}}
#'
#' @examples
#'
#' fnames <- list.files(path=paste(system.file(package="dismo"), '/ex', sep=''), pattern='grd', full.names=TRUE )
#' predictors <- raster::stack(fnames)
#' cut.raster(abio = predictors, br=T, extension = ".tif")
#'
#' @import raster
#'
#' @export
cut.raster =
  function(abio, shape.dir, extension = ".asc", br = FALSE, plot = TRUE, diretorio = "cortados", trim = FALSE){

    if (dir.exists(paste0(diretorio)) == F) {
        dir.create(paste0(diretorio))
    }

    # Definir shape para cortar
    if(missing(shape.dir)){
      # Possibilidade de cortar para o Brasil
      #importando shape do brasil
      if(br==T){
        data(wrld_simpl, package = "maptools")
        shape.dir=subset(wrld_simpl, wrld_simpl$NAME=="Brazil")
        shape=shape.dir
      }else(stop("The directory containing the shape to cut the raster was not specified"))
      }

    if(class(shape.dir)== "character") {
      shape = rgdal::readOGR(list.files(shape.dir, pattern = ".shp", full.names = T)[1])
    } else (shape = shape.dir)

    crs(shape) = crs(abio)

    # loop para cortar todos os rasters

    # sem trim
    if (trim == F) {
        ini = Sys.time()
        for (i in 1:length(names(abio))) {
            ini1 = Sys.time()

            mask(crop(abio[[i]], extent(shape)), shape, filename = paste("./Cortados", "/", names(abio)[i],
                sep = "", extension), overwrite = TRUE)

            print(Sys.time())
            cat("\n", paste("Tá indo", i))

            fim1 = Sys.time()
            cat(paste("\n", round(as.numeric(fim1 - ini1), 2), units(fim1 - ini1)))
            if (i == length(names(abio))) {
                cat("\n", "Acabou!")
                fim = Sys.time()
                fim - ini
            }
        }
    }


    # com trim
    if (trim == T) {
        ini = Sys.time()
        for (i in 1:length(names(abio))) {
            ini1 = Sys.time()

            writeRaster(trim(mask(crop(abio[[i]], extent(shape)), shape)),
                filename = paste(".\\Cortados", "\\", names(abio)[i], ".tif", sep = ""),
                format = "GTiff", overwrite = TRUE, NAflag = -9999)
            cat("\r")
            print(Sys.time())
            cat("\n", paste("T? indo", i))
            fim1 = Sys.time()
            cat(paste("\n", round(as.numeric(fim1 - ini1), 2), units(fim1 - ini1)))
            if (i == length(names(abio))) {
                cat("\n", "Acabou!", "\n")
                fim = Sys.time()
                fim - ini
            }
        }
    }

    # plotando a primeira variavel cortada
    if (plot == T) {
        plot(raster(list.files("./cortados", pattern = "", full.names = TRUE)[1], native = TRUE))
        plot(shape, add = T)
    }
}
diogosbr/modelos documentation built on May 9, 2019, 5:23 p.m.