manuscript/vignette_examples.R

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

library(sp)
library(doParallel)

example_wd <- ("~/Documents/Dropbox/current projects/moveHMM extension/momentuHMM/vignette examples/")

append.RData <- function(x, file) {
  old.objects <- load(file)
  save(list = c(old.objects, deparse(substitute(x))), file = file)
}

###################################################
### Elephant example
###################################################
#source(paste0(getwd(),"/elephantExample.R"))
load(paste0(example_wd,"elephantExample.RData"))

activityBudgets<-table(viterbi(m3))/nrow(m3$data)

save(activityBudgets,file=paste0(getwd(),"/vignette_inputs.RData"))

#latlongDat<-as.data.frame(elephantData)
#sp::coordinates(latlongDat)<-c("x","y")
#sp::proj4string(latlongDat) = sp::CRS("+proj=utm +zone=30N +datum=WGS84" )
#latlongDat<-sp::spTransform(latlongDat,sp::CRS("+proj=longlat +datum=WGS84"))
#satPlotStates<-plotSat(as.data.frame(latlongDat),zoom=8,location=c(median(rawData$lon),median(rawData$lat)),ask=FALSE,return=TRUE,states=viterbi(m3),col=c("#009E73", "#F0E442"),stateNames=c("encamped","exploratory"))

acfLag<-24*14
pdf(file=paste0(getwd(),"/plot_elephantResults%03d.pdf"),onefile=FALSE)
plot(m3,ask=FALSE,plotCI=TRUE,covs=data.frame(hour=12),cex.main=1.4,cex.lab=1.4,cex.axis=1.4,cex.legend=1.4)
acf(pr[["stepRes"]], lag.max = acfLag, na.action = na.pass, xlab="Lag (hours)", main = "",cex.lab=1.4,cex.axis=1.4)
acf(elephantData$step,na.action=na.pass,lag=acfLag,main="",xlab="Lag (hours)",cex.lab=1.4,cex.axis=1.4)
dev.off()

for(plt in seq(1,17)[-c(1,2,9,10,12,13,15,16,17)])
  unlink(paste0("plot_elephantResults0",ifelse(plt>9,"","0"),plt,".pdf"))

png(filename="elephant_plotSat.png",width=6,height=6,units="in",res=90)
plotSat(data.frame(x=rawData$lon,y=rawData$lat),zoom=8,location=c(median(rawData$lon),median(rawData$lat)),ask=FALSE)
dev.off()

rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])

###################################################
### Northern fur seal example
###################################################
#source(paste0(getwd(),"/nfsExample.R"))
load(paste0(example_wd,"nfsExample.RData"))

nfsTimeInStates<-nfsFits$miSum$Par$timeInStates

append.RData(nfsTimeInStates,file=paste0(getwd(),"/vignette_inputs.RData"))

pdf(file=paste0(getwd(),"/plot_nfsResults.pdf"))
plot(nfsFits,ask=FALSE,legend.pos="topright",cex.main=1.3,cex.lab=1.3,cex.axis=1.3,cex.legend=1.3)
dev.off()
rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])

###################################################
### Turtle example
###################################################
#source(paste0(getwd(),"/turtleExample.R"))
load(paste0(example_wd,"turtleExample.RData"))

turtle_miSum<-turtleFits$miSum[c("Par","data")]

append.RData(turtle_miSum,file=paste0(getwd(),"/vignette_inputs.RData"))

library(TeachingDemos)
library(raster)

speedrast<-speedBrick[["X2012.12.02"]]
speedrast<-projectRaster(speedrast,crs="+init=epsg:26718 +proj=utm +zone=18 +datum=NAD27 +units=km +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")
dirrast<-dirBrick[["X2012.12.02"]]
dirrast<-projectRaster(dirrast,crs="+init=epsg:26718 +proj=utm +zone=18 +datum=NAD27 +units=km +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")

rastxy<-coordinates(speedrast)
rastx0<-unique(rastxy[,1])
rasty0<-unique(rastxy[,2])
rastx<-rep(rastx0,length(rasty0))
rasty<-rep(rasty0,each=length(rastx0))

x<-turtleFits$miSum$data$x/1000
y<-turtleFits$miSum$data$y/1000

pdf(file=paste0(getwd(),"/plot_turtleResults%03d.pdf"),onefile=FALSE)
turtleFits$miSum$stateNames<-c("foraging","transit")
plot(turtleFits,plotCI=TRUE,covs=data.frame(angle_osc=1),ask=FALSE,cex.main=2,cex.lab=1.7,cex.axis=1.65)
dev.off()

for(plt in seq(1,6)[-c(1,2,4)])
  unlink(paste0("plot_turtleResults0",ifelse(plt>9,"","0"),plt,".pdf"))


pdf(file=paste0(getwd(),"/plot_turtleResults2.pdf"),width=7,height=7*50/73)
par(mar=c(4,4,0,1))
plot(speedrast,col=gray.colors(20, start=1, end=0.3),xlab="easting (km)",ylab="northing (km)",legend.args=list(text='       speed (m/s)', side=3, line=1),box=FALSE,bty="n")
my.symbols(rastx,rasty,ms.arrows, xsize=20,ysize=20,add=TRUE,angle=getValues(dirrast), r=getValues(speedrast), length=.015)
points(x,y,type="o",pch=20,col=c("#E69F00", "#56B4E9")[turtleFits$miSum$Par$states],cex=.5)
segments(x0=x[-length(x)],y0=y[-length(y)],x1=x[-1],y1=y[-1],
         col=c("#E69F00", "#56B4E9")[turtleFits$miSum$Par$states][-length(turtleFits$miSum$Par$states)],lwd=1.3)
#points(coordinates(turtleData)/1000,cex=0.15)
for(i in 1:length(turtleFits$miSum$errorEllipse))
  lines(turtleFits$miSum$errorEllipse[[i]]/1000,col=adjustcolor(c("#E69F00", "#56B4E9")[turtleFits$miSum$Par$states][i],alpha.f=0.25),cex=0.6)
dev.off()
rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])

###################################################
### Grey seal example
###################################################
#source(paste0(getwd(),"/greySealExample_TPM.R"))
load(paste0(example_wd,"greySealResults_TPM.RData"))

greySealTimeInStates<-greySealPool$Par$timeInStates

append.RData(greySealTimeInStates,file=paste0(getwd(),"/vignette_inputs.RData"))

load(paste0(getwd(),"/greySealData_TPM.RData"))
load(paste0(getwd(),"/coastUTM.RData"))

colvect<-c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2")

png(file=paste0(getwd(),"/plot_greySealResults1.png"),width=7.25,height=5,units="in",res=90)
plot(greySealPool$data$x,greySealPool$data$y,type="o", asp=1, pch=20,ylab="Latitude",xlab="Longitude",ylim=c(6060000,6289000),xlim=c(513000,925000),xaxt="n",yaxt="n",xaxs="i",yaxs="i")
axis(1,labels=c(-2,0,2,4),at=c(565312,696098,(434680-312928+696098),(565330-312928+696098)))
axis(2,labels=c(55,56,57),at=c(6094820,6206100,6317500))

points(greySealData$x,greySealData$y,pch=20)
cx=1.6

for(j in 1:5) {
  points(greySealPool$data$x[greySealPool$Par$states==j],greySealPool$data$y[greySealPool$Par$states==j],col=colvect[j],asp=1, pch=20,cex=cx)
}
points(greySealPool$data$x[length(greySealPool$data$x)],greySealPool$data$y[length(greySealPool$data$y)],col=colvect[greySealPool$Par$states[length(greySealPool$Par$states)]],asp=1, pch=20,cex=cx)

for(i in 1:(length(greySealPool$data$x)-1)) {
  polygon(greySealPool$errorEllipse[[i]],lty=2,border=rgb(red=col2rgb(colvect,alpha=FALSE)[1,greySealPool$Par$states[i]]/255,green=col2rgb(colvect,alpha=FALSE)[2,greySealPool$Par$states[i]]/255,blue=col2rgb(colvect,alpha=FALSE)[3,greySealPool$Par$states[i]]/255,alpha=0.25))
}
polygon(greySealPool$errorEllipse[[length(greySealPool$errorEllipse)]],lty=2,border=rgb(red=col2rgb(colvect,alpha=FALSE)[1,greySealPool$Par$states[length(greySealPool$Par$states)]]/255,green=col2rgb(colvect,alpha=FALSE)[2,greySealPool$Par$states[length(greySealPool$Par$states)]]/255,blue=col2rgb(colvect,alpha=FALSE)[3,greySealPool$Par$states[length(greySealPool$Par$states)]]/255,alpha=0.25))

#for(i in 1:nrow(greySealPool$Par$stateProbs$est)){
#  tmp<-sort(greySealPool$Par$stateProbs$est[i,])
#  for(j in 1:(5-1)){
#    if(tmp[j]>0.05){
#      ind<-order(greySealPool$Par$stateProbs$est[2,])[j]
#      points(greySealPool$data$x[i],greySealPool$data$y[i],col=colvect[ind],asp=1, pch=20,cex=cx/(j+1))
#    }
#  }
#}

leg=character(5)
leg[1]="Abertay haul-out state"
leg[2]="Farne Islands haul-out state"
leg[3]="Dogger Bank foraging state"
leg[4]="Low-speed exploratory state"
leg[5]="High-speed exploratory state"

legend("topright",leg,pch=rep(20,length(leg)),col=colvect,bty="n")

lines(coastUTM$Easting, coastUTM$Northing)

#add estimated locations for centres of attraction
points(centers,col="#D55E00",pch=20,cex=cx)
dev.off()

png(file=paste0(getwd(),"/plot_greySealResults2.png"),width=7.25,height=5,units="in",res=90)
plot(greySealSim$x,greySealSim$y,type="o", asp=1, pch=20,ylab="Latitude",xlab="Longitude",ylim=c(6060000,6289000),xlim=c(513000,925000),xaxt="n",yaxt="n",xaxs="i",yaxs="i")
points(greySealSim$x,greySealSim$y,col=colvect[greySealSim$states],pch=20)
axis(1,labels=c(-2,0,2,4),at=c(565312,696098,(434680-312928+696098),(565330-312928+696098)))
axis(2,labels=c(55,56,57),at=c(6094820,6206100,6317500))
legend("topright",leg,pch=rep(20,length(leg)),col=colvect,bty="n")
lines(coastUTM$Easting, coastUTM$Northing)
points(centers,col="#D55E00",pch=20,cex=1.6)
dev.off()
rm(list=ls()[-which(ls()=="example_wd" | ls()=="append.RData")])