R/LBGtype.R

In LewisAJones/LBGSim: Latitudinal biodiversity gradient simulator

#' Create probability grid matching desired latitudinal biodiversity gradient.
#' 
#' Generates a modern-, bimodal, or flat-type latitudinal biodiversity gradient. Flat-type also includes a 'null model', with a completely flat gradient.
#'
#' @param degree Integer. Steepness (values = 1:5) of latitudinal biodiveristy gradient with 1 representing the shallowest gradient and 5 the steepest. If flat-type is selected in 'type', steepness can be set to 'NULL' for a null gradient.
#' @param res Numeric. Resolution of the raster output. This should match the desired resolution of the study. Higher resolution raster allows the generation of more accurate species distributions but also increases computation time.
#' @param type Character. Define type of latitudinal biodiversity gradient desired (options: "modern", "flat", "bimodal")
#' @param restrict Logical. define whether you wish the probability gradient to be restricted. This prevents the generation of occurrences in low probability cells. Defaults to FALSE.
#' @param write Logical. If TRUE, function writes ascii grid of probability gradient and png to working directory. Defaults to FALSE.
#' @return RasterLayer. Desired probability gradient. 
#' @importFrom raster raster resample scale mosaic
#' @importFrom rasterVis levelplot
#' @importFrom NLMR nlm_edgegradient
#' @export
LBGtype <- function(degree, res, type, restrict = FALSE, write = FALSE){
  if(is.null(degree) == TRUE){
    div_grad <- raster::raster(res = res, val = 1)
    suppressWarnings(dir.create(paste("./Simulation/", sep = "")))
    png(paste("./Simulation/LBG-plot.png", sep = ""), width = 2400, height = 2000, res = 300)
    suppressWarnings(print(rasterVis::levelplot(div_grad, margin=T, scales.margin=list(x=c(0,0), y=NULL), main = NA, par.settings = rasterVis::viridisTheme())))
    dev.off()
    writeRaster(div_grad, paste("./Simulation/LBG-probability-grid.asc", sep = ""), format = "ascii", overwrite = TRUE)
    suppressWarnings(print(rasterVis::levelplot(div_grad, margin=T, scales.margin=list(x=c(0,0), y=NULL), main = NA, par.settings = rasterVis::viridisTheme())))
    return(div_grad)
  }
  if(degree > 5 | degree < 1){
    stop("Degree value not accepted. Please use a value between 1 and 5")
  }
  
  r <- raster::raster(res = res, val = 1) #resolution needs to be higher if you wish to produce many geographically unique occurrences
  
  #make LBG probability raster
  div_grad <- NLMR::nlm_edgegradient(ncol = 200, #higher resolution for more occurrences
                               nrow = 200,
                               direction = 180)
  if(type == "modern"){
    tmpdeg <- round(degree)*2
    #make LBG probability raster
    
    e <- extent(-180, 180, -90, 90) #create extent for raster 
    extent(div_grad) <- e #merge to right extent for global data 
    div_grad <- raster::resample(div_grad, r) #resample to appropiate resolution
    div_grad <- div_grad^tmpdeg #increase power to increase gradient if desired, 9 gives a reasonable modern-type LBG
    div_grad <- raster::scale(div_grad, center = FALSE, scale = TRUE)
  }
  
  else if(type == "flat"){
    tmpdeg <- round(degree)/5
    #make LBG probability raster
    
    e <- extent(-180, 180, -90, 90) #create extent for raster 
    extent(div_grad) <- e #merge to right extent for global data 
    div_grad <- raster::resample(div_grad, r) #resample to appropiate resolution
    div_grad <- div_grad^tmpdeg #increase power to narrow gradient if desired, 9 gives a reasonable modern-type LBG
    div_grad <- raster::scale(div_grad, center = FALSE, scale = TRUE)
  }
  
  
  else if(type == "bimodal"){
    tmpdeg <- round(degree)*2
    a <- div_grad
    extent(a) <- extent(-180, 0, -90, 0)
    b <- div_grad
    extent(b) <- extent(0, 180, -90, 0)
    c <- div_grad
    extent(c) <- extent(0, 180, 0, 90)
    d <- div_grad
    extent(d) <- extent(-180, 0, 0, 90)
    
    div_grad <- raster::mosaic(a,b,c,d, fun = mean)
    
    e <- extent(-180, 180, -90, 90) #create extent for raster 
    extent(div_grad) <- e #merge to right extent for global data 
    div_grad <- raster::resample(div_grad, r) #resample to appropiate resolution
    div_grad <- div_grad^tmpdeg #increase power to narrow gradient if desired, 15 gives a reasonable modern-type LBG #degree
    div_grad <- raster::scale(div_grad, center = FALSE, scale = TRUE)
  }
  
  else {
    stop("Latitudinal biodiversity gradient not recognised, use: 'modern', 'flat', or 'bimodal'.")
  }
  div_grad <- div_grad + 1
  div_grad <- (div_grad/cellStats(div_grad, 'max'))
  if(restrict == TRUE){
    mv <- minValue(div_grad)
    mv <- mv * 100
    mv <- ceiling(mv)
    mv <- mv/100
    div_grad[div_grad <= mv] <- NA
  }
  if(write == TRUE){
    suppressWarnings(dir.create(paste("./Simulation/", sep = "")))
    png(paste("./Simulation/LBG-plot.png", sep = ""), width = 2400, height = 2000, res = 300)
    suppressWarnings(print(rasterVis::levelplot(div_grad, margin=T, scales.margin=list(x=c(0,0), y=NULL), main = NA, par.settings = rasterVis::viridisTheme())))
    dev.off()
    writeRaster(div_grad, paste("./Simulation/LBG-probability-grid.asc", sep = ""), format = "ascii", overwrite = TRUE)
  }
  suppressWarnings(print(rasterVis::levelplot(div_grad, margin=T, scales.margin=list(x=c(0,0), y=NULL), main = NA, par.settings = rasterVis::viridisTheme())))
  div_grad
}