R/makeSims.R

In cjcarlson/spatExtinct: A Package for Spatial Interpolation of Extinction Dates from Sighting Data

#' Simulations, without error
#' 
# "The tiger in my tank is going to go extinct, and I'm not feeling so good myself." --Eels, "Tiger in my Tank"
#'  
#' @description
#' Returns a list of two objects: sightings (data frame of sightings with columns "year", "longitude", "latitude", "sighting", and "real") and grid (raster onto which to interpolate your results)
#'
#' @param n Dimension of simulated landscape (nxn grid)
#' @param pts Number of valid sightings per year
#'
#' @keywords simulations
#'
#' @export
#'
#' @examples
#'


makeSims <- function(n=10, pts=1) {
  options(warn=-1)
  r <- raster(matrix(nrow=n,ncol=n),xmn=0,xmx=n,ymn=0,ymx=n)
  for (i in 1:n) {
    r[i,] <- n+1-i
  }

  n <- n*5
  r <- r*5

  d <- data.frame(longitude=NA,latitude=NA,year=NA)

  for (i in 1:(n-1)) {
    tmp <- r>i
    tmp[tmp==0] <- NA
    x <- data.frame(randomPoints(tmp,n=pts))
    x$year <- i
    names(x) <- c('longitude','latitude','year')
    d <- rbind(d,x)
  }

  d <- na.omit(d)
  dx <- diff(c(xmin(r), xmax(r))) / ncol(r) / 2 # Half of horizontal width
  dy <- diff(c(ymin(r), ymax(r))) / nrow(r) / 2 # Half of vertical width
  xy <- coordinates(d[,1:2])                             # 2-column matrix of coordinates
  n <- nrow(xy)                                       # Number of sample points
  xy <- xy + c(runif(n, -dx, dx), runif(n, -dy, dy))  # Add random changes
  d[,1:2] <- xy
  d$sighting <- 1

  grid <- r*0
  options(warn=0)
  final <- list(grid=stack(r,grid),sightings=d)
  names(final$grid) <- c('real','blank')
  return(final)
}

#'Simulations, with error
#'
# "Do you know where you go? You're headed on the strings, for the E-X-T-I-N-C-T." --alt-J, "The Gospel of John Hurt"
#'  
#'@description 
#' Produces simulations of extinction dates and sightings on a known landscape with mixed-certainty, mixed-accuracy data
#'
#' Returns a list of two objects: sightings (data frame of sightings with columns "year", "longitude", "latitude", "sighting", and "real") and grid (raster onto which to interpolate your results)
#'
#'
#'
#' @param n Dimension of simulated landscape (nxn grid)
#' @param pts Number of valid sightings per year
#' @param ipts Number of invalid points total (any year)
#' @param p Proportion of valid sightings recorded as uncertain
#' 
#'
#' @keywords simulations
#'
#' @export
#'
#' @examples
#'
#' x <- errorSims(n=10,pts=2,ipts=5,p=0.2)
#' 
#' ole <- spat.OLE(x$sightings[,1:4], x$grid[[2]], k.nn=7, N.nn=10, randomize=TRUE, adaptive=FALSE, reps=100, qualitycontrol = FALSE)
#' ole.qc <- spat.OLE(x$sightings[,1:4], x$grid[[2]], k.nn=7, N.nn=10, randomize=TRUE, adaptive=FALSE, reps=100, qualitycontrol = TRUE)
#' sb1 <- spat.SB(x$sightings[,1:4], x$grid[[2]], k.nn=5, N.nn=7, randomize=TRUE, reps=100, T.bound=60, model.num=1)
#' sb2 <- spat.SB(x$sightings[,1:4], x$grid[[2]], k.nn=5, N.nn=7, randomize=TRUE, reps=100, T.bound=60, model.num=2)
#' 
#' 
#' 
#' par(mfrow=c(1,1))
#' plot(r)
#' points(x$sightings[,1:2],col=ifelse(x$sightings$real==0,"blue",ifelse(x$sightings$sighting==2,"red","black")),pch=16,cex=1.3)
#' par(mfrow=c(2,2))
#' plot(ole$Mean.OLE)
#' plot(ole.qc$Mean.OLE)
#' plot(sb1$Mean.SB)
#' plot(sb2$Mean.SB)



errorSims <- function(n=10, pts=1, ipts=2, p=0.1) {
  
  options(warn=-1)
  
  r <- raster(matrix(nrow=n,ncol=n),xmn=0,xmx=n,ymn=0,ymx=n)
  for (i in 1:n) {
    r[i,] <- n+1-i
  }
  
  n <- n*5
  r <- r*5
  
  d <- data.frame(longitude=NA,latitude=NA,year=NA)
  
  for (i in 1:(n-1)) {
    tmp <- r>i
    tmp[tmp==0] <- NA
    x <- data.frame(randomPoints(tmp,n=pts))
    x$year <- i
    names(x) <- c('longitude','latitude','year')
    d <- rbind(d,x)
  }
  
  d <- na.omit(d)
  dx <- diff(c(xmin(r), xmax(r))) / ncol(r) / 2 # Half of horizontal width
  dy <- diff(c(ymin(r), ymax(r))) / nrow(r) / 2 # Half of vertical width
  xy <- coordinates(d[,1:2])                             # 2-column matrix of coordinates
  n2 <- nrow(xy)                                       # Number of sample points
  xy <- xy + c(runif(n2, -dx, dx), runif(n2, -dy, dy))  # Add random changes
  d[,1:2] <- xy
  d$sighting <- rbinom(nrow(d),1,p)+1
  
  d$real <-1 
  
  grid <- r*0
  
  
  d2 <- data.frame(longitude=c(1:ipts),latitude=c(1:ipts),year=NA,sighting=c(2))
  
  for (i in 1:ipts) {
    d2$longitude[i] <- runif(1,xmin(r),xmax(r))
    d2$latitude[i] <- runif(1,ymin(r),ymax(r))
    d2$year[i] <- runif(1,0,n)
  }
  
  d2$real <- 0
  
  d <- rbind(d,d2)
  
  options(warn=0)
  final <- list(grid=stack(r,grid),sightings=d)
  return(final)
}