R/LatSSC.R

In LewisAJones/LBGSim: Latitudinal biodiversity gradient simulator

#' Get spatial sampling coverage
#'
#' @param data Dataframe (colnames:x,y). A dataframe of the sampling window/occurrences (x,y).
#' @param id Vector. A vector of sampling IDs (e.g. collection numbers)
#' @param bandsize Integer. Size of latitudinal bins/bands desired for analyses. If set to 180, global spatial sampling coverage will be computed.
#' @param res Numeric. Value used to convert sampling points to raster. The value of res indicates the resolution of the raster generated. A higher value generates more generous sampling, while lower represents more conservative sampling. If used in conjunction with LBGSim(), res should equal the same as used by LBGtype().
#' @return A dataframe of the spatial sampling coverage within latitudinal bands
#' @importFrom sp SpatialPointsDataFrame
#' @importFrom raster rasterize raster rowSums mask
#' @export
LatSSC <- function(data, id, bandsize, mask = NULL, write = FALSE, name = "simulation"){
  nbins <- 180/bandsize
  bin <- as.data.frame(seq(1, nbins, 1))
  max <- as.data.frame(rev(seq((-90+bandsize), 90, bandsize))) 
  min <- as.data.frame(rev(seq(-90, (90-bandsize), bandsize)))
  mid <- as.data.frame(rev(seq((-90+(bandsize/2)), (90-(bandsize/2)), bandsize)))
  hemisphere <- c("N", "S")
  hemisphere <- as.data.frame(rep(hemisphere, each = nbins/2))
  if(nrow(hemisphere) == 0){
    hemisphere <- "Global"
  }
  latbin <- cbind(bin, max, min, mid, hemisphere)
  colnames(latbin) <- c("bin", "max", "min", "mid", "hemisphere")
  x <- data$x
  y <- data$y
  z <- id
  xyz <- cbind.data.frame(x, y, z)
  xyz <- sp::SpatialPointsDataFrame(coords = cbind.data.frame(xyz$x, xyz$y), data = xyz)
  r <- raster(res = 1)
  ras <- rasterize(xyz, r, 'z', function(x, ...) length(unique(na.omit(x))))
  if(length(mask) > 0){
    mask <- raster::resample(mask, ras)
    ras <- raster::mask(ras, mask)
    msk <- raster::rowSums(mask, na.rm = TRUE)
  }
  bin <- ras
  bin[bin >= 1] <- 1
  bin[bin != 1] <- NA
  rs <- raster::rowSums(bin, na.rm = TRUE)
  ssc <- vector()
  for(i in 1:(180/bandsize)){
    min <- (bandsize * i) - (bandsize -1)
    max <- bandsize * i
    sumRAS <- sum(rs[min:max])
    if(length(mask) > 0){
      sumMSK <- sum(msk[min:max])
      sumRAS <- (sumRAS/sumMSK)*100
    }
    else{sumRAS <- (sumRAS/(ncol(bin)*bandsize))*100}
    ssc[i] <- sumRAS
  }
  latbin <- cbind.data.frame(latbin, ssc)
  plot(ssc~mid, latbin, type = "o",  main = "Latitudinal spatial sampling coverage", ylab = "Spatial sampling coverage (%)", xlab = expression(Latitude~(degree)), pch = 20, xlim = c(-90, 90), ylim = c(0, max(latbin$ssc, na.rm=TRUE)))
  if(write == TRUE){
    suppressWarnings(dir.create("./spatial_coverage/"))
    png(paste("./spatial_coverage/", name, "_SSC_plot.png", sep = ""), width = 2400, height = 2000, res = 300)
    plot(ssc~mid, latbin, type = "o",  main = "Latitudinal spatial sampling coverage", ylab = "Spatial sampling coverage (%)", xlab = expression(Latitude~(degree)), pch = 20, xlim = c(-90, 90), ylim = c(0, max(latbin$ssc, na.rm=TRUE)))
    dev.off()
    write.csv(latbin, paste("./spatial_coverage/", name, "_SSC.csv", sep = ""))
  }
  return(latbin)
}