#This is a simplified version of the workflow used in
#Bedrinana-Romano, L., Hucke-Gaete, R., Viddi, F.A., Johnson,
#D., Zerbini, A.N., Morales, J., Mate, B., Palacios, D.M., 2021.
#Defining priority areas for blue whale conservation and
#investigating overlap with vessel traffic in Chilean Patagonia,
#using a fast-fitting movement model. Scientific Reports 11, 2709.
#The model presented here is a continuous-time correlated random-walk model,
#originally developed by Devin Johnson
#(Johnson, D.S., London, J.M., Lea, M.-A., Durban, J.W., 2008.
#Continuous-time correlated random walk model for animal telemetry data.
#Ecology 89, 1208-1215.) and modified to allow the incorporation of covariates
#modulating the main movement parameters
###################Load required packages##################
library(argosfilter)
library(adehabitatLT)
library(rgeos)
library(rerddapXtracto)
library(raster)
library(rgdal)
library(TMB)
library(TMBhelper)
library(ggplot2)
library(mapdata)
library(ggpubr)
library(lubridate)
library(rmarkdown)
####################Load and preprocess ARGOS telemetry data#################
##You might need to change the names of the columns in accordance with your files
#######################################################################
#We are going to use telemetry data from a blue whale tagged on southern Chile
#see Hucke-Gaete, R., Bedrinana-Romano, L., Viddi, F.A., Ruiz, J.E., Torres-Florez, J.P.
#, Zerbini, A.N., 2018. From Chilean Patagonia to Galapagos, Ecuador:
#novel insights on blue whale migratory pathways along the Eastern South Pacific.
#PeerJ 6, e4695. for more details
data<-read.csv("84494-Locations.csv")
data$Date <- as.POSIXct(data$Date,"%H:%M:%S %d-%m-%Y",tz="GMT")
data=data[complete.cases(data$Error.Semi.major.axis),]
data$id=as.factor(data$Ptt)
#Set projection
xy=data.frame(x=data$Longitude,y=data$Latitude)
coordinates(xy)=c("x","y")
proj4string(xy) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" )
############################Plot raw data#################################################
#set limits of the map
xlimmap <- c(min(data$Longitude)-2, max(data$Longitude)+2)
ylimmap <- c(min(data$Latitude)-2, max(data$Latitude)+2)
# get outline data for map
w <- map_data("world", ylim = ylimmap, xlim = xlimmap)
# plot using ggplot
z <- ggplot(data,aes(x=Longitude, y=Latitude)) +
geom_point(aes(colour="red"), size=2.) +
scale_shape_manual(values=c(19,1))
z_1<-z + geom_polygon(data = w, aes(x=long, y = lat, group = group), fill = "black") +
theme_bw() + theme(legend.position = "none")+
coord_fixed(1.3, xlim = xlimmap, ylim = ylimmap) + ggtitle("Blue Whale #174070")
ggarrange( z_1,ncol = 1, nrow = 1)
##############################################################################
#Change to UTM
xy=spTransform(xy,CRS("+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
xy=data.frame(xy)
data$E=xy$x
data$N=xy$y
#Filter extreme locations
filtered.data<-vmask(lat=data$Latitude, lon=data$Longitude, dtime=data$Date, vmax=5)
data$filter<-filtered.data
data<-data[!data$filter=="removed",]
data<-data[!duplicated(data$Date),]
#Calculate deltas, which are time differences between locations
data_2 <- as.ltraj(xy = data[,c("Longitude","Latitude")], date = data$Date, id = data$id)
data_2<-ld(data_2)
data$delta<-data_2$dt
rm(data_2)
#Calculate variances for observational error
#Variances for modelling error in locations are derived from the
#Argos error ellipse and calculated as indicated in McClintock et al. (2015)
#This only works for tags deployed from 2013 onwards. See the paper for a solution on tags deployed earlier
data$v1 = (data$Error.Semi.major.axis/sqrt(2))^2*sin(pi*(data$Error.Ellipse.orientation/180))^2 + (data$Error.Semi.minor.axis/sqrt(2))^2*cos(pi*(data$Error.Ellipse.orientation/180))^2
data$v2 = (data$Error.Semi.major.axis/sqrt(2))^2*cos(pi*(data$Error.Ellipse.orientation/180))^2 + (data$Error.Semi.minor.axis/sqrt(2))^2*sin(pi*(data$Error.Ellipse.orientation/180))^2
#Download covariate data using rerddapXtracto
#Check out this cool tutorial on how to use rerddapXtracto
#https://coastwatch.pfeg.noaa.gov/xtracto/
chlInfo <- rerddap::info("erdMH1chlamday")#Here we will use Chlorophyll-a as covariate
#######I commented the following line of code to speed up rmarkdown process
#The function rerddapXtracto::rxtracto accesses erddap database online for downloading the requiered data
#chl<- rerddapXtracto::rxtracto(chlInfo, parameter = 'chlorophyll', xcoord = data$Longitude, ycoord = data$Latitude, tcoord = data$Date, zcoord = NULL, xlen = .01, ylen = .01)
mu.chl.cov=chl$`mean chlorophyll`
any(is.na(mu.chl.cov))#Check for NAs
## [1] TRUE
#NAs are extremely common with chlorophyll data
#just to get this going, lets fill NAs, with the values
#immediately before
for(i in 2:length(mu.chl.cov)){
if(is.na(mu.chl.cov[i])){
mu.chl.cov[i]=mu.chl.cov[i-1]
}
}
n0=length(data$E)
data$log.chl=log(mu.chl.cov)#Transform Chl data into log scale
Ialpha1=rbind(data$E[1],0)#Matrix with mean values for first speed and location in X
Ialpha2=rbind(data$N[1],0)#Matrix with mean values for first speed and location in Y
data=data[4:nrow(data),]#Crop the first three observations because no covariate values were available for them
#Load data for TMB
data.TMB<-list(x=data$E,
y=data$N,
Ialpha1=Ialpha1,#This matrix holds the first observed location and speed in x
Ialpha2=Ialpha2,#This matrix holds the first observed location and speed in y
delta=(data$delta[1:(n0-1)])/3600,
logchl=as.vector(scale(data$log.chl)),
errx=sqrt(data$v1),
erry=sqrt(data$v2),
sd1=1000,sd2=1000,
sd_beta=0.001,sd_sigma=0.001)
n=length(data.TMB$x)
#Initial values for state variables
u1=data.TMB$x
u2=data.TMB$y
V1=rep(0,length(data.TMB$y))
V2=rep(0,length(data.TMB$y))
log_beta=rep(0,length(data.TMB$y)-1)
log_sigma=rep(9,length(data.TMB$y)-1)
alpha1=rbind(u1,V1)
alpha2=rbind(u2,V2)
#########################Fit model###################################
#######I commented the following lines of code to speed up rmarkdown process
####But they need to ran to fit the model
# TMB::compile("CTCRW_matrix_cov.cpp","-O1 -g",DLLFLAGS="")
# dyn.load(dynlib("CTCRW_matrix_cov"))
# set.seed(123)
# parameters <- list(alpha1=alpha1, alpha2=alpha2,log_sigma=log_sigma,log_beta=log_beta,B0=0,B1=0,A0=9,A1=0)#
# obj <- MakeADFun(data.TMB, parameters, random=c("alpha1","alpha2","log_sigma","log_beta"),DLL="CTCRW_matrix_cov")
# obj$hessian <- T
# opt<-TMBhelper::fit_tmb(obj,getsd=T,newtonsteps=3)#lower=c(-2,-2,7,-2),upper = c(4,2,10,2)
#Lets check out our results
summary(opt$SD, "fixed", p.value = TRUE)
## Estimate Std. Error z value Pr(>|z^2|)
## A0 8.89535359 0.02383331 373.2319191 0.000000e+00
## A1 -0.35686063 0.02065009 -17.2813147 6.504829e-67
## B0 -0.01799578 0.07738357 -0.2325529 8.161086e-01
## B1 0.87203158 0.08766620 9.9471810 2.594279e-23
#As you can see there is a significant negative correlation between sigma and logChl
#and the opposite occurs with beta. This is, the whale reduces velocity
#and persistence as logChl increases
#Add estimated locations to our data.frame
est.alpha1<-matrix(c(obj$env$parList()$alpha1,obj$env$parList()$alpha2),ncol=length(data.TMB$x),nrow=2,byrow = F)
est.alpha2<-matrix(c(obj$env$parList()$alpha2,obj$env$parList()$alpha2),ncol=length(data.TMB$x),nrow=2,byrow = F)
data$est.X=est.alpha1[1,]
data$est.Y=est.alpha2[1,]
#Add long-term velocity, sigma and beta to our data.frame
data$nu=c(NA,((sqrt(pi)*exp(obj$env$parList()$log_sigma))/(2*sqrt(exp(obj$env$parList()$log_beta))))/1000)
data$sigma=c(NA,exp(obj$env$parList()$log_sigma))
data$beta=c(NA,exp(obj$env$parList()$log_beta))
#Behavioural variation can be addressed as a continuous variable nu (long term velocity)
#or it can be rescaled to range between 1 and 2 (Behavioural states)
#Here the lowest 50% of nu values are considered state 1 (ARS) and
#the rest of values are considered state 2 (transit)
data$beh=((data$nu-min(data$nu,na.rm = T))/(max(data$nu,na.rm = T)-min(data$nu,na.rm = T)))+1
data$behd=1
data$behd[data$beh>1.5]=2
###########################################Plots#################################################
#These 3 plots show the observed variation in logChl experienced by the blue whale
#and also two different ways of showing behaviroal variation
#As a continuous variable nu (long term velocity) or discretized into behavioural states post hoc
#Lets trasnform our estimted UTMs to Lat and Long
coordinates(data)=c("est.X","est.Y")
proj4string(data) <- CRS("+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
#Change to WGS84
data=spTransform(data,CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
data=data.frame(data)
# set limits of the map
xlimmap <- c(min(data$est.X)-2, max(data$est.X)+2)
ylimmap <- c(min(data$est.Y)-2, max(data$est.Y)+2)
# get outline data for map
w <- map_data("world", ylim = ylimmap, xlim = xlimmap)
# plot using ggplot
z <- ggplot(data,aes(x=est.X, y=est.Y)) +
geom_point(aes(colour=nu), size=1.5) +
scale_shape_manual(values=c(19,1))
z_1<-z + geom_polygon(data = w, aes(x=long, y = lat, group = group), fill = "black") +
theme_bw() +
scale_colour_gradient2(low = "red", mid = "purple",n.breaks=5,
high = "blue", midpoint = 14,limits=range(data$nu,na.rm = T), "nu km/h") +
coord_fixed(1.3, xlim = xlimmap, ylim = ylimmap) + ggtitle("Long-term velocity")
zb <- ggplot(data,aes(x=est.X, y=est.Y)) +
geom_point(aes(colour=log.chl), size=1.5) +
scale_shape_manual(values=c(19,1))
z_1b<-zb + geom_polygon(data = w, aes(x=long, y = lat, group = group), fill = "black") +
theme_bw() +
scale_colour_gradient2(low = "blue", mid = "green",n.breaks=5,
high = "darkgreen", midpoint = 1,limits=range(data$log.chl,na.rm = T), "log Chl") +
coord_fixed(1.3, xlim = xlimmap, ylim = ylimmap) + ggtitle("log Chl")
zc <- ggplot(data,aes(x=est.X, y=est.Y)) +
geom_point(aes(colour=behd), size=1.5) +
scale_shape_manual(values=c(19,1))
z_1c<-zc + geom_polygon(data = w, aes(x=long, y = lat, group = group), fill = "black") +
theme_bw() +
scale_colour_gradient2(low = "red", mid = "purple",n.breaks=5,
high = "blue", midpoint = 1.5,limits=range(data$beh,na.rm = T), "beh") +
coord_fixed(1.3, xlim = xlimmap, ylim = ylimmap) + ggtitle("Behavioral state")
ggarrange(z_1b,z_1,z_1c,ncol = 3, nrow = 1)
#Now lets animate this, check the original tutorial on this package
#here https://movevis.org/
library(moveVis)
library(move)
datamov=data
datamov=moveVis::df2move(df=datamov,proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",x="est.X",y="est.Y",time = "Date",track_id="id",data = datamov)
m <- moveVis::align_move(datamov, res = 24, unit = "hours")
## Warning in FUN(X[[i]], ...): CRS object has comment, which is lost in output
#######I commented the following lines of code to speed up rmarkdown process
# frames <- frames_spatial(m, path_legend=F,
# map_service = "mapbox", map_type = "satellite",map_token ="Here you need to put your personal token obtained at https://www.mapbox.com/") %>%
# add_labels(x = "Longitude", y = "Latitude") %>% # add some customizations, such as axis labels
# #add_northarrow() %>%
# #add_scalebar() %>%
# add_timestamps(m, type = "label") %>%
# add_progress()
#
#
#
# # animate frames
# animate_frames(frames, out_file = "moveVis.gif",fps = 10,overwrite=T,width = 700,height = 700,res=100)
# rm(frames)#I suggest removing the generated frames after the animation is done
#Too much information can make R unstable and produce failure in loading workspace after closing session
#rmarkdown::render("C:/Users/orcah_000/Desktop/W11/CTCRW-COVS.R")