R/vulnerability.R

In rinnan/CENFA: Climate and Ecological Niche Factor Analysis

#' Climatic vulnerability
#'
#' Calculates the climatic vulnerability of a species using a \code{cnfa} and
#' \code{departure} object.
#'
#' @param cnfa Object of class \code{cnfa}
#' @param dep Object of class \code{departure}
#' @param method character. What type of mean should be used to combine sensitivity
#'   and exposure. Choices are "arithmetic" and "geometric"
#' @param w numeric. Optional vector of length two specifying the relative weights of
#'   sensitivity and exposure. See Details
#' @param parallel logical. If \code{TRUE} then multiple cores are utilized
#' @param n numeric. Number of cores to use for calculation
#' @param filename character. Output filename (optional)
#' @param ... Additional arguments for file writing as for \code{\link[raster]{writeRaster}}
#'
#' @details
#' The values of the vulnerability raster are calculated by combining the sensitivity
#' \eqn{\sigma} and the exposure \eqn{\epsilon}. If \code{method = "arithmetic"},
#' they will be combined as
#'
#'   \eqn{\nu = (w_1\sigma + w_2\epsilon) / (\sum_i w_i).}
#'
#' If \code{method = "geometric"}, they will be combined as
#'
#'   \eqn{\nu = \sqrt(\sigma * \epsilon).}
#'
#' @references
#' Rinnan, D. Scott and Lawler, Joshua. Climate-niche factor analysis: a spatial
#' approach to quantifying species vulnerability to climate change. Ecography (2019):
#' <doi:10.1111/ecog.03937>.
#'
#' @examples
#' \donttest{
#' mod1 <- cnfa(x = climdat.hist, s.dat = ABPR, field = "CODE")
#' dep <- departure(x = climdat.hist, y = climdat.fut, s.dat = ABPR)
#' vuln <- vulnerability(cnfa = mod1, dep = dep)
#' }
#'
#' @return Returns an S4 object of class \code{vulnerability} with the following slots:
#' \describe{
#'   \item{call}{Original function call}
#'   \item{vf}{Vulnerability factor. Vector of length p that describes the amount of
#'    vulnerability in each climate variable}
#'   \item{vulnerability}{Magnitude of the vulnerability factor}
#'   \item{ras}{RasterLayer of climate vulnerability}
#'   \item{weights}{Raster layer of weights used for departure calculation}
#' }
#'
#' @seealso \code{\link{departure}}
#'
#' @export

setGeneric("vulnerability", function(cnfa, dep, method = "geometric", w, parallel = FALSE, n = 1, filename = "", ...) {
  standardGeneric("vulnerability")
})

#' @rdname vulnerability
setMethod("vulnerability",
          signature(cnfa = "cnfa", dep = "departure"),
          function(cnfa, dep, method = "geometric", w, parallel = FALSE, n = 1, filename = "", ...) {

            call <- sys.call(sys.parent())
            call <- match.call(vulnerability, call)
            if (missing(w)) w <- c(1, 1)
            if (length(w) != 2) {
              warning("more than two weights supplied; defaulting to equal weights.")
              w <- c(1, 1)
            }

            d <- dep@df + 1
            s <- cnfa@sf
            ras <- dep@ras

            if (method == "arithmetic") {
              v <- (s*w[1] + d*w[2]) / sum(w)
            } else if (method == "geometric") {
              v <- (s^w[1] * d^w[2])^(1 / sum(w))
            }
            names(v) <- names(s)
            V <- sqrt(mean(v))

            filename <- trim(filename)
            if (!canProcessInMemory(ras) && filename == '') {
              filename <- rasterTmpFile()
            }

            s.map <- sensitivity_map(cnfa, parallel = parallel, n = n)
            e.map <- exposure_map(dep, parallel = parallel, n = n)
            #e.map <- e.map + 1
            if (method == "arithmetic") {
              f1 <- function(x,y) (x*w[1] + y*w[2]) / sum(w)
              vuln.ras <- overlay(s.map, e.map, fun = f1, filename = filename, ...)
            } else if (method == "geometric") {
              if(w[1] == w[2]) {
                w <- c(1, 1)
              } else {
                w <- w / sum(w)
              }
              f1 <- function(x,y) (x^w[1] * y^w[2])^(1 / sum(w))
              vuln.ras <- overlay(s.map, e.map, fun = f1, filename = filename, ...)

            }

            vuln <- methods::new("vulnerability", call = call, vf = v, vulnerability = V, ras = vuln.ras, weights = dep@weights)
            return(vuln)
          }
)