#################################
###HOW TO
#CLASSIFY MOVEMENT MODES
#SPLIT TRAJECTORIES
#PERFORM STEP SELECTION FUNCTION
#MAKE PREDICTIONS
#MAKE SPATIAL PREDICTIONS
#ESTIMATE OVERLAP AND SPILLOVER
################################

###
#LIBRARIES
###
library(rgdal)
library(rgeos)
library(ggmap)
library(png)
library(adehabitatLT)
library(adehabitatHR)
library(animation)
library(lattice)
library(nls.multstart)
library(lubridate)
source("helper_updated.R")

#PLANAR PROJECTION 
proj="+proj=utm +zone=16 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

#LOAD SIMULATED MOOSE DATA
mo=read.csv("moose_track.csv")
mo$id=mo$GPMoose
mo=out.duplicated(mo,mo$id)
mo$year=year(dmy_hm(mo$datetime))
mo$idyear=paste(mo$id,mo$year)
mo=SpatialPointsDataFrame(mo[,c("Easting","Northing")],data=mo,proj4string=CRS(proj))
mo$datetime=dmy_hm(mo$datetime)

#LOAD SIMULATED DEER DATA
de=read.csv("deer_track.csv")
de$id=de$GPDeer
de$year=year(dmy_hm(de$datetime))
de$idyear=paste(de$id,de$year)
de=SpatialPointsDataFrame(de[,c("Easting","Northing")],data=de,proj4string=CRS(proj))
de$datetime=dmy_hm(de$datetime)

####
#RESUMING DATA
###
nd=table(de$GPDeer)
nd=data.frame(GPDeer=names(nd),n=as.vector(nd))
nm=table(mo$GPMoose)
nm=data.frame(GPMoose=names(nm),n=as.vector(nm))


###
#HERE WE EXEMPLIFY THE PROCCESS WITH DEER DATA. NOTE IN THE PAPER WE DID FOR DEER AND MOOSE
###

######
##BURN OUT FIRST LOCATION POST CAPTURE
#####
deBurn=burn.out(de,"id","datetime",cut=4)

#####
#SMOOTHING AND CLEANING TRAJECTORY
#####
deVel=clean.smoo(data=deBurn,id="id",date="datetime",quantVel=0.95,windowSize=1)

#ELIMINATE IDS WITH FEW LOCATIONS
deLong=eliminate(spdf=deVel,id=deVel$id,date="datetime",cut=180,units="days")

#RESAMPLING TRAJECTORY
library(trajr)
deLong=deLong[,c("id","idyear","Easting","Northing","datetime")]
t=lapply(split(deLong,deLong$idyear),function(x){TrajFromCoords(x@data,"Easting","Northing","datetime")})
#t=lapply(t,function(x){TrajResampleTime(x,60*60*24)})
t=lapply(t,function(x){resample.time(trj=x,stepTime=60*60*24*1,threshold=10)})

id=rep(names(t),unlist(lapply(t,nrow)))
t=do.call(rbind,t)
t$idyear=id
t=SpatialPointsDataFrame(data.frame(t[,c("x","y")]),data=data.frame(t))
t$year=strptime(t$time,format="%Y-%m-%d %H:%M:%S")$year+1900
t$id=unlist(lapply(strsplit(t$idyear," "),function(x){x[1]}))

#CALCULATE MSD AND ASSOCIATED TIME
tclean=prep_timeMig(spdf=t,id=t$id,join.idyear=t$idyear,date=t$time) #TIME BY YEAR ID
tclean=msd(coords=tclean,id="idyear",date="time",units="days")
plot.msd(from.msd=tclean,id="idyear",type="msd_time",same=F)
declean=tclean

#RUN NOMADIC MODEL
nom=nomadic.model(data=declean)
nom$coefs
nom$models

#RUN SEDENTARY MODEL
sed=sedentary.model(data=declean)
sed$coefs
sed$models

#RUN DISPERSAL MODEL
disp=dispersal.model(data=declean)
disp$coefs
disp$models

#RUN MIGRATION MODEL
mig=migration.model(data=declean)
mig$coefs
mig$models

######
##COMPARE MODELS
######
resume=data.frame(
nomadic=nom$coefs$AIC,
sedentary=sed$coefs$AIC-2,
disperser=disp$coefs$AIC,
migrater=mig$coefs$AIC)
rownames(resume)=names(nom$model)
resume

resume$mode=modeAIC(resume,delta=7)
table(resume$mode)
table(resume$mode)/sum(table(resume$mode))

plot.mig(data=declean,id="idyear",migModel=mig$models,mode=resume$mode)
plot.mig(data=declean,id="idyear",migModel=sed$models,mode=resume$mode)
plot.mig(data=declean,id="idyear",migModel=disp$models,mode=resume$mode)

##################
##EXTRACT METRICS#
##################
mig_metrics=mig.metrics(mig$coef,resume$mode)
disp_metrics=disp.metrics(disp$coef,resume$mode)
kernel.mig(mig_metrics,disp_metrics,h=20)

######################
###SPLIT TRAJECTORIES#
######################
#EXACT DATES
disp_dates=find.date(disp_metrics,type="disperser")
mig_dates=find.date(mig_metrics,type="migrater")

#disperser
dc=classify.track(data=declean,id=declean$idyear,date="time",find.dates=disp_dates,type="disperser")
d=classify.track(data=de,id=de$idyear,date="datetime",find.dates=disp_dates,type="disperser")
plot.msd(from.msd=dc,id="idyear",type="msd_time",same=F)

#migrater
mc=classify.track(data=declean,declean$idyear,"time",mig_dates,type="migrater")
m=classify.track(data=de,de$idyear,"datetime",mig_dates,type="migrater")
plot.msd(from.msd=mc,id="idyear",type="msd_time",same=F)

#sedentary
s=de[de$idyear %in% rownames(resume)[resume$mode=="sedentary"],]
s$mode=ifelse(month(s$datetime) %in%  c(4,11),"dispersal1",
       ifelse(month(s$datetime) %in% c(12,1:3),"winter1","summer"))


#########
###SSF###
#########
deC=rbind(d,m)
#deC=s
deC$mode1=deC$mode
deC$mode1[deC$mode1 %in% c("winter1","winter2")]="w"
deC$mode1[deC$mode1 %in% c("dispersal1","dispersal2")]="d"
deC$mode1[deC$mode1 %in% c("summer")]="s"
plot(deC,col=factor(deC$mode1),cex=.1,axes=T)

library(amt)
ldec=split(deC,deC$id)
ltraj=lapply(ldec,make_track,Easting,Northing,datetime,id=id,idyear=idyear,mode=mode1)
freq=do.call(rbind,lapply(ltraj,summarize_sampling_rate))
hist(freq$median)
lred=lapply(ltraj,track_resample,rate=hours(4),tolerance = minutes(30))
lstep=lapply(lred,steps_by_burst,keep_cols="start")
lrand=lapply(lstep,random_steps,n=30)

########
#LAYERS#
########
map=raster("habitat.tif")
elev=raster("elevation.tif")
wolf=raster("wolf_pressure.tif")
lext=lapply(lrand,extract_covariates,map,where="end") 
lext=lapply(lext,extract_covariates,elev,where="end") 
lext=lapply(lext,extract_covariates,wolf,where="end") 
ssf=data.frame(do.call(rbind,lext))
ssf$strata=paste(ssf$id,ssf$step_id_)
hab=data.frame(index=1:8,hab=c("wat","urb","open","df","cf","mf","wet","crop"))
ssf$hab=hab$hab[match(ssf$habitat,hab$index)]
ssf$hab=relevel(factor(ssf$hab),"cf")
ssf$mode=relevel(factor(ssf$mode),"w")
ssf$slayer=scale(ssf$elevation)
ssf$wolf_pressure=scale(ssf$wolf_pressure)
#save(ssf,file="DeerSsfData.RData")

###
#FITTING SSF MODEL USING TEMPLATE MODEL BUILDER
###
library(glmmTMB)
library(fastDummies)
ssf1=dummy_cols(ssf,select_columns="hab",ignore_na=T)
ssf1=dummy_cols(ssf1,select_columns="mode",ignore_na=T)


######
##JUST HABITAT
######
t=Sys.time()
hab.fit= glmmTMB(case_ ~ -1+ 
                slayer+
                hab_mf+ hab_df+ hab_open+ hab_urb+ hab_wet+hab_wat+ #hab_crop+ 
                (1|id/strata),family=poisson, data=na.omit(ssf1),
                start=list(theta=rep(log(1e3),2)),
                map=list(theta=factor(c(NA,NA))))
Sys.time()-t
summary(hab.fit)

######
##MODE MODULATION OF HABITAT
######
t=Sys.time()
mode.fit=glmmTMB(case_ ~ -1+ 
                #FIXED EFFECTS
                slayer+
                slayer:mode_d+slayer:mode_s+
                hab_mf+ hab_df+ hab_open+ hab_urb+ hab_wet+hab_wat+  #hab_crop+ 
                hab_mf:mode_d+ hab_df:mode_d+ hab_open:mode_d+ hab_urb:mode_d+ hab_wet:mode_d+hab_wat:mode_d+ #hab_crop:mode_d+
                hab_mf:mode_s+ hab_df:mode_s+ hab_open:mode_s+ hab_urb:mode_s+ hab_wet:mode_s+hab_wat:mode_s+ #hab_crop:mode_s+ 
                #RANDON EFFECTS
                (1|id/strata),family=poisson, data=na.omit(ssf1),
                #FIX VARIANCE OF RANDON EFFECTS (FOR STRATA AND ID)
                start=list(theta=rep(log(1e3),2)),
                map=list(theta=factor(c(NA,NA))))
Sys.time()-t
summary(mode.fit)

####################
##WOLF MODULATION OF HABITAT
#####################
t=Sys.time()
wolf.fit= glmmTMB(case_ ~ -1+ 
                #FIXED EFFECTS 
                hab_mf+ hab_df+ hab_open+ hab_urb+ hab_wet+hab_wat+  #hab_crop+
                slayer+
                wolf_pressure+
                #WOLF PRESSURE INTERACTION
                slayer:wolf_pressure+
                hab_mf:wolf_pressure+ hab_df:wolf_pressure+ hab_open:wolf_pressure+ hab_urb:wolf_pressure+ hab_wet:wolf_pressure+hab_wat:wolf_pressure+  #hab_crop:wolf_pressure+
                #RANDON EFFECTS
                (1|id/strata),family=poisson, data=na.omit(ssf1),
                #FIX VARIANCE OF RANDON EFFECTS (FOR STRATA AND ID)
                start=list(theta=rep(log(1e3),2)),
                map=list(theta=factor(c(NA,NA))))
Sys.time()-t
summary(wolf.fit)

##########
##FULL MODE AND WOLF MODULATION
##########
t=Sys.time()
full.fit= glmmTMB(case_ ~ -1+ 
                 #FIXED EFFECTS
                 slayer+
                 hab_mf+ hab_df+ hab_open+ hab_urb+ hab_wet+hab_wat+  #hab_crop+
                 wolf_pressure+
                 #MODE INTERACTION
                 slayer:mode_s+slayer:mode_d+
                 hab_mf:mode_d+ hab_df:mode_d+ hab_open:mode_d+ hab_urb:mode_d+ hab_wet:mode_d+hab_wat:mode_d+  #hab_crop:mode_d+ 
                 hab_mf:mode_s+ hab_df:mode_s+ hab_open:mode_s+ hab_urb:mode_s+ hab_wet:mode_s+hab_wat:mode_s+ #hab_crop:mode_s+
                 #wolf_pressure:mode_s+wolf_pressure:mode_d+
                 #WOLF PRESSURE INTERACTION
                 slayer:wolf_pressure+
                 hab_mf:wolf_pressure+ hab_df:wolf_pressure+ hab_open:wolf_pressure+ hab_urb:wolf_pressure+ hab_wet:wolf_pressure+hab_wat:wolf_pressure+ #hab_crop:wolf_pressure+ 
                 #RANDON EFFECTS - INTERCEPTS ONLY
                 (1|id/strata),family=poisson, data=na.omit(ssf1),
                 #FIX VARIANCE OF RANDON EFFECTS (FOR STRATA AND ID)
                 start=list(theta=rep(log(1e3),2)),
                 map=list(theta=factor(c(NA,NA))))
Sys.time()-t
summary(full.fit)

###
#MODEL RANKING
###
AIC(hab.fit)
AIC(mode.fit)
AIC(wolf.fit)
AIC(full.fit)
library(AICcmodavg)
aictab(list(hab=hab.fit,mode=mode.fit,wolf=wolf.fit,full=full.fit))
ssf.deer=list(hab.fit,mode.fit,wolf.fit,full.fit)


###########################################################
##HERE FINISHED THE DATA PREPARATION TO RUN SSF FOR DEER
##REPETE THE SAME PROCEDURE FOR MOOSE
############################################################

#LOAD HERE THE RESULTS FOR MOOSE - IT IS A SHORTCUT IF YOU DO NOT WANT TO RUN EVERITHING ABOVE AGAIN
load("ssf_moose.RData")


###
##PREDICTION FOR DEER AND MOOSE HABITAT SELECTION
###
#PLOT RESULTS FOR OBSERVED WOLF PREDATION SCENARIO
pw=plot.ssfMode(ssf.deer[[4]],ssf.moose[[4]],wp="actual")
#PLOT RESULTS FOR MINIMUM WOLF PREDATION SCENARIO
pnw=plot.ssfMode(ssf.deer[[4]],ssf.moose[[4]],wp="min")

##
#SPATIAL PREDICTION
##
load("Elevation.RData")
wolf=raster("wolf_pressure.tif")
map=raster("map_USA-CAN.tif")


map=raster("habitat.tif")
elev=raster("elevation.tif")
wolf=raster("wolf_pressure.tif")

#ADJUSTS RASTERS EXTENSION
elev=resample(elev,map,method="ngb")
wolf=resample(wolf,map,method="ngb")

#MASK GEOGRAPHIC ARAE FOR PREDICTION
join=rbind(de[,"id"],mo[,"id"])
mcp=mcp(join[,-1],percent=100)
mask=gBuffer(mcp,width=10^4)
elev=mask(elev,mask)
map=mask(map,mask)
wolf=mask(wolf,mask)

#MASK
load("Deer.RData")
load("Moose.Rdata")
co=rbind(de[,"Month"],mo[,"Month"])
mask=mcp(co[,-1],100)
mask=gBuffer(mask,width=10^4)
elev=mask(elev,mask)
map=mask(map,mask)
wolf=mask(wolf,mask)

maps=stack(map,elev,wolf)
names(maps)=c("hab","slayer","wolf_pressure")
hab=data.frame(index=1:8,hab=c("wat","urb","open","df","cf","mf","wet","crop"))

#################################
##SPATIAL PREDICTION FOR DEER####
#################################
##DEER FOR OBSERVED SCENARIO
#WINTER
dumhab=dummy.raster(maps,hab,season="w",scale=T,predator="actual")
t=Sys.time()
pw=predict(dumhab,ssf.deer[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pw))
#DISPERSAL
dumhab=dummy.raster(maps,hab,season="d",scale=T,predator="actual")
t=Sys.time()
pd=predict(dumhab,ssf.deer[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pd))
#SUMMER
dumhab=dummy.raster(maps,hab,season="s",scale=T,predator="actual")
t=Sys.time()
ps=predict(dumhab,ssf.deer[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(ps))

map.deerW=stack(rastToIso(pw),rastToIso(pd),rastToIso(ps))
names(map.deerW)=c("dw","dd","ds")

##DEER FOR MIMINUM SCENARIO
#WINTER
dumhab=dummy.raster(maps,hab,season="w",scale=T,predator="min")
t=Sys.time()
pw=predict(dumhab,ssf.deer[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pw))
#DISPERSAL
dumhab=dummy.raster(maps,hab,season="d",scale=T,predator="min")
t=Sys.time()
pd=predict(dumhab,ssf.deer[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pd))
#SUMMER
dumhab=dummy.raster(maps,hab,season="s",scale=T,predator="min")
t=Sys.time()
ps=predict(dumhab,ssf.deer[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(ps))

map.deerNW=stack(rastToIso(pw),rastToIso(pd),rastToIso(ps))
names(map.deerNW)=c("dw","dd","ds")


##########################
##SPATIAL PREDICTION FOR MOOSE####
##########################
##MOOSE FOR OBSERVED SCENARIO
#WINTER
dumhab=dummy.raster(maps,hab,season="w",scale=T,predator="actual")
t=Sys.time()
pw=predict(dumhab,ssf.moose[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pw))
#DISPERSAL
dumhab=dummy.raster(maps,hab,season="d",scale=T,predator="actual")
t=Sys.time()
pd=predict(dumhab,ssf.moose[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pd))
#SUMMER
dumhab=dummy.raster(maps,hab,season="s",scale=T,predator="actual")
t=Sys.time()
ps=predict(dumhab,ssf.moose[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(ps))
map.mooseW=stack(rastToIso(pw),rastToIso(pd),rastToIso(ps))
names(map.mooseW)=c("mw","md","ms")

##MOOSE FOR MIMINUM SCENARIO
#WINTER
dumhab=dummy.raster(maps,hab,season="w",scale=T,predator="min")
t=Sys.time()
pw=predict(dumhab,ssf.moose[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pw))
#DISPERSAL
dumhab=dummy.raster(maps,hab,season="d",scale=T,predator="min")
t=Sys.time()
pd=predict(dumhab,ssf.moose[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(pd))
#SUMMER
dumhab=dummy.raster(maps,hab,season="s",scale=T,predator="min")
t=Sys.time()
ps=predict(dumhab,ssf.moose[[4]],na.rm=F,type="response")
Sys.time()-t
plot(1-rastToIso(ps))
map.mooseNW=stack(rastToIso(pw),rastToIso(pd),rastToIso(ps))
names(map.mooseNW)=c("mw","md","ms")


###
#PLOTTING OVERLAP IN OBSERVED AND MINIMUM SCENARIO
##
require(raster)
md=(1-map.deerW)>.4
plot(md)
mm=(1-map.mooseW)>.4
plot(mm)
mm[mm==1]=2
plot(md+mm)
overW=md+mm

md=(1-map.deerNW)>.4
plot(md)
mm=(1-map.mooseNW)>.4
plot(mm)
mm[mm==1]=2
plot(md+mm)
overNW=md+mm

layout(matrix(1:6,3,2),heights=c(1.3,1,1.3),widths=c(1,.9))
plot.space(overW,col.wolf="gray30",col=
             c(NA,adjustcolor("mediumpurple4",.9),adjustcolor("steelblue4",.9),"orange"),par=F,axis=T,wolf=T,sea=T)
plot.space(overNW,col.wolf="gray30",col=
             c(NA,adjustcolor("mediumpurple4",.9),adjustcolor("steelblue4",.9),"orange"),par=F,axis=F,wolf=F,sea=F)