R/GRSex.R

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

#' @title Germplasm representativeness score estimation (Ex-situ conservation indicators).
#' @name grs_exsitu
#' @description Performs an estimation of germplasm representativeness score for ex-situ gap analysis (GRSex) using Ramirez-Villegas et al., (2010) methodology
#' This function uses a germplasm buffer raster file (e.g. CA50) and a thresholded species distribution model.
#'  \deqn{GRSex = min(100,(Germplasm buffer area/species distribution in 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 data frame file saved at gap_analysis folder with four columns:
#'
#' \tabular{lcc}{
#'  ID \tab Species name \cr
#'  SPP_AREA_km2 \tab Area occupied by the species using as input a SDM thresholded file in tiff format \cr
#'  G_AREA_km2 \tab Area occupied by the germplasm accessions in a species distribution model \cr
#'  GRS \tab GRSex result \cr
#' }
#'
#' @examples grs_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 <- grs_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 GRSex. It loads occurrences if they exist, then
# loads the presence/absence surface, creates the G buffer (i.e. CA50) and finally
# outputs the GRS and areas 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 GRS and areas of G buffer (i.e. CA50)
#                       and of the presence/absence surface.

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

suppressMessages(require(rgdal))
suppressMessages(require(raster))

#importFrom("methods", "as")
#importFrom("stats", "complete.cases", "filter", "median")
#importFrom("utils", "data", "memory.limit", "read.csv", "write.csv")
if(missing(bufferDistance)){
  bufferDistance <- 50000
}
# create a dataframe to hold the components
df <- data.frame(matrix(ncol = 2, nrow = length(species_list)))
colnames(df) <- c("species", "GRSex")

  for(i in 1:length(sort(species_list))){
  # select species G occurrences
  occData <- occurrenceData %>%
    dplyr::filter(taxon == species_list[i]) %>%
    dplyr::filter(type == "G")%>%
    dplyr::select(longitude,latitude)

    sp::coordinates(occData) <- ~longitude+latitude
    sp::proj4string(occData) <- 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 = occData,
                                       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
  # calculate area of buffer
  cell_size<-raster::area(buffer_rs, na.rm=TRUE, weights=FALSE)
  cell_size<-cell_size[!is.na(cell_size)]
  gBufferRas_area<-length(cell_size)*median(cell_size)

  # calculate area of the threshold model
  cell_size<- raster::area(sdmMask, na.rm=TRUE, weights=FALSE)
  cell_size<- cell_size[!is.na(cell_size)]
  pa_spp_area <<-length(cell_size)*median(cell_size)
  # calculate GRSex
  grs <- min(c(100, gBufferRas_area/pa_spp_area*100))

  df$species[i] <- as.character(species_list[i])
  df$GRSex[i] <- grs
    }

  return(df)
}