R/plot_crwData.R

In momentuHMM: Maximum Likelihood Analysis of Animal Movement Behavior Using Multivariate Hidden Markov Models

#' Plot \code{crwData}
#' 
#' Plot observed locations, error ellipses (if applicable), predicted locations, and prediction intervals from \code{\link{crwData}} or \code{\link{crwHierData}} object.
#'
#' @param x An object \code{crwData} or \code{crwHierData} (as returned by \code{\link{crawlWrap}}).
#' @param animals Vector of indices or IDs of animals for which information will be plotted.
#' Default: \code{NULL} ; all animals are plotted.
#' @param compact \code{TRUE} for a compact plot (all individuals at once), \code{FALSE} otherwise
#' (default -- one individual at a time). Ignored unless \code{crwPredictPlot=FALSE}.
#' @param ask If \code{TRUE}, the execution pauses between each plot.
#' @param plotEllipse If \code{TRUE} (the default) then error ellipses are plotted (if applicable). Ignored unless \code{crwPredictPlot=FALSE}.
#' @param crawlPlot Logical indicating whether or not to create individual plots using \code{\link[crawl]{crwPredictPlot}}. See \code{\link[crawl]{crwPredictPlot}} for details.
#' @param ... Further arguments for passing to \code{\link[crawl]{crwPredictPlot}}
#' 
#' @details 
#' In order for error ellipses to be plotted, the names for the semi-major axis, semi-minor axis, and 
#' orientation in \code{x$crwPredict} must respectively be \code{error_semimajor_axis}, \code{error_semiminor_axis},     
#' and \code{error_ellipse_orientation}. 
#' 
#' If the \code{crwData} (or \code{crwHierData}) object was created using data generated by
#' \code{\link{simData}} (or \code{\link{simHierData}}) or \code{\link{simObsData}}, then the true locations (\code{mux},\code{muy}) are also plotted.
#' 
#' @seealso \code{\link[crawl]{crwPredictPlot}}
#'
#' @examples
#' \dontrun{
#' # extract simulated obsData from example data
#' obsData <- miExample$obsData
#' 
#' # error ellipse model
#' err.model <- list(x= ~ ln.sd.x - 1, y =  ~ ln.sd.y - 1, rho =  ~ error.corr)
#'
#' # create crwData object
#' crwOut <- crawlWrap(obsData=obsData,
#'          theta=c(4,0),fixPar=c(1,1,NA,NA),
#'          err.model=err.model)
#'
#' plot(crwOut,compact=TRUE,ask=FALSE,plotEllipse=FALSE)
#' }
#'
#' @importFrom grDevices gray adjustcolor
#' @importFrom graphics abline axis hist mtext par plot points legend title
# @importFrom conicfit calculateEllipse
#' @importFrom crawl crwPredictPlot

#' @method plot crwData
#' @aliases plot.crwHierData
#' @export
plot.crwData <- function(x,animals=NULL,compact=FALSE,ask=TRUE,plotEllipse=TRUE,crawlPlot=FALSE,...)
{
  # check arguments
  if(plotEllipse){
    if (!requireNamespace("conicfit", quietly = TRUE)) {
      stop("Package \"conicfit\" needed for this function to work. Please install it.",
           call. = FALSE)
    }
  }
  if(length(x)<1 | any(dim(x)<1))
    stop("The data input is empty.")
  
  if(is.null(x$crwFits) |  is.null(x$crwPredict))
    stop("Missing field(s) in data.")
  
  data <- x$crwPredict
  model_fits <- x$crwFits
  
  convFits<-which(unlist(lapply(model_fits,function(x) inherits(x,"crwFit"))))
  if(!length(convFits)) stop('crwData object contains no valid crawl::crwMLE model fits')
  model_fits<-model_fits[convFits]
  
  coordNames <- model_fits[[1]]$coord
  data$x<-data[[coordNames[1]]]
  data$y<-data[[coordNames[2]]]
  
  nbAnimals <- length(unique(data$ID))
  
  ##################################
  ## Define animals to be plotted ##
  ##################################
  if(is.null(animals)) # all animals are plotted
    animalsInd <- 1:nbAnimals
  else {
    if(is.character(animals)) { # animals' IDs provided
      animalsInd <- NULL
      for(zoo in 1:length(animals)) {
        if(length(which(unique(data$ID)==animals[zoo]))==0) # ID not found
          stop("Check animals argument.")
        
        animalsInd <- c(animalsInd,which(unique(data$ID)==animals[zoo]))
      }
    }
    
    if(is.numeric(animals)) { # animals' indices provided
      if(length(which(animals<1))>0 | length(which(animals>nbAnimals))>0) # index out of bounds
        stop("Check animals argument.")
      
      animalsInd <- animals
    }
  }
  
  if(nbAnimals>7) {
    # to make sure that all colours are distinct (emulate ggplot default palette)
    hues <- seq(15, 375, length = nbAnimals + 1)
    colors <- hcl(h = hues, l = 65, c = 100)[1:nbAnimals]
  } else
    colors <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
  
  # graphical parameters
  par(mar=c(5,4,4,2)-c(0,0,2,1)) # bottom, left, top, right
  par(ask=ask)
  
  if(compact & !crawlPlot) {
    ################################
    ## Map of all animals' tracks ##
    ################################
    par(mfrow = c(1,1))
    
    # determine bounds
    ind <- which(data$ID %in% unique(data$ID)[animalsInd])
    x <- data$x[ind]
    y <- data$y[ind]
    mu.x <- data$mu.x[ind]
    mu.y <- data$mu.y[ind]
    truex <- data$mux[ind]
    truey <- data$muy[ind]
    
    xmin <- min(c(x,mu.x,truex),na.rm=T)
    xmax <- max(c(x,mu.x,truex),na.rm=T)
    ymin <- min(c(y,mu.y,truey),na.rm=T)
    ymax <- max(c(y,mu.y,truey),na.rm=T)
    # make sure that x and y have same scale
    if(xmax-xmin>ymax-ymin) {
      ymid <- (ymax+ymin)/2
      ymax <- ymid+(xmax-xmin)/2
      ymin <- ymid-(xmax-xmin)/2
    } else {
      xmid <- (xmax+xmin)/2
      xmax <- xmid+(ymax-ymin)/2
      xmin <- xmid-(ymax-ymin)/2
    }
    
    ID <- unique(data$ID)[animalsInd[1]]
    ind <- which(data$ID==ID)
    x <- data$x[ind]
    y <- data$y[ind]
    mu.x <- data$mu.x[ind]
    mu.y <- data$mu.y[ind]
    truex <- data$mux[ind]
    truey <- data$muy[ind]

    # plot the first animal's track
    if(!is.null(truex) & !is.null(truey)){
      if(plotEllipse & !is.null(data$error_semimajor_axis) & !is.null(data$error_semiminor_axis) & !is.null(data$error_ellipse_orientation)){
        s<-which(!is.na(x))[1]
        errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
        plot(errorEllipse[,1],errorEllipse[,2],type="l",xlab="x",ylab="y",xlim=c(xmin,xmax),ylim=c(ymin,ymax),lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1),asp=1)
        for(s in which(!is.na(x))[-1]){
          errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
          lines(errorEllipse[,1],errorEllipse[,2],lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1))
        }
        points(truex,truey,type="l",lwd=3,col=gray(.75))
      } else {
        plot(truex,truey,type="l",lwd=3,xlab="x",ylab="y",xlim=c(xmin,xmax),ylim=c(ymin,ymax),col=gray(.75),asp=1)
      }
      points(mu.x,mu.y,type="o",lwd=1.3,xlab="x",ylab="y",pch=20,col=colors[animalsInd[1]],cex=0.5)
    } else {
      if(plotEllipse & !is.null(data$error_semimajor_axis) & !is.null(data$error_semiminor_axis) & !is.null(data$error_ellipse_orientation)){
        s<-which(!is.na(x))[1]
        errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
        plot(errorEllipse[,1],errorEllipse[,2],type="l",xlab="x",ylab="y",xlim=c(xmin,xmax),ylim=c(ymin,ymax),lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1),asp=1)
        for(s in which(!is.na(x))[-1]){
          errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
          lines(errorEllipse[,1],errorEllipse[,2],lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1))
        }
        points(mu.x,mu.y,type="o",lwd=1.3,pch=20,col=colors[animalsInd[1]],cex=0.5)
      } else {
        plot(mu.x,mu.y,type="o",lwd=1.3,xlab="x",ylab="y",pch=20,xlim=c(xmin,xmax),ylim=c(ymin,ymax),col=colors[animalsInd[1]],cex=0.5,asp=1)
      }
    }
    points(x[!is.na(x)],y[!is.na(x)],xlab="x",ylab="y",pch=20,cex=0.5)
    
    # add each other animal's track to the map
    for(zoo in animalsInd[-1]) {
      ID <- unique(data$ID)[zoo]
      ind <- which(data$ID==ID)
      x <- data$x[ind]
      y <- data$y[ind]
      mu.x <- data$mu.x[ind]
      mu.y <- data$mu.y[ind]
      truex <- data$mux[ind]
      truey <- data$muy[ind]
      
      if(plotEllipse & !is.null(data$error_semimajor_axis) & !is.null(data$error_semiminor_axis) & !is.null(data$error_ellipse_orientation)){
        for(s in which(!is.na(x))){
          errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
          lines(errorEllipse[,1],errorEllipse[,2],lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1))
        }
      }
      if(!is.null(truex) & !is.null(truey)) points(truex,truey,type="l",lwd=3,col=gray(.75))
      points(mu.x,mu.y,type="o",lwd=1.3,pch=20,col=colors[zoo],cex=0.5)
      points(x[!is.na(x)],y[!is.na(x)],xlab="x",ylab="y",pch=20,cex=0.5)
    }
    
    legendtxt<-c("observed")
    legendcol<-c(1)
    legendpch<-c(20)
    legendlty<-c(NA)
    legendlwd<-c(NA)
    if(!is.null(truex) & !is.null(truey)){
      legendtxt<-c(legendtxt,"true")
      legendcol<-c(legendcol,gray(.75))
      legendpch<-c(legendpch,NA)
      legendlty<-c(legendlty,1)
      legendlwd<-c(legendlwd,3)
    }
    legendtxt<-c(legendtxt,unique(data$ID)[animalsInd])
    legendcol<-c(legendcol,colors[animalsInd])
    legendpch<-c(legendpch,rep(NA,length(animalsInd)))
    legendlty<-c(legendlty,rep(1,length(animalsInd)))
    legendlwd<-c(legendlwd,rep(3,length(animalsInd)))

    legend("topright",legendtxt,col=legendcol,pch=legendpch,lty=legendlty,lwd=legendlwd,bty="n")
    
  } else {
    for(zoo in animalsInd) {
      
      ID <- unique(data$ID)[zoo]
      ind <- which(data$ID==ID)
        
      if(!crawlPlot){
        
        x <- data$x[ind]
        y <- data$y[ind]
        mu.x <- data$mu.x[ind]
        mu.y <- data$mu.y[ind]
        truex <- data$mux[ind]
        truey <- data$muy[ind]
        
        ################################
        ## Map of each animal's track ##
        ################################
        par(mfrow=c(1,1))
        # determine bounds
        xmin <- min(c(x,mu.x,truex),na.rm=T)
        xmax <- max(c(x,mu.x,truex),na.rm=T)
        ymin <- min(c(y,mu.y,truey),na.rm=T)
        ymax <- max(c(y,mu.y,truey),na.rm=T)
        # make sure that x and y have same scale
        if(xmax-xmin>ymax-ymin) {
          ymid <- (ymax+ymin)/2
          ymax <- ymid+(xmax-xmin)/2
          ymin <- ymid-(xmax-xmin)/2
        } else {
          xmid <- (xmax+xmin)/2
          xmax <- xmid+(ymax-ymin)/2
          xmin <- xmid-(ymax-ymin)/2
        }
        # map of the animal's track
        
        legendtxt<-c("observed")
        legendcol<-1
        legendpch<-20
        legendlty<-NA
        legendlwd<-NA
        if(!is.null(truex) & !is.null(truey)){
          if(plotEllipse & !is.null(data$error_semimajor_axis) & !is.null(data$error_semiminor_axis) & !is.null(data$error_ellipse_orientation)){
            s<-which(!is.na(x))[1]
            errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
            plot(errorEllipse[,1],errorEllipse[,2],type="l",xlab="x",ylab="y",xlim=c(xmin,xmax),ylim=c(ymin,ymax),lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1),asp=1)
            for(s in which(!is.na(x))[-1]){
              errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
              lines(errorEllipse[,1],errorEllipse[,2],lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1))
            }
            points(truex,truey,type="l",lwd=3,col=gray(.75))
          } else {
            plot(truex,truey,type="l",lwd=3,xlab="x",ylab="y",xlim=c(xmin,xmax),ylim=c(ymin,ymax),col=gray(.75),asp=1)
          }
          points(mu.x,mu.y,type="o",lwd=1.3,pch=20,col=colors[zoo])
          legendtxt<-c(legendtxt,"true")
          legendcol<-c(legendcol,gray(.75))
          legendpch<-c(legendpch,NA)
          legendlty<-c(legendlty,1)
          legendlwd<-c(legendlwd,3)
        } else {
          if(plotEllipse & !is.null(data$error_semimajor_axis) & !is.null(data$error_semiminor_axis) & !is.null(data$error_ellipse_orientation)){
            s<-which(!is.na(x))[1]
            errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
            plot(errorEllipse[,1],errorEllipse[,2],type="l",xlab="x",ylab="y",xlim=c(xmin,xmax),ylim=c(ymin,ymax),lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1),asp=1)
            for(s in which(!is.na(x))[-1]){
              errorEllipse<-getEllipse(x[s],y[s],data$error_semimajor_axis[ind][s],data$error_semiminor_axis[ind][s],data$error_ellipse_orientation[ind][s])
              lines(errorEllipse[,1],errorEllipse[,2],lty=1,cex=0.5,col=adjustcolor(1,alpha.f=0.1))
            }
            points(mu.x,mu.y,type="o",lwd=1.3,pch=20,col=colors[zoo])
          } else {
            plot(mu.x,mu.y,type="o",lwd=1.3,xlab="x",ylab="y",pch=20,xlim=c(xmin,xmax),ylim=c(ymin,ymax),col=colors[zoo],asp=1)
          }
        }
        points(x[!is.na(x)],y[!is.na(x)],pch=20,cex=0.5)

        mtext(paste("ID",ID),side=3,outer=TRUE,padj=2)
        legendtxt<-c(legendtxt,"predicted")
        legendcol<-c(legendcol,colors[zoo])
        legendpch<-c(legendpch,NA)
        legendlty<-c(legendlty,1)
        legendlwd<-c(legendlwd,3)
        legend("topright",legendtxt,col=legendcol,pch=legendpch,lty=legendlty,lwd=legendlwd,bty="n")
      } else {
        crawl::crwPredictPlot(data[ind,],...)
        titleInd <- TRUE
        if(!missing(...)) {
          if("main" %in% names(list(...))) titleInd <- FALSE
        }
        if(titleInd) title(main=paste("ID",ID))
      }
    }      
  }
  

  
  # set graphical parameters back to default
  par(ask=FALSE)
  par(mfrow=c(1,1))
  par(mar=c(5,4,4,2))
}

#' @method plot crwHierData
#' @export
plot.crwHierData <- function(x,animals=NULL,compact=FALSE,ask=TRUE,plotEllipse=TRUE,crawlPlot=FALSE,...)
{
  data <- x
  data$crwPredict <- data$crwPredict[which(x$crwPredict$level==attr(x$crwPredict,"coordLevel")),]
  class(data$crwPredict) <- append("crwPredict",class(data$crwPredict))
  plot(crwData(data),animals=animals,compact=compact,ask=ask,plotEllipse=plotEllipse,crawlPlot=crawlPlot,...)
}

getEllipse<-function(x,y,M,m,r){
  conicfit::calculateEllipse(x,y,M,m,90-r) #calculateEllipse angle runs counterclockwise from east but Argos error ellipse runs clockwise from north
}