R/ersIn.R

In dcarver1/gapAnalysisR: An R package to calculate conservation indicators using species distribution models

#' @title Germplasm representativeness score estimation (In-situ conservation indicators).
#' @name ers_insitu
#' @description Performs an estimation of germplasm representativeness score for in-situ gap analysis (GRSin) using Khoury et al., (2019) methodology
#' This function uses a germplasm buffer raster file (e.g. CA50), a thresholded species distribution model, and a raster file of protected areas
#'  \deqn{ERSin = min(100,(Number of ecoregions where germplasm accessions in protected areas/
#' Number of ecoregions where species distribution is intersected with protected areas)*100)}
#'
#' @param species A name species compiled using '_'  to call the raster files (SDM and germplasm buffer)
#'  from Workspace/parameters/inputs folder and occurrences files from Workspace/parameter/occurrences folder
#' @param Workspace A forder where the pipeline will be executed
#' @param  run_version The version of the analysis used (e.g 'v1')
#'
#' @return It returns a raster file with the species distribution model restricted to the protected areas raster file provided.
#' Also, this function returns a data frame file saved at gap_analysis folder with four columns:
#'
#' \tabular{lcc}{
#'  ID \tab Species name \cr
#'  SPP_N_ECO \tab Number of ecosystems where species distribution model is within protected areas \cr
#'  SPP_WITHIN_PA_N_ECO \tab Number of ecosystems where germplasm accessions are within protected areas \cr
#'  ERS \tab ERSex result \cr
#' }
#'
#'
#' @examples ers_insitu('Cucurbita_digitata',Workspace,'v1')
#'  \dontrun
#'  Workspace  <-  'E:/CIAT/workspace/Workspace_test/workspace'
#'  run_version  <- 'v1'
#'  species_list <- c('Cucurbita_cordata',
#'   'Cucurbita_digitata',
#'   'Cucurbita_foetidissima',
#'   'Cucurbita_palmata')
#'
#'   run_version <-'v1'
#
#'  lapply(1:length(species_list),function(i){
#'     species <- species_list[[i]]
#'     x <- ers_insitu(species,Workspace,run_version)
#'     print(paste0(species,' DONE!'))
#'  })
#'
#'@references
#'
#'Ramírez-Villegas, J., Khoury, C., Jarvis, A., Debouck, D. G., & Guarino, L. (2010).
#'A Gap Analysis Methodology for Collecting Crop Genepools: A Case Study with Phaseolus Beans.
#'PLOS ONE, 5(10), e13497. Retrieved from https://doi.org/10.1371/journal.pone.0013497
#'
#' Khoury, C. K., Amariles, D., Soto, J. S., Diaz, M. V., Sotelo, S., Sosa, C. C., … Jarvis, A. (2019).
#' Comprehensiveness of conservation of useful wild plants: An operational indicator for biodiversity
#' and sustainable development targets. Ecological Indicators. https://doi.org/10.1016/j.ecolind.2018.11.016
#'
#' @export
ers_insitu <- function(species_list,occurrenceData,raster_list) {

  suppressMessages(require(sp))
  suppressMessages(require(raster))
  suppressMessages(require(dplyr))
  suppressMessages(require(tidyr))
  suppressMessages(require(rgdal))
  suppressMessages(require(tmap))
  suppressMessages(require(fasterize))
  suppressMessages(require(sf))

  #importFrom("methods", "as")
  #importFrom("stats", "complete.cases", "filter", "median")
  #importFrom("utils", "data", "memory.limit", "read.csv", "write.csv")

  ## ERSin analyzes how well protected areas cover the maxent model with regard to ecosystems covered.
  df <- data.frame(matrix(ncol=2, nrow = length(species_list)))
  colnames(df) <- c("species", "ERSin")
  # load in protect area raster
  proArea <- raster(system.file("data/protectedArea/wdpa_reclass.tif",
                                package = "gapAnalysisR"))
  # Load in ecoregions shp
  ecoReg <- raster::shapefile(system.file("data/ecoRegion/tnc_terr_ecoregions.shp", package = "gapAnalysisR"))


  for(i in 1:length(species_list)){
    # select threshold map for a given species
    for(j in 1:length(raster_list)){
      if(grepl(j, i, ignore.case = TRUE)){
        sdm <- raster_list[[j]]
      }
    }
    # select occurrence data for the given species
    occData1 <- occurrenceData %>%
      dplyr::filter(taxon == species_list[i])
    sp::coordinates(occData1) <- ~longitude+latitude
    raster::crs(occData1) <- raster::crs(ecoReg)
    # extract the ecoregion values to the points
    ecoVal <- data.frame(over(x = occData1, y = ecoReg))%>%
      dplyr::select(ECO_NUM )%>% #ECO_ID
      dplyr::distinct() %>%
      tidyr::drop_na() %>%
      dplyr::filter(ECO_NUM > 0) # -9998 are lakes #ECO_ID != -9998

    # number of ecoregions in modeling area
    ecoInSDM <- nrow(ecoVal)

    # mask protected areas to threshold
    proArea1 <- raster::crop(x = proArea, y=sdm)
    sdm[sdm == 0] <- NA
    proArea1 <- sdm * proArea1

    #convert protect area to points
    protectPoints <- sp::SpatialPoints(raster::rasterToPoints(proArea1))
    # extract values from ecoregions to points
    crs(protectPoints) <- crs(ecoReg)
    ecoValsPro <- sp::over(x = protectPoints, y = ecoReg) %>%
      dplyr::select(ECO_NUM )%>% #ECO_ID
      dplyr::distinct(ECO_NUM ) %>% #ECO_ID
      tidyr::drop_na()%>%
      dplyr::filter(ECO_NUM > 0) #ECO_ID
    # subset ecoVal in protect area so only values that are in the SDM area included.
    ecoInProt <- nrow(ecoValsPro)

    #clause for 9 in protected area
    if(ecoInProt == 0){
      df$species[i] <- as.character(species_list[i])
      df$ERSin[i] <- 0
    }else{
      # calculate ers
      ers <- min(c(100, (ecoInProt/ecoInSDM)*100))
      df$species[i] <- as.character(species_list[i])
      df$ERSin[i] <- ers
    }
  }
  return(df)
}