R/plotRD.r

In HMMoce: Improved Analysis of Marine Animal Movement Data Using Hidden Markov Models

#' Plot RD results
#'
#' \code{plotRD} uses HMM output via \code{calc.track} to calculate and plot residency distributions for each behavior state and combined data
#'
#' @param distr is output array from \code{hmm.smoother}
#' @param track is output dataframe from \code{calc.track}
#' @param ptt is individual identification number
#' @param known is 3 column data frame containing date, lat, lon of known
#'   movement track. This is only useful for comparing HMMoce results to known
#'   track collected by SPOT or GPS, for example. Default is NULL.
#' @param g grid from \code{\link{setup.grid}}
#' @param xlims a vector of length 2 indicating longitude limits (-180 to 180).
#' @param ylims a vector of length 2 indicating latitude limits
#' @param makePlot is logical indicating whether to make a plot of the RD
#' @param save.plot is logical indicating whether you want the plot written to
#'   disk using \code{pdf}.
#' @importFrom fields "world"
#'
#' @return a list of one raster layer for the combined RD and one raster brick (2 layers) for the individual behavior RDs.
#' @export

plotRD <- function(distr, track, ptt, known=NULL, g, xlims, ylims, makePlot = TRUE, save.plot=FALSE){
  
  crs <- "+proj=longlat +datum=WGS84 +ellps=WGS84"
  rd.cols <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, 'RdYlGn')))
  
  p.1 <- apply(distr[1,,,], 1, sum)
  p.2 <- apply(distr[2,,,], 1, sum)
  sv <- -(apply(distr[1,,,], 1, sum) > apply(distr[2,,,], 1, sum)) + 1
  #sv[sv == 0] <- NA
  
  # normalize for each behav at each time point
  norm <- array(NA, dim=dim(distr)[c(3,4,2,1)])
  for (i in 1:dim(norm)[3]){
    norm[,,i,1] <- (distr[1,i,,] / max(distr[1,i,,], na.rm=T)) * p.1[i]
    norm[,,i,2] <- (distr[2,i,,] / max(distr[2,i,,], na.rm=T)) * p.2[i]
  }
  
  all <- apply(norm, 1:2, FUN=function(x) sum(x, na.rm=T))
  all <- all / max(all, na.rm=T)
  
  # sum across time and normalize final surface
  norm <- apply(norm, c(1:2,4), FUN=function(x) sum(x, na.rm=T))
  norm[,,1] <- norm[,,1] / max(norm[,,1], na.rm=T)
  norm[,,2] <- norm[,,2] / max(norm[,,2], na.rm=T)
  
  # get pdf as raster
  norm <- raster::flip(raster::brick(norm, xmn=min(g$lon), xmx=max(g$lon),
                                     ymn=min(g$lat), ymx=max(g$lat), crs=crs, transpose=T), 2)
  all <- raster::flip(raster::raster(t(all), xmn=min(g$lon), xmx=max(g$lon),
                                     ymn=min(g$lat), ymx=max(g$lat), crs=crs), 2)
  
  if(makePlot){
    # set up plot brks and cols
    #norm.rnge <- cellStats(norm, 'range')
    norm.breaks = seq(.001, 1, length.out=201)
    norm.mid = norm.breaks[1:(length(norm.breaks)-1)]
    plot.rd.col <- rd.cols(length(norm.breaks)-1)
    
    # set axes
    x.at <- pretty(xlims, 5)
    x.labels <- parse(text=paste(abs(x.at), "*degree~W", sep=""))
    y.at <- pretty(ylims, 5)
    y.labels <- parse(text=paste(abs(y.at), "*degree~N", sep=""))
    
    # set margins
    old.mar <- par()$mar
    mar.default <- par('mar')
    mar.bottom <- mar.default
    mar.bottom[1] <- 0
    
    # get some plot indices
    idx1 <- which((sv+1) == 1)
    idx2 <- which((sv+1) == 2)
    TT <- length(track$lon)
    
    if(save.plot) grDevices::pdf(paste(ptt, '_RD.pdf', sep = ''), width = 7, height = 12)
    
    # BUILD IT
    nf <- layout(matrix(c(1,2,
                          3,4,
                          5,6), 3, 2, byrow=T), widths=c(7,2), heights=c(4,4,4))
    #layout.show(nf)
    
    # COMBINED PANEL
    par (mar=mar.bottom)
    
    raster::image(all, maxpixels=raster::ncell(all), xlim=xlims, ylim=ylims, axes=F,
                  col=plot.rd.col, breaks=norm.breaks, xlab='', ylab='')#, main='All')#, zlim=c(.05,1))#, axes=F)#, main=dt.idx) #, breaks=zbreaks
    axis(1, at=x.at, labels=FALSE)
    axis(2, at=y.at, labels=y.labels);
    world(add=T, fill=T, col='grey60')
    box()
    if(!is.null(known)) lines(known$lon, known$lat)
    
    # COLORBAR
    image(1, norm.mid, t(as.matrix(norm.mid)), breaks=norm.breaks, col=plot.rd.col, axes=FALSE, xlab="",
          ylab='Expected Proportion of Time Spent')
    axis(2, at=c(.05, .25, .5, .75, .95), labels=paste(rev(c(5, 25, 50, 75, 95)), '%'));box();
    
    # MIGRATORY PANEL
    par (mar=mar.bottom)
    raster::image(norm[[1]], maxpixels=raster::ncell(norm[[1]]), xlim=xlims, ylim=ylims, axes=F,
                  col=plot.rd.col, breaks=norm.breaks, xlab='', ylab='')#, main='Migratory')#, zlim=c(.05,1))#, axes=F)#, main=dt.idx) #, breaks=zbreaks
    axis(1, at=x.at, labels=FALSE)
    axis(2, at=y.at, labels=y.labels)
    world(add=T, fill=T, col='grey60')
    box()
    mtext('Migratory behavior', 4, line=2)
    graphics::points(track$lon[idx1], track$lat[idx1], bg='white', pch=21, cex=1.25)
    graphics::points(track$lon[1], track$lat[1], bg = 'green', pch = 21, cex=1.75)
    graphics::points(track$lon[TT], track$lat[TT], bg = 'red', pch = 21, cex=1.75)
    
    # FILL COLORBAR SPACE
    plot.new()
    
    # RESIDENT PANEL
    par (mar=mar.default)
    raster::image(norm[[2]], maxpixels=raster::ncell(norm[[2]]), xlim=xlims, ylim=ylims, axes=F,
                  col=plot.rd.col, breaks=norm.breaks, xlab='', ylab='')#, main='Resident')#, zlim=c(.05,1))#, axes=F)#, main=dt.idx) #, breaks=zbreaks
    axis(1, at=x.at, labels=x.labels)
    axis(2, at=y.at, labels=y.labels);
    world(add=T, fill=T, col='grey60')
    box()
    mtext('Resident behavior', 4, line=2)
    graphics::points(track$lon[idx2], track$lat[idx2], bg='grey60', pch=21, cex=1.25)
    graphics::points(track$lon[1], track$lat[1], bg = 'green', pch = 21, cex=1.75)
    graphics::points(track$lon[TT], track$lat[TT], bg = 'red', pch = 21, cex=1.75)
    
    # FILL COLORBAR SPACE
    plot.new()
    
    if(save.plot) grDevices::dev.off()
  }
  
  return(list(allRD = all, behavRD = norm))
  
}