R/ERSin.R

In GapAnalysis: Conservation Indicators Using Spatial Information

#' @title Ecological representativeness score in situ
#' @name ERSin
#' @description The ERSin process provides an ecological measurement of the proportion of a species range
#'  that can be considered to be conserved in protected areas. The ERSin calculates the proportion of ecoregions
#'  encompassed within the range of the taxon located inside protected areas to the ecoregions encompassed
#'  within the total area of the distribution model, considering comprehensive conservation to have been accomplished
#'  only when every ecoregion potentially inhabited by a species is included within the distribution of the species
#'  located within a protected area.
#'
#' @param taxon A character object that defines the name of the species as listed in the occurrence dataset
#' @param sdm a terra rast object that represented the expected distribution of the species
#' @param occurrenceData a data frame of values containing columns for the taxon, latitude, longitude, and type. Coordinates are assumed to be in the WGS84 (EPSG:4326) coordinate reference system.
#' @param protectedAreas A terra rast object the contian spatial location of protected areas.
#' @param ecoregions A terra vect object the contains spatial information on all ecoregions of interests
#' @param idColumn A character vector that notes what column within the ecoregions object should be used as a unique ID
#' @param limitByPoints A boolean parameter (TRUE/FALSE) to determine if you want to limit the ecoregions considered to those with observations present.
#' TRUE will exclude all ecoregions with no points within. FALSE will include all ecoregions.
#' This was implemented to prevent edge effects where pixels from the distribution extend into
#' neighboring ecoregions as a product of differences in raster/vector geometry rather than being predicted there directly.
#'
#' @return A list object containing
#' 1. results : a data frames of values summarizing the results of the function
#' 2. missingEcos : a terra vect object showing all the ecoregions within the distribution with no protected areas present
#' 3. map : a leaflet object showing the spatial results of the function
#'
#'
#' @examples
#' ##Obtaining occurrences from example
#' data(CucurbitaData)
#' ##Obtaining Raster_list
#' data(CucurbitaRasts)
#' ##Obtaining protected areas raster
#' data(ProtectedAreas)
#' ## ecoregion features
#' data(ecoregions)
#'
#' # convert the dataset for function
#' taxon <- "Cucurbita_cordata"
#' sdm <- terra::unwrap(CucurbitaRasts)$cordata
#' protectedAreas <- terra::unwrap(ProtectedAreas)
#' ecoregions <- terra::vect(ecoregions)
#'
#' #Running ERSin
#' ers_insitu <- ERSin(taxon = taxon,
#'                     sdm = sdm,
#'                     occurrenceData = CucurbitaData,
#'                     protectedAreas = protectedAreas,
#'                     ecoregions = ecoregions,
#'                     idColumn = "ECO_NAME",
#'                     limitByPoints = FALSE
#'                     )
#'
#'
#' @references
#' Khoury et al. (2019) Ecological Indicators 98:420-429. \doi{10.1016/j.ecolind.2018.11.016}
#' Carver et al. (2021) GapAnalysis: an R package to calculate conservation indicators using spatial information
#' @importFrom dplyr tibble pull
#' @importFrom terra crop aggregate zonal
#' @importFrom leaflet addTiles addPolygons addLegend addRasterImage addCircleMarkers
#' @export

ERSin <- function(
  taxon,
  sdm,
  occurrenceData,
  protectedAreas,
  ecoregions,
  idColumn,
  limitByPoints = FALSE
) {
  # filter the occurrence data to the species of interest
  d1 <- occurrenceData |>
    dplyr::filter(occurrenceData$species == taxon) |>
    terra::vect(
      geom = c("longitude", "latitude"),
      crs = "+proj=longlat +datum=WGS84"
    )
  # add color
  d1$color <- ifelse(d1$type == "H", yes = "#1184d4", no = "#6300f0")
  # limit ecoregions to point locations
  if (isTRUE(limitByPoints)) {
    ecoregions <- ecoregions[d1, ]
  }
  # set id column for easier indexing
  ecoregions$id_column <- as.data.frame(ecoregions)[[idColumn]]

  # aggregate spatial features
  # guard: limitByPoints can leave no ecoregions when the taxon has no
  # usable coordinates, and terra::aggregate errors on an empty SpatVector
  if (nrow(ecoregions) > 0) {
    ecoregions <- terra::aggregate(x = ecoregions, by = "id_column")
  }

  # detect no-model case
  noModel <- !inherits(sdm, "SpatRaster") || terra::nlyr(sdm) == 0

  if (noModel) {
    nEcoModel <- 0
    nProModel <- 0
    ers <- 0
    selectedEcos <- ecoregions[0, ]
    protectedEcos <- ecoregions[0, ]
    missingEcos <- ecoregions[0, ]
    proMask <- NULL
  } else {
    # crop protected areas to sdm
    pro <- terra::crop(protectedAreas, sdm)

    # mask to model
    proMask <- pro * sdm

    # crop ecos to sdm
    eco <- terra::crop(ecoregions, sdm)

    # get ecoregions in sdm
    eco$totEco <- terra::zonal(
      x = sdm,
      z = eco,
      fun = "sum",
      na.rm = TRUE
    ) |>
      dplyr::pull()

    selectedEcos <- eco[eco$totEco > 0, ]
    nEcoModel <- nrow(selectedEcos)

    # get ecoregions in protected areas
    eco$totPro <- terra::zonal(
      x = proMask,
      z = eco,
      fun = "sum",
      na.rm = TRUE
    ) |>
      dplyr::pull()

    protectedEcos <- eco[eco$totPro > 0, ]
    nProModel <- nrow(protectedEcos)

    # get missing ecos
    missingEcos <- selectedEcos[!selectedEcos$id_column %in% protectedEcos$id_column, ]

    # calculate ERS
    if (nProModel == 0) {
      ers <- 0
    } else {
      ers <- (nProModel / nEcoModel) * 100
    }
  }

  # results
  df <- dplyr::tibble(
    Taxon = taxon,
    "Ecoregions within model" = nEcoModel,
    "Ecoregions with protected areas" = nProModel,
    "ERS insitu" = ers
  )

  # generate the base map
  map_title <- "<h3 style='text-align:center; background-color:rgba(255,255,255,0.7); padding:2px;'>Ecoregions within the SDM without Protected Area</h3>"

  map <- leaflet::leaflet() |>
    leaflet::addTiles()

  if (nrow(selectedEcos) > 0) {
    map <- map |>
      leaflet::addPolygons(
        data = selectedEcos,
        color = "#444444",
        weight = 1,
        opacity = 1.0,
        popup = ~id_column,
        fillOpacity = 0.5,
        fillColor = "#44444420"
      )
  }

  if (nrow(missingEcos) > 0) {
    map <- map |>
      leaflet::addPolygons(
        data = missingEcos,
        color = "#444444",
        weight = 1,
        opacity = 1.0,
        popup = ~id_column,
        fillOpacity = 0.5,
        fillColor = "#f0a01f"
      )
  }

  if (!noModel) {
    map <- map |>
      leaflet::addRasterImage(
        x = sdm,
        colors = "#47ae24"
      ) |>
      leaflet::addRasterImage(
        x = proMask,
        colors = "#746fae"
      )
  }

  map <- map |>
    leaflet::addLegend(
      position = "topright",
      title = "ERS in situ",
      colors = c("#47ae24", "#746fae", "#f0a01f", "#44444440"),
      labels = c("Distribution", "Protected Areas", "Eco gaps", "All Ecos"),
      opacity = 1
    ) |>
    leaflet::addControl(html = map_title, position = "bottomleft")

  # output
  output <- list(
    results = df,
    missingEcos = missingEcos,
    map = map
  )

  return(output)
}