R/ersEx.R

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

#' @title Environmental representativeness score estimation (Ex-situ conservation indicators).
#' @name ers_exsitu
#' @description Performs an estimation of enviromental representativeness score for ex-situ gap analysis (ERSex) using Ramirez-Villegas et al., (2010) methodology,
#' this function uses a csv with coordinates, a germplasm buffer raster file and a vectorial file of eccoregions
#'  \deqn{ERSex = min(100,(Number of ecoregions where germplasm accessions are available/
#' Number of ecoregions where species is available)*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 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 the species was located \cr
#'  G_N_ECO \tab Number of ecosystems where germplasm accessions were available \cr
#'  ERS \tab ERSex result \cr
#' }
#'
#'
#'
#'
#' @examples ers_exsitu('Cucurbita_digitata',Workspace,'v1')
#'
#' 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_exsitu(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

##########################################   Start Functions    ###############################################
# This function calculates the ERSex It loads occurrences if they exist, then
# loads the presence/absence surface, creates the G buffer (i.e. CA50) and finally
# outputs the ERS and # eco classes in a data.frame (which is written into a file).
# @param (string) species: species ID
# @param (logical) debug: whether to save or not the intermediate raster outputs
# @return (data.frame): This function returns a data frame with ERS, # eco classes
#                       of G buffer (i.e. CA50) and of the presence/absence surface.

ers_exsitu <- function(species_list,occurrenceData, raster_list, bufferDistance) {

#load packages
suppressMessages(require(sp))
suppressMessages(require(raster))
suppressMessages(require(dplyr))
suppressMessages(require(tidyr))

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

  # load in ecoregion dataset
  ecoReg <- raster::shapefile(system.file("data/ecoRegion/tnc_terr_ecoregions.shp", package = "gapAnalysisR"))
  # maybe this directly downloads an element from the dataverse

  # generate a dataframe to store the output values
  df <- data.frame(matrix(ncol = 2, nrow = length(species_list)))
  colnames(df) <- c("species", "ERSex")

  # loop through all species
  for(i in 1:length(species_list)){
    speciesOcc <- occurrenceData %>%
      dplyr::filter(taxon == species_list[i])
    if(length(speciesOcc$type == "G") == 0){
      df$species[i] <- species_list[i]
      df$ERSex[i] <- 0
      }else{
        occDataG <- speciesOcc  %>%
          dplyr::filter(type == "G")%>%
          dplyr::select(longitude,latitude)
          sp::coordinates(occDataG) <- ~longitude+latitude
          sp::proj4string(occDataG) <- CRS("+proj=longlat +datum=WGS84")
        # select raster with species name
          for(j in 1:length(raster_list)){
            if(grepl(j, i, ignore.case = TRUE)){
              sdm <- raster_list[[j]]
            }
          }
        # convert SDM from binary to 1-NA for mask and area
        sdmMask <- sdm
        sdmMask[sdmMask == 0] <- NA

        # buffer G points
        buffer <- geobuffer::geobuffer_pts(xy = occDataG,
                                             dist_m = bufferDistance,
                                             output = 'sf')
        # rasterizing and making it into a mask
        buffer_rs <- fasterize::fasterize(buffer, sdm)
        buffer_rs[!is.na(buffer_rs[])] <- 1
        buffer_rs <- buffer_rs * sdmMask

        gPoints <- sp::SpatialPoints(raster::rasterToPoints(buffer_rs))
        # extract values from ecoregions to points
        raster::crs(gPoints) <- raster::crs(ecoReg)
        ecoValsG <- sp::over(x = gPoints, y = ecoReg) %>%
          dplyr::distinct(ECO_NUM)%>%
          tidyr::drop_na(ECO_NUM) %>% #ECO_ID
          dplyr::filter(ECO_NUM > 0) #ECO_ID

        # create point object from all the occurence data for the species
        sp::coordinates(speciesOcc) <- ~longitude+latitude
        sp::proj4string(speciesOcc) <- CRS("+proj=longlat +datum=WGS84")
        raster::crs(speciesOcc) <- raster::crs(ecoReg)

        # number of ecoregions present in model
        ecoVal <- data.frame(over(x = speciesOcc, y = ecoReg))%>%
            dplyr::select(ECO_NUM )%>% #ECO_ID
            dplyr::distinct() %>%
            tidyr::drop_na() %>%
            dplyr::filter(ECO_NUM > 0) # -9998 are lakes #ECO_ID != -9998

        #calculate ERS
        ers <- min(c(100, (nrow(ecoValsG)/nrow(ecoVal))*100))
        # assign values to df
        df$species[i] <- as.character(species_list[i])
        df$ERSex[i] <- ers
      }
  }
  return(df)

}