R/CDF.R

In jacpete/jpfxns: JP Functions

#' @title Spatial Cummulative Distribution Function
#'
#' @description
#' Calculates a spatial cummulative distribution function for measuring selection
#' of resources in a study area.
#'
#' @details
#' If \code{replicationMethod == "features"}, features must be of type SpatialPoints*.
#'
#' Add more details like so.
#'
#' @param relocations SpatialPoints* object.
#' @param features Spatial* object (see \code{details} for restrictions).
#' @param studyArea SpatialPolygons* object.
#' @param replicationMethod Do you want to use pseudoreplications of relocations
#' or features. The default is set to relocations.
#' @param iterations The number of random samples to draw from the pseudoreplicated data.
#' The default is set to 30.
#' @param breakValue The distance (m) that the lags will be set to. The defalut is set
#' to 50.
#' @param randomMultiplier Determines the size of the random sample that each
#' pseudoreplication will be drawn from. Default is set to 30.
#'
#' @return An object of class CDF with a list of 5 vectors and 2 attributes. More detail.
# will need to look into defining a class in more detail
# @examples need to find test data and add examples
#' @export
CDF <- function(relocations, features, studyArea, replicationMethod = c("relocations", "features"), iterations = 30, breakValue = 50, randomMultiplier = 30) {

  ###Input Error Handling
  if (!grepl("SpatialPoints", class(relocations))) {
    stop("\"relocations\" input not of a \'SpatialPoints\' class")
  } else if (!grepl("Spatial", class(features))) {
    stop("\"features\" input not of a \'Spatial\' class")
  } else if (!grepl("SpatialPolygons", class(studyArea))) {
    stop("\"studyArea\" input not of a \'SpatialPolygons\' class")
  } else if (!(class(iterations) %in% c("numeric","integer"))) {
    stop("\"iterations\" input not a number")
  } else if (!(class(breakValue) %in% c("numeric","integer"))) {
    stop("\"breakValue\" input not a number")
  } else if (!(class(randomMultiplier) %in% c("numeric","integer"))) {
    stop("\"randomMultiplier\" input not a number")
  }
  replicationMethod <- match.arg(replicationMethod)


  ###Ensure the "relocations" and "feature" inputs are clipped to the studyArea
  relocations <- jpfxns::splitClip(relocations, studyArea, sep = 1000)
  features <- jpfxns::splitClip(features, studyArea, sep = 1000)
  if (is.null(relocations) | is.null(features)) {
    warning("No relocations and/or features in the given study area.")
    return(NULL)
  }
  relocations <- jpfxns::SPtoSPDF(relocations)
  features <- jpfxns::SPtoSPDF(features)

  ###Set which input will get replicated using replicationMethod argument
  if (replicationMethod == "relocations") {
    pseudoInput <- relocations
  } else if (replicationMethod == "features") {
    if (!grepl("SpatialPoints", class(features))){
      stop("\"features\" input not of a \'SpatialPoints\' class and replicationMethod set to \'features\'.
           Cannot replicate features types other than SpatialPoints")
    }
    pseudoInput <- features
  }


  ##############################
  ##############################
  ##############################
  ##############################
  ######This section would be changed to accommodate a density aware random sampling method


  LargePseudoSample <- SPtoSPDF(sp::spsample(studyArea, (randomMultiplier * nrow(pseudoInput)), type = "random"))

  pseudoSamplelist <- lapply(1:iterations, function(x) {
    LargePseudoSample[sample(length(LargePseudoSample), nrow(pseudoInput)),]
  })


  ##############################
  ##############################
  ##############################
  ##############################

  actualNearFeature <- splitNearDist(relocations, features, sep = 1000, output = "nearest")

  if (replicationMethod == "relocations") {
    pseudoNearFeatureList <- lapply(1:iterations, function(x) {
      splitNearDist(pseudoSamplelist[[x]], features, sep = 1000, output = "nearest")
    })
  } else if (replicationMethod == "features") {
    pseudoNearFeatureList <- lapply(1:iterations, function(x) {
      splitNearDist(relocations, pseudoSamplelist[[x]], sep = 1000, output = "nearest")
    })
  }

  pseudoNearFeatureDF <- do.call(cbind.data.frame, pseudoNearFeatureList)
  names(pseudoNearFeatureDF) <- as.character(seq_along(pseudoNearFeatureList))

  ##############################
  ##############################
  ##############################
  ##############################

  breaks <- seq(0, (max(c(max(actualNearFeature), max(pseudoNearFeatureDF))) + breakValue), breakValue)

  actualFrequency_by_breaks <- sapply(seq_along(breaks), function(x) {
    if (x != length(breaks)) {
      out <- length(which(actualNearFeature > breaks[x] & actualNearFeature <= breaks[x+1]))
    } else {
      out <- length(which(actualNearFeature > breaks[x]))
    }
    return(out)
  })

  actualCumSumFrequency <- cumsum(actualFrequency_by_breaks)
  actualCDF <- sapply(actualCumSumFrequency, function(x) x/sum(actualFrequency_by_breaks))


  pseudoCDFList <- lapply(1:ncol(pseudoNearFeatureDF), function(y) {
    pseudoFrequency_by_breaks <- sapply(seq_along(breaks), function(x) {
      if (x != length(breaks)) {
        out <- length(which(pseudoNearFeatureDF[,y] > breaks[x] & pseudoNearFeatureDF[,y] <= breaks[x+1]))
      } else {
        out <- length(which(pseudoNearFeatureDF[,y] > breaks[x]))
      }
      return(out)
    })
    pseudoCumSumFrequency <- cumsum(pseudoFrequency_by_breaks)
    pseudoCDF <- sapply(pseudoCumSumFrequency, function(x) x/sum(pseudoFrequency_by_breaks))
  })

  pseudoCDF_DF <- do.call(cbind.data.frame, pseudoCDFList)
  names(pseudoCDF_DF) <- as.character(seq_along(pseudoCDFList))

  pseudoCDF_Mean <- rowMeans(pseudoCDF_DF)


  #Selection Neutral Avoidance Points from CDF's
  CDF_SNA <- actualCDF -  pseudoCDF_Mean


  StandardError <- sapply(seq(1, nrow(pseudoCDF_DF), 1), function(x) {
    sd(pseudoCDF_DF[x,]) / sqrt(ncol(pseudoCDF_DF))
  })
  #######################################
  #######################################
  #######################################
  #######################################

  ###Output
  output <- list(breaks, actualCDF, pseudoCDF_Mean, CDF_SNA, StandardError)
  names(output) <- c("Breaks", "ActualCDF", "PseudoCDF", "CDF_SNA", "SE")
  class(output) <- "CDF"
  attr(output, 'iterations') <- iterations

  return(output)
}