R/plothistslice.R

In stranda/holoSimCell: Linked demographic/genetic simulation

###
### a plothist like function but one where you can specify a time point for the landscape history
### and plot all of the colonizations up to that point (as well as the habitat suitability at that point)
###
###

#' plot a slice of the landscape history
#'
#' Plots the history of colonization events from a spatially-explicit forward simulation of population expansion up to a specified time step in the simulation.
#'
#' @param timeslice the layer of the landscape that corresponds to a time step (generation) in the simulation
#' @param ph population history object (direct output from \code{gepophist2.cells})
#' @param landscape the landscape object that drives the simulation (see \code{ashSetupLandscape} for a full description of the landscape object)
#' @param window (optional, default = c(0,timeslice)) a two-element vector that defines the window of time for colonization events shown in the plot
#'
#' @details 
#' \itemize{
#' This function takes a fine-grained spatio-temporal history of a species generated by the forward-time simulator and produces a map that shows population spread through time (up to a specified time step). The input ph object includes the following named components:
#' \item{\code{pophist}} {a data frame with population ID, cell row, cell column, timing of colonization, and source population associated with colonists.}
#' \item{\code{Nvecs}} {a data frame with number of generations columns and number of populations rows that provides population size through time from the forward demographic simulation.}
#' \item{\code{tmat}} {a matrix of pairwise migration rates between populations, used in subsequent coalescent simulations.}
#' \item{\code{struct}} {a vector summarizing the structure of the simulated landscape from the forward simulation, x and y dimensions of the grid, number of generations, disperal parameters, etc.}
#' \item{\code{hab_suit}} {a matrix with number of generations rows and number of populations columns that provides habitat suitability through time.}
#' \item{\code{coalhist}} {a data frame specifying historical events for the coalescent simulation, includes time of event (time), source (src) and sink (snk) IDs, and fractional colonization (prop).}
#' \item{\code{popslst}} {a list object with one element per population from the forward-time demographic simulation.  Each element includes information on location on the landscape, colonization time, source populations, and population size through time.}
#' }
#'
#' @return
#' Does not return an object, but creates a plot illustrating the history of colonization of cells on the landscape over a specified time in the simulation.  Arrows indicate the direction of colonization and points are sized relative to population size in the forward simulation.  Habitat suitability (at generation \code{timeslice}) is represented by a color gradient (green as highest suitability).
#'
#' @examples
#' library(holoSimCell)
#' parms <- drawParms(control = system.file("extdata/ashpaper","Ash_priors.csv",package="holoSimCell"))
#' load(file=paste0(system.file(package="holoSimCell"),"/extdata/landscapes/",pollenPulls[[1]]$file))
#' refpops <- pollenPulls[[1]]$refs
#' avgCellsz <- mean(c(res(landscape$sumrast)))
#'
#' ph = getpophist2.cells(h = landscape$details$ncells, xdim = landscape$details$x.dim, ydim = landscape$details$y.dim,
#'                        landscape=landscape,
#'                        refs=refpops,   
#'                        refsz=parms$ref_Ne,
#'                        lambda=parms$lambda,
#'                        mix=parms$mix,  
#'                        shortscale=parms$shortscale*avgCellsz,  
#'                        shortshape=parms$shortshape, 
#'                        longmean=parms$longmean*avgCellsz,  
#'                        ysz=res(landscape$sumrast)[2], 
#'                        xsz=res(landscape$sumrast)[1], 
#'                        K = parms$Ne) 
#'
#' plotHistSlice(timeslice = 20, ph = ph, landscape = landscape)
#' plotHistSlice(timeslice = 200, ph = ph, landscape = landscape, window = c(150,200))
#'
#' @seealso \code{\link{getpophist2.cells}}, \code{\link{plothist}}, \code{\link{ashSetupLandscape}}
#'
#' @export
plotHistSlice <- function(timeslice,ph,landscape,window=c(0,timeslice))
{
    if ((nrow(landscape$hab_suit)<timeslice) |( timeslice<1))
        stop ("need to specify a timeslice within the range of time clicks for this landscape")
    r <- landscape$sumrast
    values(r) <- landscape$hab_suit[timeslice,]
    pophist <- ph$pophist
    nv <- ph$Nvec
    xy <- data.frame(xyFromCell(r,unique(sort(pophist$pop))))
    names(xy) <- c("plotx","ploty")
    xy$pop <- unique(sort(pophist$pop))
    xy$ploty <- c(t(t(matrix(xy$ploty,nrow=dim(r)[2],ncol=dim(r)[1]))[dim(r)[1]:1,]))
    pophist <- merge(pophist,xy,att.x=T)
    pophist <- pophist[((!is.na(pophist$arrive)) & (pophist$arrive<=timeslice)),]
    pophist <- pophist[order(pophist$arrive),]

    opar=par()
    par(mar=c(2,2,2,0.5)+0.1)
    sp::plot(flip(r,2), legend=F)

    szfunc <- function(x) {log(x+1)/6}
    
    for (i in 1:nrow(pophist)) #plot population locations
    {
        
        points(x=pophist$plotx[i],y=pophist$ploty[i],pch=16,col=1,cex=szfunc(nv[pophist$pop[i],timeslice]))
    }

    for (i in 1:nrow(pophist)) #plot arrows
        if ((pophist$arrive[i]>=min(window)) & (pophist$arrive[i]<=max(window)))
        {
        x1=pophist$plotx[i]
        y1=pophist$ploty[i]

        x0=pophist$plotx[pophist$pop==pophist$source[i]]
        y0=pophist$ploty[pophist$pop==pophist$source[i]]

        arrows(x0,y0,x1,y1, length=0.07, angle=15, code=2)
        }
    suppressWarnings(par(opar))
}