R/proximity.R

In juoe/spatialtoolbox: Collection of functions for doing spatial stuff.

#' Calculate distance raster via gdal_proximity using RQGIS.
#'
#' @param x Raster object or path to raster layer.
#' @param values Numeric indicating target values for distance computation. Default is all non-NA values.
#' @param filename Optional filename for writing output. If no filename is supplied, a temporary file will be created.
#' @param units Compute distance in pixels (0) or geographic units of x (1)
#' @param in_meters Whether to force computation in meters. Default is FALSE.
#' @return Distance raster
#' @export

proximity <- function(x, values = NULL, units = 1, in_meters = FALSE, filename = NULL) {

  if (is.null(filename)) {
    filename <- tempfile(fileext='.tif')
  }

  if (is.null(values)){
    x <- raster::mask(x, is.na(x), inverse = TRUE, updatevalue = 1)
    values <- 1
  }

  if (in_meters){

    in_crs <- raster::crs(x)

    web_mercator <- raster::crs("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs")

    # maybe find faster alternative
    x <- raster::projectRaster(x, crs = web_mercator)

    units <- 1

    out <- RQGIS::run_qgis(alg = "gdalogr:proximity",
                           INPUT = x,
                           OUTPUT = filename,
                           VALUES = values,
                           UNITS = units,
                           RTYPE = 6,
                           load_output = TRUE)

    out <- raster::projectRaster(out, crs = in_crs)

  } else {

  out <- RQGIS::run_qgis(alg = "gdalogr:proximity",
                         INPUT = x,
                         OUTPUT = filename,
                         VALUES = values,
                         UNITS = units,
                         load_output = TRUE)
  }

  return(out)

}