R/SimDist.R

In LewisAJones/LBGSim: Latitudinal biodiversity gradient simulator

#' Randomly simulates the distribution of a species based on probability raster and distance gradient raster.
#'
#' @param n Integer. Number of species to simulate distributions for. Defaults to 10.
#' @param tmplayer RasterLayer. Probability raster of study area. Typically, raster should be masked to suitable environmental conditions (e.g. <200 m depth)
#' @param circ_div_grad RasterLayer. Distance gradient for populating species distribution.
#' @param parallel Logical. Implement parallel processing. Defaults to FALSE. Parallel processing is currently only supported for "unix" operating systems.
#' @return A dataframe of species occurrences (x,y).  
#' @importFrom pbmcapply pbmclapply
#' @importFrom parallel detectCores
#' @importFrom virtualspecies sampleOccurrences
#' @importFrom raster extent mask resample res cellStats getValues
#'
simDist <- function(n = 10, tmplayer, circ_div_grad, parallel = FALSE){
  if(parallel == TRUE){
    data <- pbmcapply::pbmclapply(1:n, function(i) {
      MaxN <- raster::res(tmplayer)[1] #generate number of maximum points based on resolution of LBG raster
      MaxN <- 150/MaxN #generate number of maximum points based on resolution of LBG raster
      geoxx <- seq(1,180,1)
      geoxx <- dgeom(geoxx-1, prob = 1/24)
      geoyy <- seq(1,90,1)
      geoyy <- dgeom(geoyy-1, prob = 1/24)
      geord <- seq(1,MaxN,1)
      geord <- dgeom(geord-1, prob = 1/10)
      mxx <- sample(1:180, 1, prob = geoxx) #maximum longitude range size in degree
      mxy <- sample(1:90, 1, prob = geoyy) #maximum latitude range size in degree
      tmprp <- virtualspecies::sampleOccurrences(x = tmplayer, n = 1, bias = "manual", weights = tmplayer, replacement = FALSE, plot = FALSE)
      tmprp <- as.matrix(tmprp$sample.points[,1:2])
      xmax <- tmprp[1]+mxx #set up extent based on initial random point
      xmin <- tmprp[1]-mxx #set up extent based on initial random point
      ymax <- tmprp[2]+mxy #set up extent based on initial random point
      ymin <- tmprp[2]-mxy #set up extent based on initial random point
      ex <- raster::extent(xmin, xmax, ymin, ymax) #set up extent based on initial random point
      tmp_grad <- circ_div_grad #set up distance gradient for random point generation
      raster::extent(tmp_grad) <- ex #set up distance gradient for random point generation
      tmp_grad <- raster::resample(tmp_grad, tmplayer, method = "ngb", na.rm = TRUE) #alter distance gradient raster to same extent and resolution
      tmp_grad <- raster::mask(tmp_grad, tmplayer) #mask probability raster to shallow water area
      tmp_grad <- (tmp_grad/raster::cellStats(tmp_grad, 'max'))
      rdv <- sample(1:MaxN, 1, prob = geord) #generate number of random occurrences
      if(length(unique(raster::getValues(tmp_grad))) < 3){
        x <- NA
        y <- NA
        rp <- cbind.data.frame(x,y)  }
      else if(sum(!is.na(raster::getValues(tmp_grad))) <= 2){
        x <- NA
        y <- NA
        rp <- cbind.data.frame(x,y)
      }
      else{rp <- virtualspecies::sampleOccurrences(x = tmp_grad, n = rdv, bias = "manual", weights = tmp_grad, replacement = TRUE, plot = FALSE)
      rp <- as.data.frame(rp$sample.points[,1:2])}
      rp}, mc.cores = parallel::detectCores(), mc.preschedule = TRUE, mc.cleanup = TRUE)
  } else {
    data <- lapply(1:n, function(i){
    MaxN <- raster::res(tmplayer)[1] #generate number of maximum points based on resolution of LBG raster
    MaxN <- 150/MaxN #generate number of maximum points based on resolution of LBG raster
    geoxx <- seq(1,180,1)
    geoxx <- dgeom(geoxx-1, prob = 1/24)
    geoyy <- seq(1,90,1)
    geoyy <- dgeom(geoyy-1, prob = 1/24)
    geord <- seq(1,MaxN,1)
    geord <- stats::dgeom(geord-1, prob = 1/10)
    mxx <- sample(1:180, 1, prob = geoxx) #maximum longitude range size in degree
    mxy <- sample(1:90, 1, prob = geoyy) #maximum latitude range size in degree
    #tmprp <- dismo::randomPoints(mask = tmplayer, n = 1, prob = TRUE, lonlatCorrection = TRUE) #create random point for initial source population
    tmprp <- virtualspecies::sampleOccurrences(x = tmplayer, n = 1, bias = "manual", weights = tmplayer, replacement = FALSE, plot = FALSE)
    tmprp <- as.matrix(tmprp$sample.points[,1:2])
    xmax <- tmprp[1]+mxx #set up extent based on initial random point
    xmin <- tmprp[1]-mxx #set up extent based on initial random point
    ymax <- tmprp[2]+mxy #set up extent based on initial random point
    ymin <- tmprp[2]-mxy #set up extent based on initial random point
    ex <- raster::extent(xmin, xmax, ymin, ymax) #set up extent based on initial random point
    tmp_grad <- circ_div_grad #set up distance gradient for random point generation
    raster::extent(tmp_grad) <- ex #set up distance gradient for random point generation
    tmp_grad <- raster::resample(tmp_grad, tmplayer, method = "ngb", na.rm = TRUE) #alter distance gradient raster to same extent and resolution
    tmp_grad <- raster::mask(tmp_grad, tmplayer) #mask probability raster to shallow water area
    tmp_grad <- (tmp_grad/raster::cellStats(tmp_grad, 'max'))
    rdv <- sample(1:MaxN, 1, prob = geord) #generate number of random occurrences
    if(length(unique(raster::getValues(tmp_grad))) < 3){
      x <- NA
      y <- NA
      rp <- cbind.data.frame(x,y)  }
    else if(sum(!is.na(getValues(tmp_grad))) <= 2){
      x <- NA
      y <- NA
      rp <- cbind.data.frame(x,y)
    }
    else{rp <- virtualspecies::sampleOccurrences(x = tmp_grad, n = rdv, bias = "manual", weights = tmp_grad, replacement = TRUE, plot = FALSE)
    rp <- as.data.frame(rp$sample.points[,1:2])}
    rp})}
  return(data)
}