R/spClassifyRast.R

In tfrescino/FIESTA: Forest Inventory Estimation and Analysis

#' Data - Reclass raster.
#'
#' Wrapper to reclass a raster using a vector of cut breaks.
#'
#'
#' @param rastfn String. Path name of raster to classify.
#' @param cutbreaks Integer vector. Breaks to use for classifying (e.g.,
#' c(0,50,75) uses function in calc: 'ifelse (A >= 0 & A < 50, 1, 
#' ifelse (A >= 50 & A < 75, 2, ifelse (A >= 75, 3, 255)))'
#' @param bnd sf R object or String. Boundary to clip raster (optional).
#' Can be a spatial sf object, full pathname to a shapefile, or name of
#' a layer within a database.
#' @param bnd_dsn String. Name of data source name with bnd_layer, if in a
#' database.
#' @param bnd.filter String. Optional filter of bnd_layer.
#' @param buffdist Number. The distance to buffer the polygon before clipping
#' raster, in units of raster.
#' @param savedata_opts List. See help(savedata_options()) for a list
#' of options. Only used when savedata = TRUE.  
#' @return Data.
#' @author Tracey S. Frescino
#' @keywords data
#' @export
spClassifyRast <- function(rastfn,
                           cutbreaks,
                           bnd = NULL, 
                           bnd_dsn = NULL, 
                           bnd.filter = NULL, 
						   buffdist = NULL,
                           savedata_opts = NULL) {


  ## Set global variables
  gui <- FALSE
  savedata <- TRUE
  savetmp <- TRUE

  ##################################################################
  ## CHECK PARAMETER NAMES
  ##################################################################
  input.params <- names(as.list(match.call()))[-1]
  formallst <- names(formals(spClassifyRast))
  if (!all(input.params %in% formallst)) {
    miss <- input.params[!input.params %in% formallst]
    stop("invalid parameter: ", toString(miss))
  }

  ## Check parameter lists
  pcheck.params(input.params, savedata_opts=savedata_opts)


  ## Set savedata defaults
  savedata_defaults_list <- formals(savedata_options)[-length(formals(savedata_options))]
  
  for (i in 1:length(savedata_defaults_list)) {
    assign(names(savedata_defaults_list)[[i]], savedata_defaults_list[[i]])
  }
  
  ## Set user-supplied savedata values
  if (length(savedata_opts) > 0) {
    if (!savedata) {
      message("savedata=FALSE with savedata parameters... no data are saved")
    }
    for (i in 1:length(savedata_opts)) {
      assign(names(savedata_opts)[[i]], savedata_opts[[i]])
    }
  }

  ## Check overwrite, outfn.date, outfolder, outfn
  ########################################################
  if (savedata) {
    outlst <- pcheck.output(out_dsn=out_dsn, out_fmt=out_fmt, 
                  outfolder=outfolder, outfn.pre=outfn.pre, 
                  outfn.date=outfn.date, overwrite_dsn=overwrite_dsn, 
                  overwrite_layer=overwrite_layer, 
                  append_layer=append_layer, createSQLite=FALSE, gui=gui)
    out_dsn <- outlst$out_dsn
    outfolder <- outlst$outfolder
    out_fmt <- outlst$out_fmt
    overwrite_layer <- outlst$overwrite_layer
    overwrite_dsn <- outlst$overwrite_dsn
    append_layer <- outlst$append_layer
  } 
  if (is.null(out_layer)) {
    out_layer <- "rastcl"
  }

  ## Verify ref raster
  ########################################################
  rast <- getrastlst.rgdal(rastfn)

  ## Get crs of reference raster
  rast.crs <- rasterInfo(rast)$crs
  rast.nodata <- rasterInfo(rast)$nodata_value
  rast.datatype <- rasterInfo(rast)$datatype


  ## Check cutbreaks
  if (!is.vector(cutbreaks) || !is.numeric(cutbreaks)) {
    stop("cutbreaks must be a numeric vector")
    if (!identical(sort(cutbreaks),cutbreaks)) {
	  stop("cutbreaks must be in sequential order")
	}
  }
 
  ####################################################################
  ## Get boundary info
  ## If a boundary is input, clip the ref raster to the boundary mask
  ####################################################################

  ## Import boundary
  bndx <- pcheck.spatial(layer=bnd, dsn=bnd_dsn, caption="boundary")
  if (!is.null(bndx)) {
    bndx <- datFilter(bndx, xfilter=bnd.filter, stopifnull=TRUE)$xf

    ## Compare crs of reference raster and reproject if different
    bndx <- crsCompare(bndx, rast.crs)$x	

    ## Check buffdist
    if (!is.null(buffdist)) {
      if (!is.numeric(buffdist)) stop("invalid buffdist... must be numeric")
    }
	
    if (!is.null(buffdist)) {
      ## This will buffer the polygon 1 pixel to include all pixels inside boundary
      bndx <- sf::st_buffer(bndx, dist=buffdist)
    }
 	
					  
	## Clip reference raster to polygon boundary 
    ###########################################################################
	if (savetmp) {
	  cliprastfn <- file.path(outfolder, "clip_rast.tif")
	} else {
	  cliprastfn <- tempfile("clip_rast", fileext="tif")
    }
	
	## Get nodata value of ref raster.
	## If the nodata values was not assigned (NA), then use a default value.
	if (is.na(rast.nodata)) {
	  nodata <- getDefaultNodata(rast.datatype)
	} else {
	  nodata <- rast.nodata
	}
	  
   ####################################################################
   ## Clip raster
   ####################################################################
   clipRaster(src = bndx, 
             srcfile = rast, 
			 src_band = 1, 
			 dstfile = cliprastfn,			 
             fmt = "GTiff", 
			 init = nodata,
			 dstnodata = nodata,
			 maskByPolygons = TRUE, 
             options = "COMPRESS=LZW")
  } else {
  
    cliprastfn <- rast
  }
  
  ## Create expression based on cut breaks
  expr <- ""
  for (i in 1:(length(cutbreaks))) {
    if (i == length(cutbreaks)) {
	  expr <- paste0(expr, "ifelse (A >= ", cutbreaks[i], ", ", i, ", ")
	} else {
      expr <- paste0(expr, "ifelse (A >= ", cutbreaks[i], 
					" & A < ", cutbreaks[i + 1], ", ", i, ", ") 
	}
  }
  expr <- paste0(expr, "255", paste(rep(")", length(cutbreaks)), collapse=""))  
  message(expr)
  
  lut <- data.frame(MIN = cutbreaks, 
                    MAX = c(cutbreaks[-1], paste0(cutbreaks[length(cutbreaks)], "+")), 
					CLASS = seq(1:length(cutbreaks))) 
  
  
  ####################################################################
  ## Reclass raster
  ####################################################################
  outfn <- paste0(file.path(outfolder, out_layer), ".tif")
  gdalraster::calc(expr, 
	               rasterfiles = cliprastfn,
				   var.names = "A",
				   dstfile = outfn,
				   dtName = "Byte", 
				   fmt = "GTiff",
				   options = c("COMPRESS=DEFLATE"),
                   nodata_value = 255,
				   setRasterNodataValue = TRUE,
				   write_mode = "overwrite"
				   )


   ####################################################################
   ## Save lookup table to outfolder
   ####################################################################
   datExportData(lut, 
          savedata_opts = list(outfolder = outfolder, 
		                       out_fmt = "csv", 
							   out_layer = paste0(out_layer, "_lut")))


  return(list(outfn=outfn, lut=lut))
}