R/extractPAstatistics.R

In cjcampbell/SDMetrics: Metrics To Compare Temporally Explicit Species Distribution Models

#' Extract comparative statistics from two presence/absence RasterLayers.
#'
#'
#' Extract metrics of range size change, seasonality, and centroid distance
#' from two presence/absence surfaces. Such surfaces are often generated by
#' thresholding a continuous probability-of-origin surface for conversion to a
#' binary surface.
#'
#' RasterLayers must have identical resolutions and extents.
#'
#' @param PA1 Presence/absence RasterLayer from season 1
#' @param PA2 Presence/absence RasterLayer from season 2
#' @param ... Additional arguments to append to the returned data.frame
#' @param appendPA True/false to paste character string "PA_" to the column name of the output metrics.
#'
#' @importFrom raster compareRaster
#' @importFrom raster reclassify
#' @importFrom stars st_as_stars
#' @importFrom sf st_as_sf
#' @importFrom sf st_combine
#' @importFrom sf st_make_valid
#' @importFrom sf st_intersection
#' @importFrom sf st_difference
#' @importFrom sf st_union
#'
#' @return A data.frame containing presence-absence surface based metrics
#'
#' @export
extractPAstatistics <- function(PA1, PA2, ..., appendPA = TRUE) {
  stopifnot(
    {"Class of PA1 and PA2 must be the same" = class(PA1) == class(PA2) },
    {"Class of PA arguments must be RasterLayer" = class(PA1) == c("RasterLayer")},
    {"PA1 and PA2 must have identical extent and resolutions" =
      raster::compareRaster(PA1, PA2) }
    )

  PA_sf1 <- PA1 %>%
    raster::reclassify( matrix(c(-Inf,0.5,NA), ncol = 3, byrow = T)) %>%
    stars::st_as_stars() %>%
    sf::st_as_sf(merge = T) %>%
    sf::st_combine() %>%
    sf::st_make_valid()
  PA_sf2 <- PA2 %>%
    raster::reclassify( matrix(c(-Inf,0.5,NA), ncol = 3, byrow = T)) %>%
    stars::st_as_stars() %>%
    sf::st_as_sf(merge = T) %>%
    sf::st_combine() %>%
    sf::st_make_valid()

  df <- statsFromPolygons(PA_sf1, PA_sf2, ...)

  if(appendPA) {
    df <- dplyr::rename_with(df, .fn = ~paste0("PA_", .), cols = vars(c(
      "rangeSizeChange", "yrRoundVseasonal", "centroidDist_km" )))
  }

  # data.frame with returned details.
  return(df)
}