R/volumeUD.R

In aspillaga/fishtrack3d: Analysis of Passive Acoustic Telemetry Data in 3D

#' Calculate utilization distribution probability volumes
#'
#' This function calculates the utilization distribution probability volumes
#' from 2D or 3D UD values.
#'
#' @param ud a \code{RasterLayer} (2D), \code{RasterStack} or
#'     \code{RasterBrick} (3D) object with UD values.
#' @param ind.layer logical. If \code{TRUE}, the UD volume is calculated for
#'     each layer separately (each layer in the raster belongs to a different
#'     individual or a different time-period). If \code{FALSE} (the default),
#'     UD volume is calculated taking into account all the layers (for UD-3D,
#'     where all the layers correspond to different depth-intervals for the
#'     same individual and time-period).
#'
#' @return a \code{RasterLayer} or a \code{RasterStack} object with UD
#'     probability volumes.
#'
#' @export
#'
#'
volumeUD <- function(ud, ind.layer = FALSE) {

  if (ind.layer & nlayers(ud) > 1) {
    return(stack(lapply(unstack(ud), volumeUD)))
  }

  # Check if arguments are correct =============================================
  if (is.null(ud) | !class(ud) %in% c("RasterLayer", "RasterStack",
                                      "RasterBrick")) {
    stop(paste("Utilization distributions ('ud') must be in a 'RasterLayer',",
               "'RasterStack' or 'RasterBrick' object."), call. = FALSE)
  }

  if (round(sum(values(ud), na.rm = TRUE), 7) != 1) {
    stop("All the UDs must sum 1.")
  }

  names <- names(ud)
  rank <- (1:length(raster::values(ud)))[rank(raster::values(ud))]
  raster::values(ud) <- 1 - cumsum(sort(raster::values(ud)))[rank]
  names(ud) <- names

  return(ud)

}



#' Predict utilization distribution (UD) values from a \code{kde} object
#'
#' This function takes a \code{kde} object and predicts the utilization
#' distribution (UD) values for a 3D grid defined by a raster (horizontal
#' resolution) and a sequence of depths (vertical resolution.)
#'
#' @param kde a kde object.
#' @param raster raster to extract the horizontal coordinates of the 3D grid.
#' @param depths vector with depths values defining the vertical resolution of
#'     the 3D grid.
#'
#' @return A \code{RasterStack} object with the UD volumes of a different depth
#'     interval in each layer.
#'
#' @import raster
#' @import sp
#'
#' @export
#'
#'
predictKde <- function(kde, raster, depths) {

  # Check if arguments are correct =============================================
  if (is.null(kde) | class(kde) != "kde") {
    stop("The 'kde' object must be a 'kde' object from the 'ks' package.",
         call. = FALSE)
  }

  if (is.null(raster) | class(raster) != "RasterLayer") {
    stop("The 'raster' object must be a 'RasterLayer' object.", call. = FALSE)
  }

  pred <- lapply(depths, function(d) {
    rast.t <- raster::raster(raster)
    raster::values(rast.t) <- predict(kde,
                                      x = data.frame(sp::coordinates(rast.t),
                                                     z = -d))
    return(rast.t)
  })

  pred <- raster::stack(pred)
  pred <- pred / sum(raster::values(pred))
  names(pred) <- paste0("d", depths)

  return(pred)

}



#' Compute 3D contour mesh
#'
#' This function takes a a \code{RasterStack} object with UD volumes
#' (generated by the \code{predictKde} function) and generates a 3D mesh to
#' be ploted.
#'
#' @param rast.vol a \code{RasterStack} or \code{RasterBrick} object with
#'     UD contour volumes.
#' @param levels vector with the probabilities of the contours to generate.
#'
#' @return A \code{list} with the 3D contour meshes.
#'
#' @export
#'
#'
ud3dmesh <- function(rast.vol, levels = c(0.5, 0.95)) {

  # Check if arguments are correct =============================================
  if (is.null(rast.vol) | !class(rast.vol) %in% c("RasterBrick",
                                                  "RasterStack")) {
    stop(paste("The 'raster.vol' object must be a 'RasterBrick' or",
               "'RasterStack' object."), call. = FALSE)
  }

  vol <- plyr::laply(unstack(rast.vol), as.matrix)

  m <- list(x = sort(unique(coordinates(rast.vol)[, "x"])),
            y = rev(sort(unique(coordinates(rast.vol)[, "y"]))),
            z = as.numeric(substr(names(rast.vol), 2, 6)),
            v = aperm(vol, c(3:1)))

  contours <- lapply(levels, function(l) {
    contour <- misc3d::computeContour3d(m$v, level = l, x = m$x, y = m$y, z = m$z)
    indx <- matrix(0:(nrow(contour) - 1), ncol = 3, byrow = TRUE)
    return(list(contour = contour, indx = indx))
  })

  contours <- as.list(contours)
  names(contours) <- levels

  return(contours)

}