R/walker.R

In MAnalytics/stppSim: Spatiotemporal Point Patterns Simulation

#' @title A landscape walker
#' @description A dynamic object capable of moving
#' and avoiding obstacles on a landscape.
#' @param n number of events
#' to be generated by a walker within a
#' temporal bin.
#' @param s_threshold defines the spatial
#' perception range of a walker at a given
#' location. Default: \code{250} (in the same
#' linear unit
#' as the `poly` - polygon shapefile).
#' @param step_length the maximum step taken
#' by a walker from one point to the next.
#' @param poly (An sf or S4 object)
#' a polygon shapefile defining the
#' extent of the landscape
#' @param restriction_feat (An S4 object) optional
#' shapefile containing features
#' in which walkers cannot walk through.
#' Default: \code{NULL}.
#' @param field a number in the range of \code{[0-1]}
#' (i.e. restriction values) assigned
#' to all features; or
#' the name of a numeric field to extract such
#' restriction values for different classes of
#' feature.
#' Restriction value `0` and `1` indicate the
#' lowest and the highest obstructions, respectively.
#' Default: \code{NULL}.
#' @param coords a vector of the form c(x, y) giving the
#' initial coordinates of a walker (i.e., coordinates
#' of origins).
#' Default value is \code{c(0,0)} for an
#' arbitrary square space.
#' @param pt_itx To check whether any of the
#' specified initial origin coordinates
#' falls outside the boundary.
#' Default: \code{TRUE}.
#' @param show.plot (TRUE or False) To show the time series
#' plot. Default is \code{FALSE}.
#' @usage walker(n = 5, s_threshold = 250, step_length = 20,
#' poly = NULL, restriction_feat=NULL, field = NA, coords=c(0,0),
#' pt_itx = TRUE, show.plot = FALSE)
#' @examples
#' #load boundary of Camden
#' load(file = system.file("extdata", "camden.rda",
#' package="stppSim"))
#' boundary = camden$boundary # get boundary
#' walkerpath <- walker(n = 5, s_threshold = 250, step_length = 20,
#' poly = boundary, restriction_feat=NULL, field = NULL,
#' coords = c(0,0), pt_itx = TRUE, show.plot = FALSE)
#' #plot(walkerpath)
#' @details
#' A walker is propelled by an in-built stochastic
#' transition matrix
#' and a specified set of spatial and temporal
#' parameters. The transition
#' matrix defines two states, namely; the `exploratory`
#' and a `performative` states. A walker is capable
#' of avoiding obstructions (i.e., `restriction_feat`)
#' if included. The resulting number of events may be
#' slightly different from the value `n` because of the
#' stochastic process involved.
#' @return Returns a trace of walker's path, and the
#' resulting events.
#' @references
#' Quaglietta L, Porto M (2019). SiMRiv: Individual-Based,
#' Spatially-Explicit Simulation and Analysis of Multi-State
#' Movements in River Networks and Heterogeneous Landscapes_.
#' R package version 1.0.4, <URL:
#' https://CRAN.R-project.org/package=SiMRiv>.
#' @importFrom dplyr select filter
#' @importFrom SiMRiv species transitionMatrix
#' state.CRW simulate
#' @importFrom chron chron
#' @importFrom stats time
#' @importFrom sf st_intersects st_as_sf st_centroid
#' st_point
#' @importFrom stringr str_split str_remove_all
#' @export

walker <- function(n = 5, s_threshold = 250, step_length = 20,
                   poly=NULL, restriction_feat=NULL,
                   field = NA,
                   coords = c(0,0), pt_itx = TRUE, show.plot = FALSE){

  #output holder
  output <- list()

  points <- text <- sn <- x <- y <- X3 <- NULL

  #test polygon geometry
  if(!is.null(poly)){
    #-----
    poly_tester(poly)
    #-----
  }

  if(pt_itx == TRUE){
    if(is.null(poly)  & coords[1]!=0 & coords[2]!=0){
      stop("`poly` argument is NULL!")
    }

  #poly defined but not coords, then get centroid
  if(!is.null(poly) & coords[1]==0 & coords[2]==0){ #extract centroid coords
    cent_xy <- as.character(st_centroid(st_as_sf(poly)$geometry))
    cent_xy <- str_remove_all(cent_xy, "[(c,)]")
    cent_xy_x <-  as.numeric(str_split(cent_xy, " ", 2)[[1]][1])
    cent_xy_y <-  as.numeric(str_split(cent_xy, " ", 2)[[1]][2])
    cent_xy <- c(cent_xy_x, cent_xy_y)
    #set coords
    coords <- cent_xy
    }
  }

  #configure walker
  # Walker <- species(
  # state.CRW(0.005) + state.CRW(0.99),
  # transitionMatrix(exp(-6.25679 + 1.26863*log(n)), #
  #                  0.70))

  Walker <- species(
    state.CRW(0.005) + state.CRW(0.99),
    transitionMatrix(exp(-6.25679 + 1.26863*log(n)), #
                     0.70))

  #-6.25679 + 1.26863*log(n) is the
  #power regression that relate x and y (see 'calibra..R')

  #define resolution (must be greater than step_length)
  min_res <- ((step_length * 2) + (500/20)) + 10#plus some values

  #create the landscape restriction raster
  #create boundary
  #test polygon geometry

  if(!is.null(poly)){
    landscape <- space_restriction(shp = poly,
                                   res = min_res, binary = TRUE)

    #update the landscape restriction map
    if(!is.null(restriction_feat)){
      landscape <- space_restriction(shp = restriction_feat,
                                     baseMap = landscape,
                                    res = min_res, field = field)
    }
    #plot(landscape)
    #check point polygon intersection
    if(pt_itx == TRUE){
      st_int <- st_intersects(st_as_sf(poly),
                              st_point(coords, dim="XY"))
      st_int_yes <- as.numeric(st_int)
    }
  }

  #plot(landscape)
  #does the point fall within the boundary?

  #plot(poly)
  #plot(st_point(x = coords, dim = "XY"))

  #meaning 1-step/hrs
  #s_threshold <- 1000
  Walker <- (Walker + (step_length + (s_threshold/10))) #20 assumed value
  #Walker <- (Walker + 20) * s_threshold

  if(is.null(landscape)){
    sim <- simulate(Walker, time=200, coords)
    #plot(sim[,1], sim[,2])
  }

  if(!is.null(landscape)){
    sim <- simulate(Walker, time=200, resist = landscape, coords)
    #plot(landscape)
    #points(sim[,1], sim[,2], add=TRUE)
  }

  #200 is the no of time.steps to be simulated
  #extract event locations
  sim_events <- data.frame(sim)
  colnames(sim_events) <- c("X1","X2","X3")
  sim_events <- sim_events %>%
    filter(X3 == 1)
  #200 (number of steps per origin)..
  #was used in the calibration ##nrow(sim) ##length(which(sim[,3]==1))

  #create hour sequence
  # hourly unit
  hm <- merge(0:23, seq(0, 0, by = 0))
  hour_seq <- chron(times. = paste(hm$x, ':', hm$y, ':', 0))

    #if event is present
    if(nrow(sim_events) >= 1){

      #assign random (but progressing) time to events
      hr_sample <- sample(hour_seq, nrow(sim_events), replace = TRUE)
      hr_sample <- hr_sample[order(hr_sample)]

      sim_events_ <- cbind(sim_events, 1:nrow(sim_events), hr_sample)

      if(show.plot == TRUE){
        plot(sim, type="l", asp=1, col="gray80")
        points(sim_events_, col="red")
        text(sim_events_[,1], sim_events_[,2],
         labels=sim_events_[,4], cex= 0.7, pos=3)
      }

      colnames(sim_events_) <- c("x","y","yes","sn", "time")

      sim_events_ <- data.frame(sim_events_) %>%
        dplyr::select(sn, x, y, time)


    #output$intersection <- st_int_yes  #poly-point intersection
    #output$p_events <- sim_events_
    output <- sim_events_

    return(output)
    }

}