########################
#AUXILIARY FUNCTIONS####
########################

############
#GET GMAPS##
############
goomap=function(center,zoom=15,nosignal=0,crs=CRS("+proj=utm"),adjlat=0,adjlon=0){
register_google(key = "AIzaSyCmuIHFjiFb0-11Aht2Ko2Ue2cuO_Ol6rA")
#AIzaSyABdKk1XuC6u9CYWy7P8GBz-q-Kq57IHnA
nhumi=get_map(location=coordinates(center)+c(adjlon,adjlat),maptype="satellite",zoom =zoom)
png("nhumi.png")
par(mar=c(0,0,0,0))
plot(nhumi)
dev.off()
n=readPNG("nhumi.png")
bb=data.frame(attr(nhumi,"bb"))
names(bb)=rep(c("a","b"),2)
par(mar=c(4,4,0,1))
bbox=SpatialPoints(rbind(bb[2:1],bb[4:3]),proj4string=CRS("+datum=WGS84 +proj=longlat"))
bbox=data.frame(spTransform(bbox,CRS=crs))
plot(bbox[,1],bbox[,2]+c(50,0),type="n",xlab="Longitude",ylab="Latitude",cex.axis=0.8,cex.main=0.8)
rasterImage(n,bbox[1,1],bbox[1,2],bbox[2,1],bbox[2,2])
rect(bbox[1,1],bbox[1,2]+1,bbox[2,1],bbox[1,2]+nosignal,col="white",border="white")
#north.arrow(749400,7730850,len=15,lab="N",col=1)
#GISTools::map.scale(x=748150, y=7729650, len=200,ndivs=2,units="100m") 
}

####################
#PLOT TRAJECTORIES##
####################
plot.traj=function(coords,id,col=NULL,add=T,type="p",cex=.5){
n=length(unique(id))
if(is.null(col)){col=sample(colours(),n)}
else{col=col}
names(col)=unique(id)
coords$col=col[match(id,names(col))]
plot(coords,col=coords$col,pch=16,cex=cex,add=add)
l=split(coords,id)
if(type=="l"){invisible(lapply(l,function(x){lines(coordinates(x),col=x$col,lty=3)}))}
}

########################
#TRACKING THROUGH TIME##
########################
tracktime=function(list){
list=list[order(as.numeric(names(list)))]
#list=list[order(names(list))]
plot(range(do.call("c",list),na.rm=T),c(1,length(list)),type="n",ylab="",xlab="Dates",yaxt="n")
mtext("IDS",2)
for(i in 1:length(list)){
lines(c(list[[i]][1],list[[i]][2]),c(i,i),lwd=2)
text(min(c(list[[i]][1],list[[i]][2])),i+.3,names(list)[i],cex=.6)
}
}

############################
#TRACKING THROUGH TIME DOT##
########################
tracktimedot=function(list){
list=list[order(as.numeric(names(list)))]
plot(range(do.call(c,list),na.rm=T),c(1,length(list)),type="n",ylab="",xlab="Dates",yaxt="n")
mtext("IDS",2)
for(i in 1:length(list)){
lines(c(list[[i]][1],list[[i]][2]),c(i,i),lwd=2)
text(min(c(list[[i]][1],list[[i]][2])),i+.3,names(list)[i],cex=.6)
}
}

###################################
##MEASURING MOVED DISTANCE BY HOUR##
###################################
act=function(coords,id,datetime){
traj=as.ltraj(coordinates(coords),id=id,date=datetime)
coords$dist=ld(traj)$dist
coords$h=strptime(trunc(datetime,"hour"),format="%Y-%m-%d %H:%M:%S")$h
coords$day=trunc(datetime,"days")
s=split(coords,id)
ss=lapply(s,function(x){tapply(x$dist,list(x$day,hour=x$h),sum,na.rm=T)})
sss=lapply(ss,function(x){apply(x,2,median,na.rm=T)})
s4=do.call(rbind,sss)
#plot(c(0,23),range(s4),type="n")
plot(c(0,23),c(0,1),type="n",ylab="Activity level",xlab="TIme of day",xaxp=c(0,24,12))
apply(s4,1,function(x){lines(0:23,x/max(x),,lty=3)})
lines(0:23,apply(s4,2,median)/max(apply(s4,2,median)),lty=1,col="red",lwd=10)
return(s4)
}

############################
#MEAN SQUARED DISPLACEMENT##
############################
msd=function(coords,id="id",date="datetime",units="days"){
datetime=coords@data[,date]
id=coords@data[,id]
coords=coords[order(id,datetime),]
id=id[order(id,datetime)]
datetime=datetime[order(id,datetime)]
s=split(coords,factor(id))
#s=lapply(s,function(x){x[order(x$datetime),]})
date=lapply(s,function(x){difftime(x@data[,date],x@data[,date][1],units=units)})
s=lapply(s,function(x){gDistance(x[1,],x,byid=T)})
#s=data.frame(id,datetime,time=as.vector(do.call("c",date)),msd=as.vector(do.call(rbind,s)))
#return(s)
coords@data=cbind(coords@data,time0=as.vector(do.call("c",date)),msd=as.vector(do.call(rbind,s)))
return(coords)
}

####################################
#PLOTING MEAN SQUARED DISPLACEMENT##
####################################
plot.msd=function(from.msd,id,type="time0",same=T){
from.msd$t=from.msd@data[,type]
s=split(from.msd,from.msd@data[,id])
n=sort(unique(factor(as.vector(from.msd@data[,id]))))
par(mfrow=c(ceiling(sqrt(length(n))),floor(sqrt(length(n)))+1),mar=c(0,0,0,0))
#par(mar=c(4,4,0,1))
for(i in 1:length(n)){
temp=from.msd[from.msd@data[,id]==as.vector(n[i]),]
if(same==T){plot(range(temp$t),range(from.msd$msd),type="n",yaxt="n",xaxt="n",xlab="",ylab="")}
else{plot(range(from.msd$t),range(from.msd$msd),type="n",yaxt="n",xaxt="n",xlab="",ylab="")}
lines(temp$t,temp$msd)
mtext(n[i],3,-1,cex=.7)
}
}

#######################
#TRAJECTORY CLEANING###
#######################
out.duplicated=function(coords,id){
s=split(coords,id)
s=lapply(s,function(x){x[!duplicated(x$datetime),]})
return(do.call(rbind,s))
}

####################################
#PLOTING MEAN SQUARED DISPLACEMENT##
####################################
plot.msd1=function(from.msd,type="time",same=T,who=1){
if(type=="time"){from.msd$t=from.msd$time}
if(type=="datetime"){from.msd$t=from.msd$datetime}
s=split(from.msd,from.msd$id)
n=sort(unique(as.numeric(as.vector(from.msd$id))))
#par(mfrow=c(ceiling(sqrt(length(n))),floor(sqrt(length(n)))),mar=c(0,0,0,0))
m=c(rep(1,5),2,3,4,5,rep(1,5),6,7,8,9,rep(1,5),10,11,12,13,rep(1,5),14,15,16,17,rep(1,5),18:66)
#windows(80,80)
layout(matrix(m,9,10))
par(mar=c(4,4,0,1))
temp=from.msd[from.msd$id==who,]
if(same==T){plot(range(temp$t),range(from.msd$msd),type="n",ylab="Mean Squared Displacement",xlab="Date")}
else{plot(range(from.msd$t),range(from.msd$msd),type="n",ylab="Mean Squared Displacement",xlab="Date")}
lines(temp$t,temp$msd)
mtext(who,3,-1,cex=.7)
n=n[-which(n==who)]
for(i in 1:length(n)){
temp=from.msd[from.msd$id==n[i],]
par(mar=c(0,0,0,0))
if(same==T){plot(range(temp$t),range(from.msd$msd),type="n",yaxt="n",xaxt="n",xlab="",ylab="")}
else{plot(range(from.msd$t),range(from.msd$msd),type="n",yaxt="n",xaxt="n",xlab="",ylab="")}
lines(temp$t,temp$msd)
mtext(n[i],3,-1,cex=.7)
}
}

########################
###GIFS ANIMATION#######
########################
anima=function(x,ids=NULL,idn="GPDeer",zoom=10,crs=CRS(proj4string(x)),adjlat=0,adjlon=0,nosignal=2*10^3,delay=6,jump=1,col=NULL){
cm=SpatialPoints(matrix(colMeans(coordinates(x)),1,2),proj4string=CRS("+proj=utm +zone=16 +north"))
c=spTransform(cm,CRS=CRS("+proj=longlat"))
register_google(key = "AIzaSyCmuIHFjiFb0-11Aht2Ko2Ue2cuO_Ol6rA")
#AIzaSyABdKk1XuC6u9CYWy7P8GBz-q-Kq57IHnA
nhumi=get_map(location=coordinates(c)+c(adjlon,adjlat),maptype="satellite",zoom =zoom)
png("nhumi.png")
par(mar=c(0,0,0,0))
plot(nhumi)
dev.off()
n=readPNG("nhumi.png")
bb=data.frame(attr(nhumi,"bb"))
names(bb)=rep(c("a","b"),2)
par(mar=c(4,4,0,1))
bbox=SpatialPoints(rbind(bb[2:1],bb[4:3]),proj4string=CRS("+datum=WGS84 +proj=longlat"))
bbox=data.frame(spTransform(bbox,CRS=crs))
if(!is.null(ids)){x=x[x@data[,idn] %in% ids,]}
ids=unique(x@data[,idn])

jump=jump
trans=.1
delay=delay
start=delay+1
if(is.null(col)){col=sample(colours(),length(ids))}
plotmapply=function(x,y){if(length(x)>delay){plot(x[length(x),],col=y,add=T,pch=16,cex=1.5)}}
plotmapply1=function(x,y,z,pch){
	if(length(x)>delay){
	lines=Lines(list(Line(list(data.frame(x=coordinates(x)[,1][(length(x)-delay):length(x)],y=coordinates(x)[,2][(length(x)-delay):length(x)])))),ID=1)
	line=SpatialLines(list(lines))
	points=spsample(line,z,type="regular")
	trans=seq(0.1,1,l=z)
	col1=mapply(adjustcolor,col=rep(y,z),alpha.f=trans)
	plot(points,col=col1,add=T,pch=pch,cex=seq(0.1,1.5,l=z))
	}
	if(length(x)==1){plot(x,add=T,pch=pch,col=1)}
}
linemapply=function(x,y){if(length(x)>delay){lines(coordinates(x)[,1][(length(x)-delay):length(x)],coordinates(x)[,2][(length(x)-delay):length(x)],col=y)}}

long=seq(start,length(x),jump)
for(i in 1:length(long)){
plot(bbox[,1],bbox[,2]+c(50,0),type="n",xlab="Longitude",ylab="Latitude",cex.axis=0.8,cex.main=0.8)
rasterImage(n,bbox[1,1],bbox[1,2],bbox[2,1],bbox[2,2])
#Sys.sleep(.3)
rect(bbox[1,1],bbox[2,2]-1000,bbox[2,1]+40,bbox[2,2]+10^6,col="white",border=NA)

temp=x[1:long[i],]
s=split(temp,temp@data[,idn])

#FOR ELIMINATE DIED ANIMALS
if(i>1){
temp0=x[1:long[i-1],]
s0=split(temp0,temp0@data[,idn])
check=unlist(lapply(s0,length))
current=unlist(lapply(s,length))
s=mapply(function(d,ch,cu){if(ch==cu){d[length(d),]}
else{d}},d=s,ch=check,cu=current)
pch=ifelse(unlist(lapply(s,length))>1,16,8)
}

name=unlist(lapply(s,function(x){unique(x@data[,idn])}))
mapply(plotmapply1,s,col[1:length(s)],z=50,pch=pch)
#mapply(linemapply,s,col)
mtext(paste(x@data$datetime[long[i]]),3,-1)
}
}

#################################################
##WHO MIGRATES? NONLINEAR MIXED MODEL APPROACH###
#################################################

#######################
#WHO HAVE ENOUGH DATA##
#######################
eliminate0=function(spdf,id,date="datetime",cut=180,units="days"){
s=split(spdf@data[,date],id)
is=lapply(s,function(x){a=range(x);d=difftime(a[2],a[1],units=units);return(d>cut)})
spdf=split(spdf,id)[unlist(is)]
return(do.call(rbind,spdf))
}

eliminate=function(spdf,id,date="datetime",cut=180,units="days",lag=30,nyearTime="31-01"){
spdf$id=id
spdf$year=strptime(spdf@data[,date],format="%Y-%m-%d %H:%M:%S")$year+1900
spdf$idyear=paste(id,spdf$year)
lid=split(spdf,spdf$id)
#SPLIT YEAR, CONFIRM AN JOIN IF USEFULL
elim=function(x){
lyear=split(x,x$year)
keep=unlist(lapply(lyear,function(x){a=range(x@data[,date]);d=difftime(a[2],a[1],units=units);return(d>cut)}))
res=list()
for(i in 1:length(keep)){
if(i<length(keep)){
if(keep[i]==F){res[[i]]=lyear[[i]]}
if(keep[i]==T & keep[i+1]==T){res[[i]]=lyear[[i]]}
  if(keep[i]==T & keep[i+1]==F & 
     difftime(lyear[[i+1]]@data[1,date],lyear[[i]]@data[length(lyear[[i]]),date],units=units)<lag){
  nyt=as.POSIXct(paste0(nyearTime,"-",lyear[[i+1]]$year[1]),format="%d-%m-%Y")
  res[[i]]=rbind(lyear[[i]],lyear[[i+1]][lyear[[i+1]]@data[,date]<nyt,])
  res[[i]]$idyear=lyear[[i]]$idyear[1]
  }
  else{res[[i]]=lyear[[i]]}
}
else{res[[i]]=lyear[[i]]}
}
res=do.call(rbind,res[keep])
return(res)
}
res=lapply(lid,elim)
res=res[unlist(lapply(res,length))>0]
res=do.call(rbind,res)
return(res)
}

########################################
##PREPARING TIME FOR MIGRATION ANALYSES#
########################################
prep_timeMig=function(spdf,id,join.idyear,date){
y=unlist(lapply(strsplit(join.idyear," "),function(x){x[2]}))
begin=as.POSIXct(paste0(y,"-01-01 00:00:00"),format="%Y-%m-%d %H:%M:%S")
msd_time=difftime(date,begin,units="days")
spdf@data=cbind(spdf@data,msd_time=as.vector(msd_time))
return(spdf)
}


####
#NOMADISM
####
nomadic.model=function(data){
nid=unique(data$idyear)
lnoma=list()
for(i in 1: length(nid)){
z=data[data$idyear==nid[i],]
lnoma[[i]]<- nls(msd ~ 4*D*msd_time,
                    data = z@data,
                    start = c(D = 100))
}
names(lnoma)=nid
nom=data.frame(do.call(rbind,lapply(lnoma,coefficients)))
nom$AIC=unlist(lapply(lnoma,AIC))
return(list(coefs=nom,models=lnoma))
}


####
#SENDENTARY
####
sedentary.model=function(data,xlim=c(0,600)){
nid=unique(data$idyear)
lsede=list()
for(i in 1: length(nid)){
z=data[data$idyear==nid[i],]
lsede[[i]]=try(nls_multstart(msd ~ a+Asym*(1-exp(D*msd_time)),
         data = z@data,
			   iter=500,
			   start_lower=c(a=-10000,Asym=min(z$msd),D=-.1),
         start_upper=c(a=0,Asym=max(z$msd),D=0.1),
			   lower=c(a=-100000,Asym=0,D=-1000),
			   upper=c(a=0,Asym=max(z$msd),D=0),
			   supp_errors="Y"))
print(i)
plot(z@data$msd~z@data$msd_time,cex=.1,xlim=xlim,xlab="Days since first day of the year")
lines(predict(lsede[[i]],data.frame(msd_time=0:365),level=0),col=2)
}
names(lsede)=nid
sed=data.frame(do.call(rbind,lapply(lsede,coefficients)))
sed$AIC=unlist(lapply(lsede,AIC))
return(list(coefs=sed,models=lsede))
}


##############
##DISPERSAL###
##############
dispersal.model=function(data,xlim=c(0,600)){
nid=unique(data$idyear)
ldisp=list()
for(i in 1: length(nid)){
z=data[data$idyear==nid[i],]
ldisp[[i]]=try(nls_multstart(msd ~ Asym/(1+exp((xmid-msd_time)/scal)),
         data = z@data,
			   iter=500,
			   start_lower=c(Asym=mean(z$msd),xmid = 0,scal=0),
         start_upper=c(Asym=max(z$msd),xmid = 365,scal2=100),
			   lower=c(Asym=0,xmid=0,scal=0),
			   upper=c(Asym=max(z$msd),xmid=365,scal=5),
			   supp_errors="Y"))
print(i)
plot(z@data$msd~z@data$msd_time,cex=.1,xlim=xlim,xlab="Days since first day of the year")
lines(predict(ldisp[[i]],data.frame(msd_time=0:365),level=0),col=2)
}
names(ldisp)=nid
disp=data.frame(do.call(rbind,lapply(lapply(ldisp,coefficients),function(x){x[order(names(x))]})))
disp$AIC=unlist(lapply(ldisp,AIC))
return(list(coefs=disp,models=ldisp))
}

###############
##MIGRATION####
###############
migration.model=function(data,xlim=c(0,600)){
nid=unique(data$idyear)
lmig=list()
for(i in 1: length(nid)){
z=data[data$idyear==nid[i],]
lmig[[i]]=try(nls_multstart(msd ~ Asym/(1+exp((xmidA-msd_time)/scal1)) + (-Asym/(1+exp((xmidB-msd_time)/scal2))),
         data = z@data,
			   iter=500,
			   start_lower=c(Asym=mean(z$msd),xmidA = 0,xmidB = 0, scal1 = 0,scal2=0),
                     start_upper=c(Asym=max(z$msd),xmidA = 365, scal1 = 100,xmidB = 365,scal2=100),
			   lower=c(Asym=0,xmidA = 0,xmidB = 0, scal1 = 0,scal2=0),
			   upper=c(Asym=max(z$msd),xmidA = 365,scal1 = 5, xmidB = 365,scal2=5),
			   supp_errors="Y"))
print(i)
plot(z@data$msd~z@data$msd_time,cex=.1,xlim=xlim,xlab="Days since first day of the year")
lines(predict(lmig[[i]],data.frame(msd_time=0:450),level=0),col=2)
}
names(lmig)=nid
mig=data.frame(do.call(rbind,lapply(lapply(lmig,coefficients),function(x){x[order(names(x))]})))
mig$AIC=unlist(lapply(lmig,AIC))
return(list(coefs=mig,models=lmig))
}

#####################
##MIXED-MIGRATION####
#####################
mixed_migration.model=function(data){
nid=unique(data$idyear)
lmig=list()
for(i in 1: length(nid)){
z=data[data$idyear==nid[i],]
lmig[[i]]=try(nls_multstart(msd ~ Asym1/(1+exp((xmidA-msd_time)/scal1)) + (-Asym2/(1+exp((xmidB-msd_time)/scal2))),
         data = z@data,
			   iter=500,
			   start_lower=c(Asym1=mean(z$msd),Asym2=0,xmidA = 0,xmidB = 0, scal1 = 0,scal2=0),
         		   start_upper=c(Asym1=max(z$msd),Asym2=max(z$msd),xmidA = 365, scal1 = 100,xmidB = 365,scal2=100),
			   #lower=c(Asym1=-10^6,Asym2=-10^6,xmidA = 0,xmidB = 0, scal1 = 0,scal2=0),
			   #upper=c(Asym1=0,Asym2=0,xmidA = 365,scal1 = 5, xmidB = 365,scal2=5),
			   supp_errors="Y"))
print(i)
plot(z@data$msd~z@data$msd_time,cex=.1,xlim=c(0,600),xlab="Days since first day of the year")
lines(predict(lmig[[i]],data.frame(msd_time=0:450),level=0),col=2)
}
names(lmig)=nid
mig=data.frame(do.call(rbind,lapply(lapply(lmig,coefficients),function(x){x[order(names(x))]})))
mig$AIC=unlist(lapply(lmig,AIC))
return(list(coefs=mig,models=lmig))
}



###
#PLOT MIGRATION
###
plot.mig=function(data,id,migModel,mode=NULL){
data=data[order(data@data[,id]),]
n=names(migModel)
z=data[data@data[,id] %in% n,]
z=split(z,as.vector(z@data[,id]))
migmodel=migModel[order(names(migModel))]
coef=lapply(migModel,coefficients)
par(mfrow=c(ceiling(sqrt(length(n))),floor(sqrt(length(n)))+1),mar=c(0,0,2,0))
pl=function(x,y,co,mo){
plot(x@data$msd~x@data$msd_time,cex=.1,xaxt="n",yaxt="n",ylab="",xlab="",xlim=c(0,400),ylim=c(0,max(x@data$msd)))
lines(predict(y,data.frame(msd_time=0:600),level=0),col=2)
mtext(co,3,cex=.5)
#mtext(toString(paste(round(co,0),sep="=")),3,cex=.5)#names(co),
mtext(mo,4,-1,cex=.5)
}
mapply(pl,z,migModel,names(migModel),mode)
}


#####################################
#SMOOTHING AND CLEANING TRAJECTORY###
#####################################
clean.smoo=function(data,id="id",date="datetime",quantVel=0.99,windowSize=3){
ids=unique(data@data[,id])
require(zoo)
idsL=list()
for(j in 1:length(ids)){
t=data[data@data[,id]==ids[j],]
clean=as.ltraj(coordinates(t),date=t@data[,date],id=as.vector(t@data[,id]))
clean=clean[[1]]$dist/clean[[1]]$dt<quantile(clean[[1]]$dist/clean[[1]]$dt,quantVel,na.rm=T)
clean[is.na(clean)]=TRUE
cleaned=t[clean==F,]
t1=t[clean,]
xy=data.frame(x=rollapply(coordinates(t1)[,1], windowSize, mean, na.rm = TRUE, partial=T),y=rollapply(coordinates(t1)[,2], windowSize, mean, na.rm = TRUE, partial=T))
idsL[[j]]=SpatialPointsDataFrame(xy,data=data.frame(t1),proj4string=CRS(proj4string(data)))
print(ids[j])
}
return(do.call(rbind,idsL))
}

#####################################
#RESAMPLE TRAJECTORY - TIME##########
#####################################
resample.time=function (trj, stepTime, newFps = NULL, threshold=NULL){
  print(trj$id[1])
require(lubridate)
      if(is.null(threshold)){
	resul=TrajResampleTime(trj,stepTime)
	}
	if(!is.null(threshold) & length(which(as.numeric(diff(trj$time),units="days")>threshold))>0){
	no=which(as.numeric(diff(trj$time),units="days")>threshold)
	times <- seq(from = min(trj$time), to = max(trj$time), by = stepTime)
	#cut=trj$time[c(no,no+1)]
	lcut=lapply(no,function(x){trj$time[c(x,x+1)]})
	lout=lapply(lcut,function(x){times %within% interval(x[1],x[2])})
  keep=do.call(cbind,lout)==F
	kepp=apply(keep,1,sum)==ncol(keep)
	#out=times %within% interval(cut[1],cut[2])
	x <- stats::approx(trj$time, trj$x, times)$y[keep]
	y <- stats::approx(trj$time, trj$y, times)$y[keep]
	resul=TrajFromCoords(data.frame(x, y, times[keep]), timeCol = 3, fps = newFps, 
	spatialUnits = TrajGetUnits(trj), timeUnits = TrajGetTimeUnits(trj))
	}
	else{resul=TrajResampleTime(trj,stepTime)}
  return(resul) 
}

########################
##SELECT MOVEMENT MODE##
########################
#BY AIC
modeAIC=function(resume,delta=NULL){
modes=colnames(resume)
mode_complexity=1:4
names(mode_complexity)=modes
	if(is.null(delta)){mo=apply(resume,1,function(x){modes[which(x==min(x))]})}
	if(is.numeric(delta)){
	mo=apply(resume,1,function(x){
	tie=unlist(x[order(x)])-(unlist(x[order(x)])[1])>delta
	if(all(tie)){choice=modes[order(x)][1]}
	else{choice=names(which(mode_complexity[names(which(tie==F))]==min(mode_complexity[names(which(tie==F))])))}
	return(choice)
	})
	}
return(unlist(mo))
}

#BY CONCORDANCE
concordance=function(data,y,model){
SQres=sum((data@data[,var]-predict(model))^2)
SQtotal=sum((data@data[,var]-mean(data@data[,var]))^2)
SQmodel=sum((predict(model)-mean(predict(model)))^2)
SQreg=sum((data@data[,var]-mean(predict(model)))^2)
n=length(data)
cc=SQres/(SQtotal+SQmodel+n*SQreg)
}

###
#SELECT BY CC
###
#cc_noma=mapply(concordance,split(de2,as.vector(de2$id2)),y="msd",lnoma)
#cc_sede=mapply(concordance,split(de2,as.vector(de2$id2)),y="msd",lsede)
#cc_disp=mapply(concordance,split(de2,as.vector(de2$id2)),y="msd",lmig)
#cc_mig=mapply(concordance,split(de2,as.vector(de2$id2)),y="msd",ldisp)
#cbind(cc_noma,cc_sede,cc_disp,cc_mig)
#modecc=apply(cbind(cc_noma,cc_sede,cc_disp,cc_mig),1,function(x){c("nomadism","sedentary","disperser","migratory")[which(x==min(x))]})

###################################
##EXTRACT DERIVATE METRICS#########
###################################
mig.metrics=function(mig.coef,mode){
mig.coef=mig.coef[mode=="migrater",]
mig.coef$tdep=mig.coef$xmidA-3*mig.coef$scal1
mig.coef$tset=mig.coef$xmidA+3*mig.coef$scal1
mig.coef$tret=mig.coef$xmidB-3*mig.coef$scal2
mig.coef$tarr=mig.coef$xmidB+3*mig.coef$scal2
mig.coef$travelWay=mig.coef$tset-mig.coef$tdep
mig.coef$travelBack=mig.coef$tarr-mig.coef$tret
mig.coef$id=unlist(lapply(strsplit(rownames(mig.coef)," "),function(x){x[1]}))
mig.coef$year=unlist(lapply(strsplit(rownames(mig.coef)," "),function(x){x[2]}))
return(mig.coef)
}

disp.metrics=function(disp.coef,mode){
disp.coef=disp.coef[mode=="disperser",]
disp.coef$tdep=disp.coef$xmid-3*disp.coef$scal
disp.coef$tset=disp.coef$xmid+3*disp.coef$scal
disp.coef$travelWay=disp.coef$tset-disp.coef$tdep
disp.coef$id=unlist(lapply(strsplit(rownames(disp.coef)," "),function(x){x[1]}))
disp.coef$year=unlist(lapply(strsplit(rownames(disp.coef)," "),function(x){x[2]}))
return(disp.coef)
}


###
#FIND JULIAN DATES TO SPLIT TRAJECTORY
###
find.date=function(modeID,type="disperser"){
id=unlist(lapply(strsplit(rownames(modeID)," "),function(x){x[1]}))
year=unlist(lapply(strsplit(rownames(modeID)," "),function(x){x[2]}))
	if(type=="disperser"){
	begin=as.POSIXct(paste0("01-01-",year),format="%d-%m-%Y")
	wend=begin+modeID$tdep*60*60*24
	sbegin=begin+modeID$tset*60*60*24
	resul=data.frame(id,year,wend,sbegin)
	}
	if(type=="migrater"){
	begin=as.POSIXct(paste0("01-01-",year),format="%d-%m-%Y")
	wend=begin+modeID$tdep*60*60*24
	sbegin=begin+modeID$tset*60*60*24
	send=begin+modeID$tret*60*60*24
	nwbegin=begin+modeID$tarr*60*60*24
	resul=data.frame(id,year,wend,sbegin,send,nwbegin)
	}
return(resul)
}


######################
#CLASSIFY TRAJECTORY##
######################
classify.track=function(data,id,date,find.dates,type="disperser"){
data=data[id %in% paste(find.dates$id,find.dates$year),]
id=id[id %in% paste(find.dates$id,find.dates$year)]
ldata=split(data,as.vector(id))
ldate=split(find.dates,paste(find.dates$id,find.dates$year))
if(type=="disperser"){
class=function(x,dates){
x$mode=NA
x$mode[x@data[,date]<dates$wend]="winter"
x$mode[x@data[,date]>=dates$wend & x@data[,date]<=dates$sbegin]="dispersal"
x$mode[x@data[,date]>dates$sbegin]="summer"
return(x)
}
}
if(type=="migrater"){
class=function(x,dates){
x$mode=NA
x$mode[x@data[,date]<dates$wend]="winter1"
x$mode[x@data[,date]>=dates$wend & x@data[,date]<=dates$sbegin]="dispersal1"
x$mode[x@data[,date]>dates$sbegin & x@data[,date]<dates$send]="summer"
x$mode[x@data[,date]>=dates$send & x@data[,date]<=dates$nwbegin]="dispersal2"
x$mode[x@data[,date]>dates$nwbegin]="winter2"
return(x)
}
}
resul=mapply(class,ldata,ldate)
return(do.call(rbind,resul))
}

##########################
#PLOTING MIGRATION DATES##
##########################
kernel.mig=function(mig.metric,disp.metric=NULL,h=10){
require(ks)
#DATES
begin=as.POSIXct("01-01-2020",format="%d-%m-%Y")
end=as.POSIXct("01-01-2021",format="%d-%m-%Y")
departure=begin+mig.metric$xmidA*60*60*24
arrival=begin+mig.metric$xmidB*60*60*42
#ADD "DISPERSAL" DATES
if(!is.null(disp.metric)){departureD=begin+disp.metric$xmid*60*60*24}

#PLOT
departureALL=c(mig.metric$xmidA,disp.metric$xmid)
kdeD=kde(departureALL,h=h,xmin=0,xmax=450)
plot(kdeD,xaxt="n",ylab="Density probability",xlab="Julian day")
points(jitter(rep(0,length(departureALL)),0.02)~departureALL,pch=1,cex=.8)
kdeA=kde(mig.metric$xmidB,h=h,xmin=0,xmax=450)
plot(kdeA,add=T,col=2)
points(jitter(rep(0,length(mig.metric$xmidB)),0.02)~mig.metric$xmidB,pch=1,col=2,cex=.8)
label=toupper(months(seq(begin,end,by="month"),ab=T))
at=difftime(seq(begin,end,by="month"),begin,units="days")
axis(1,label=label,at=at,cex.axis=.8,las=2)

#ADD LINES
pos=seq(max(kdeD$estimate)/2,max(kdeD$estimate),l=nrow(mig.metric))
pmig=function(d,a,p){
lines(c(p,p)~c(d,a),lwd=2,col="gray")
}
mapply(pmig,d=mig.metric$xmidA,a=mig.metric$xmidB,p=pos)
#ADD POINTS
pos=seq(max(kdeD$estimate)/4,(max(kdeD$estimate)/2)-diff(pos)[1],l=nrow(disp.metric))
pmig=function(d,p){
points(p~d,lwd=2,col="gray",pch=16,cex=.8)
}
mapply(pmig,d=disp.metric$xmid,p=pos)

#ISOPLETHS
kiso(kdeD,.95,plot.it=T,col=1)
kiso(kdeD,.50,plot.it=T,cex=2,adj=.5)
kiso(kdeA,.95,plot.it=T,col=2,adj=.3)
kiso(kdeA,.50,plot.it=T,col=2,adj.5)
}

##########################
#NEW PLOTING MIGRATION DATES##
##########################
kernel.mig1=function(mig.metric,disp.metric=NULL,h=10){
require(ks)
#DATES
begin=as.POSIXct("01-01-2020",format="%d-%m-%Y")
end=as.POSIXct("01-01-2021",format="%d-%m-%Y")
departure=begin+mig.metric$tdep*60*60*24
arrival=begin+mig.metric$tret*60*60*42
#ADD "DISPERSAL" DATES
if(!is.null(disp.metric)){departureD=begin+disp.metric$tdep*60*60*24}

#PLOT
departureALL=c(mig.metric$tdep,disp.metric$tdep)
kdeD=kde(departureALL,h=h,xmin=0,xmax=450)
plot(kdeD,xaxt="n",ylab="Density probability",xlab="Julian day")
points(jitter(rep(0,length(departureALL)),0.02)~departureALL,pch=1,cex=.8)
kdeA=kde(mig.metric$tret,h=h,xmin=0,xmax=450)
plot(kdeA,add=T,col=2)
points(jitter(rep(0,length(mig.metric$tret)),0.02)~mig.metric$tret,pch=1,col=2,cex=.8)
label=toupper(months(seq(begin,end,by="month"),ab=T))
at=difftime(seq(begin,end,by="month"),begin,units="days")
axis(1,label=label,at=at,cex.axis=.8,las=2)

#ADD LINES
pos=seq(max(kdeD$estimate)/2,max(kdeD$estimate),l=nrow(mig.metric))
pmig=function(d,a,p){
lines(c(p,p)~c(d,a),lwd=2,col="gray")
}
mapply(pmig,d=mig.metric$tdep,a=mig.metric$tret,p=pos)
#ADD POINTS
pos=seq(max(kdeD$estimate)/4,(max(kdeD$estimate)/2)-diff(pos)[1],l=nrow(disp.metric))
pmig=function(d,p){
points(p~d,lwd=2,col="gray",pch=16,cex=.8)
}
mapply(pmig,d=disp.metric$tdep,p=pos)

#ISOPLETHS
kiso(kdeD,.95,plot.it=T,col=1)
kiso(kdeD,.50,plot.it=T,cex=2,adj=.5)
kiso(kdeA,.95,plot.it=T,col=2,adj=.3)
kiso(kdeA,.50,plot.it=T,col=2,adj.5)

leg=c("Departure winter range","Departure summer range","Dispersal events","Time in summer range")
legend("topright",bty="n",legend=leg,pch=c(1,1,16,NA),lty=c(NA,NA,NA,1),col=c(1,2,8,8),cex=.8)
}

##########################
#NEW PLOTING MIGRATION DATES WITH COLROS##
##########################
#mig.metric=dmig_metrics;disp.metric=ddisp_metrics
kernel.mig.col=function(mig.metric,disp.metric=NULL,h=10,col1="purple4",col2="plum3"){
  require(ks)
  #DATES
  begin=as.POSIXct("01-01-2020",format="%d-%m-%Y")
  end=as.POSIXct("01-01-2021",format="%d-%m-%Y")
  departure=begin+mig.metric$tdep*60*60*24
  arrival=begin+mig.metric$tret*60*60*42
  #ADD "DISPERSAL" DATES
  if(!is.null(disp.metric)){departureD=begin+disp.metric$tdep*60*60*24}
  
  #PLOT
  departureALL=c(mig.metric$tdep,disp.metric$tdep)
  kdeD=kde(departureALL,h=h,xmin=0,xmax=450)
  plot(kdeD,xaxt="n",ylab="Density probability",xlab="Julian day",col=col1)
  points(jitter(rep(0,length(departureALL)),0.02)~departureALL,pch=16,cex=.8,col=col1)
  kdeA=kde(mig.metric$tret,h=h,xmin=0,xmax=450)
  plot(kdeA,add=T,col=col2)
  points(jitter(rep(0,length(mig.metric$tret)),0.02)~mig.metric$tret,pch=16,col=col2,cex=.8)
  label=toupper(months(seq(begin,end,by="month"),ab=T))
  at=difftime(seq(begin,end,by="month"),begin,units="days")
  axis(1,label=label,at=at,cex.axis=.8,las=2)
  
  #ADD LINES
  pos=seq(max(kdeD$estimate)/2,max(kdeD$estimate),l=nrow(mig.metric))
  pmig=function(d,a,p){
    lines(c(p,p)~c(d,a),lwd=2,col=col1,lty=1)
  }
  mapply(pmig,d=mig.metric$tdep,a=mig.metric$tret,p=pos)
  #ADD POINTS
  pos=seq(max(kdeD$estimate)/4,(max(kdeD$estimate)/2)-diff(pos)[1],l=nrow(disp.metric))
  pmig=function(d,p){
    points(p~d,lwd=2,col=col1,pch=1,cex=.8)
  }
  mapply(pmig,d=disp.metric$tdep,p=pos)
  
  #ISOPLETHS
  #kiso(kdeD,.95,plot.it=T,col=1)
  kiso(kdeD,.50,plot.it=T,cex=2,adj=.5,col=col1)
  #kiso(kdeA,.95,plot.it=T,col=2,adj=.3)
  kiso(kdeA,.50,plot.it=T,adj=.5,col=col2)
  
  #leg=c("Departure winter range","Departure summer range","Dispersal events","Time in summer range")
  #legend("topright",bty="n",legend=leg,pch=c(1,1,16,NA),lty=c(NA,NA,NA,1),col=c(1,2,8,8),cex=.8)
}

###############
###CIRCULAR PLOT DATES
######
plot.dates=function(obj,lwd=4,h=20){
layout(matrix(1:5,5,1),heights=c(1,.2,.2,.2,.35))
##DATES PART
par(mar=c(0,4.1,0,2.1))
begin=as.POSIXct("01-01-2020",format="%d-%m-%Y")
end=as.POSIXct("01-01-2021",format="%d-%m-%Y")
o=obj
obj=obj[,c("tdep","tset","tret","tarr")]*60*60*24
plot(seq(begin,end,by="month"),seq(1,nrow(obj),l=13),yaxt="n",
type="n",ylab="",xlab="",xaxt="n")
mtext("Individuals",2,1)
alpha=1
for(i in 1:nrow(obj)){
w1=c(begin+0,begin+obj$tdep[i])
segments(w1[1],i,w1[2],i,lwd=lwd,col=adjustcolor("lightblue",alpha))
t1=c(begin+obj$tdep[i],begin+obj$tset[i])
segments(t1[1],i,t1[2],i,lwd=lwd,col=adjustcolor(gray(.7),alpha))
s=c(begin+obj$tset[i],begin+obj$tret[i])
segments(s[1],i,s[2],i,lwd=lwd,col=adjustcolor("lightgreen",alpha))
t2=c(begin+obj$tret[i],begin+obj$tarr[i])
segments(t2[1],i,t2[2],i,lwd=lwd,col=grey(.2))
w2=c(begin+obj$tarr[i],begin+365*60*60*24)
if(w2[2]>w2[1]){segments(w2[1],i,w2[2],i,lwd=lwd,col=adjustcolor("lightblue",alpha))}
}
##KDE PART
obj=o
kdep=kde(obj$tdep,h=h,xmin=0,xmax=450)
kset=kde(obj$tset,h=h,xmin=0,xmax=450)
kret=kde(obj$tret,h=h,xmin=0,xmax=450)
karr=kde(obj$tarr,h=h,xmin=0,xmax=450)
#cumden=function(vector){
#oe=cumsum(vector[order(vector,decreasing=
#/sum(cumsum(vector[order(vector)]))
#okeep=(1:length(vector))[order(vector,decreasing=T)]
#return(oe[okeep])
#}
par(mar=c(0,4.1,0,2.1))
##DEPARTURE
eval=begin+kdep$eval.points*60*60*24
plot(kdep$estimate~eval,xaxt="n",ylab="",xlab="",type="n",col=gray(.7),yaxt="n",xlim=c(begin,end))
polygon(c(0,eval,0),c(0,kdep$estimate,0),col=gray(.7,.5),border=NA)
points(jitter(rep(0,length(obj$tdep)),0.02)~c(begin+obj$tdep*60*60*24),pch=1,cex=.8,xaxt="n",col=gray(.7))
mtext("Departure winter range",3,-1.5,adj=.03,col=gray(.7))
isod=kiso(kdep,.5,plot.it=T,col=gray(.7))
polygon(begin+isod[[1]]$x*60*60*24,isod[[1]]$y,col=gray(.7,.5),border=NA)
##SETLEMENT
plot(kset$estimate~eval,xaxt="n",ylab="",xlab="",type="n",col="lightgreen",yaxt="n",xlim=c(begin,end))
polygon(c(0,eval,0),c(0,kset$estimate,0),col=adjustcolor("lightgreen",.5),border=NA)
points(jitter(rep(0,length(obj$tset)),0.02)~c(begin+obj$tset*60*60*24),pch=1,cex=.8,xaxt="n",col="lightgreen")
mtext("Settling summer range",3,-1.5,adj=.03,col="lightgreen")
isos=kiso(kset,.5,plot.it=T,col=gray(.7))
polygon(begin+isos[[1]]$x*60*60*24,isos[[1]]$y,col=adjustcolor("lightgreen",.5),border=NA)
mtext("Kernel density probability",2,1)

##RETURN
plot(kret$estimate~eval,xaxt="n",ylab="",xlab="",type="n",col=grey(0.2),yaxt="n",xlim=c(begin,end))
polygon(c(0,eval,0),c(0,kret$estimate,0),col=grey(0.2,0.7),border=NA)
points(jitter(rep(0,length(obj$tret)),0.02)~c(begin+obj$tret*60*60*24),pch=1,cex=.8,xaxt="n",col=grey(0.2))
isor=kiso(kret,.5,plot.it=T,col=gray(.7))
polygon(begin+isor[[1]]$x*60*60*24,isor[[1]]$y,col=grey(.2),border=NA)
mtext("Departure summer range",3,-1.5,adj=.03,col=gray(.2))

##ARRIVAL
par(mar=c(4,4.1,0,2.1))
plot(karr$estimate~eval,xaxt="n",ylab="",xlab="",type="n",col=adjustcolor("lightblue",.5),yaxt="n",xlim=c(begin,end))
polygon(c(0,eval,0),c(0,karr$estimate,0),col=adjustcolor("lightblue",.5),border=NA)
points(jitter(rep(0,length(obj$tar)),0.02)~c(begin+obj$tar*60*60*24),pch=1,cex=1,xaxt="n",col="lightblue")
axis(1,at=seq(begin,end,"month"),label=c(month.abb,month.abb[1]))
isoa=kiso(karr,.5,plot.it=T,col=gray(.7))
polygon(begin+isoa[[1]]$x*60*60*24,isoa[[1]]$y,col=adjustcolor("lightblue",.5),border=NA)
mtext("Julian Dates",1,2.5)
mtext("Resettling winter range",3,-1.5,adj=.03,col="lightblue")
}



####
##CIRCULAR
####
plot.dates.circ=function(obj,inset=c(-.05,-.05),cex=0.8){
require(circlize)
circos.clear()
factors = letters[1]
circos.par(start.degree=90, cell.padding = c(0, 0, 0, 0))
circos.initialize(factors, xlim = c(0, 365))
par(xpd=TRUE,mar=c(0,0,0,0))
begin=as.POSIXct("01-01-2020",format="%d-%m-%Y")
end=as.POSIXct("01-01-2021",format="%d-%m-%Y")
circos.clear()
circos.par("start.degree" = 90, cell.padding = c(0, 0, 0, 0))
circos.initialize(1, xlim = c(0, 365)) # 'a` just means there is one sector
circos.track(ylim = c(1,nrow(obj)), track.height = 0.8, bg.border = NA)
obj=mig_metrics
#LABELS AND SECTORS
at=which(!duplicated(month(seq(begin,end,"day"))))-1
ata=c(at,365)
ata=(360/365)*at
c(ata[4:1],ata[12:3])
ata=for(i in 1:(length(ata))){
draw.sector(c(ata)[i],c(ata,0)[-1][i], clock.wise = F,border=8)
}
circos.axis("top",labels=month.abb,major.at=at,labels.col=grey(.2),major.tick=F)
alpha=1
for(i in 1:nrow(obj)) {
#circos.track(ylim = c(0,10),bg.border = 1)
w1=c(0,obj$tdep[i])
circos.segments(w1[1],i,w1[2],i,lwd=4,col=adjustcolor("lightblue",alpha))
t1=c(obj$tdep[i],obj$tset[i])
circos.segments(t1[1],i,t1[2],i,lwd=4,col=adjustcolor(grey(.7),alpha))
s=c(obj$tset[i],obj$tret[i])
circos.segments(s[1],i,s[2],i,lwd=4,col=adjustcolor("lightgreen",alpha))
t2=c(obj$tret[i],obj$tarr[i])
circos.segments(t2[1],i,t2[2],i,lwd=4,col=adjustcolor(grey(.2)))
w2=c(obj$tarr[i],365)
circos.segments(w2[1],i,w2[2],i,lwd=4,col=adjustcolor("lightblue",alpha))
}
labs=c("Winter range","Spring migration","Summer range","Fall migration")
legend("bottomright",legend=labs,lty=1,lwd=4,col=c("lightblue",grey(.7),"lightgreen",grey(.2)),
cex=cex,bty="n",inset=inset)
title("Migration phenology")
}

####
##ESTIMATE AND PLOT DIRECTIONS
####
plot.dir=function(obj,dist,cex=1){
require(circular)
s=split(obj,obj$idyear)
cw1=do.call(rbind,lapply(s,function(x){colMeans(coordinates(x[x$mode %in% c("winter1"),]),na.rm=T)}))
##FIXING NON CLASSIFIED CW1
w1NA=rownames(cw1[is.na(cw1[,1]),])
swNA=s[w1NA]
cw1NA=do.call(rbind,lapply(swNA,function(x){coordinates(x[x$mode %in% c("dispersal1"),])[1,]}))
cw1[w1NA,]=cw1NA
cs=do.call(rbind,lapply(s,function(x){colMeans(coordinates(x[x$mode %in% c("summer"),]))}))
cw2=do.call(rbind,lapply(s,function(x){colMeans(coordinates(x[x$mode %in% c("winter2"),]))}))
##FIXING NON CLASSIFIED CW2
w2NA=rownames(cw2[is.na(cw2[,1]),])
swNA=s[w2NA]
cw2NA=do.call(rbind,lapply(swNA,function(x){coordinates(x)[nrow(coordinates(x)),]}))
cw2[w2NA,]=cw2NA

#spmig=rbind(cw1,cs)
#spmig=data.frame(spmig,idyear=rownames(spmig))
#rad=unlist(lapply(as.ltraj(spmig[,1:2],id=spmig$idyear,typeII=F),function(x){x$abs.angle[1]}))
#dir=conversion.circular(as.circular(rad,modulo="2pi",type="directions"),units="degree",template="geographics")
rad=atan2((cs-cw1)[,2],(cs-cw1)[,1])
dirs=conversion.circular(as.circular(rad,modulo="2pi",type="directions"),units="degree",template="geographics")
rad=atan2((cw2-cs)[,2],(cw2-cs)[,1])
dirw=conversion.circular(as.circular(rad,modulo="2pi",type="directions"),units="degree",template="geographics")

layout(matrix(c(1,1,2,2,3,4,5,6),4,2))
cws=rbind(cw1,cs)
par(mar=c(4.1,4.1,1,1))
plot(cws,col=rep(c("lightblue","lightgreen"),e=nrow(cws)/2),pch=16,xlab="Easting",ylab="Northing",cex=cex)
for(i in 1:nrow(cw1)){arrows(cw1[i,1],cw1[i,2],cs[i,1],cs[i,2],lty=1,length=.10,col=8)}
#plot(mcp(SpatialPoints(cw1),percent=100),add=T)
#plot(getverticeshr(kernelUD(SpatialPoints(cw1)),95),add=T,border=NA,col=adjustcolor("lightblue",.3))
#plot(getverticeshr(kernelUD(SpatialPoints(cs)),95),add=T,border=NA,col=adjustcolor("lightgreen",.3))
mtext("Spring migration",3,-1.3,adj=0.05,col="lightgreen")

csw=rbind(cs,cw2)
plot(csw,col=rep(c("lightgreen","lightblue"),e=nrow(csw)/2),pch=16,xlab="Easting",ylab="Northing",cex=cex)
for(i in 1:nrow(cw1)){arrows(cs[i,1],cs[i,2],cw2[i,1],cw2[i,2],lty=1,length=.10,col=8)}
mtext("Fall migration",3,-1.3,adj=0.05,col="lightblue")

par(mar=c(4.1,4.1,1,1))
hist(c(mig_metrics$Asym,disp_metrics$Asym)/1000,xlab="Migration distance (km)",
     main="",col=adjustcolor("lightgreen",.5),freq=F,border="lightgreen")
par(mar=c(0,0,0,0))
rose.diag(dirs, bins = 18,col=adjustcolor("lightgreen",.5),border="lightgreen",shrink=1)
points(dirs,col="lightgreen",cex=cex)

par(mar=c(4.1,4.1,1,1))
hist(dist/1000,xlab="Migration distance (km)",
     main="",col=adjustcolor("lightblue",.5),freq=F,border="lightblue")
par(mar=c(0,0,0,0))
rose.diag(dirw, bins = 18,col=adjustcolor("lightblue",.5),border="lightblue",tol=0)
points(dirw,col="lightblue",cex=cex)
}

####
##ESTIMATE AND PLOT DIRECTIONS TWO SPECIES
####
#col1=c("mediumpurple4","plum4")
#col2=c("steelblue4","steelblue3")
plot.dir2=function(obj,cex=1,col1=c(1,gray(.5)),col2=c(1,gray(.5))){
  require(circular)
  s=split(obj,obj$idyear)
  cw1=do.call(rbind,lapply(s,function(x){colMeans(coordinates(x[x$mode %in% c("winter1"),]),na.rm=T)}))
  ##FIXING NON CLASSIFIED CW1
  w1NA=rownames(cw1[is.na(cw1[,1]),])
  swNA=s[w1NA]
  cw1NA=do.call(rbind,lapply(swNA,function(x){coordinates(x[x$mode %in% c("dispersal1"),])[1,]}))
  cw1[w1NA,]=cw1NA
  cs=do.call(rbind,lapply(s,function(x){colMeans(coordinates(x[x$mode %in% c("summer"),]))}))
  cw2=do.call(rbind,lapply(s,function(x){colMeans(coordinates(x[x$mode %in% c("winter2"),]))}))
  ##FIXING NON CLASSIFIED CW2
  w2NA=rownames(cw2[is.na(cw2[,1]),])
  swNA=s[w2NA]
  cw2NA=do.call(rbind,lapply(swNA,function(x){coordinates(x)[nrow(coordinates(x)),]}))
  cw2[w2NA,]=cw2NA
  
  #spmig=rbind(cw1,cs)
  #spmig=data.frame(spmig,idyear=rownames(spmig))
  #rad=unlist(lapply(as.ltraj(spmig[,1:2],id=spmig$idyear,typeII=F),function(x){x$abs.angle[1]}))
  #dir=conversion.circular(as.circular(rad,modulo="2pi",type="directions"),units="degree",template="geographics")
  rad=atan2((cs-cw1)[,2],(cs-cw1)[,1])
  dirs=conversion.circular(as.circular(rad,modulo="2pi",type="directions"),units="degree",template="geographics")
  rad=atan2((cw2-cs)[,2],(cw2-cs)[,1])
  dirw=conversion.circular(as.circular(rad,modulo="2pi",type="directions"),units="degree",template="geographics")
  
  #layout(matrix(c(1,1,2,2,3,4,5,6),4,2))
  cws=rbind(cw1,cs)
  par(mar=c(4.1,4.1,1,1))
  idsp=unique(obj@data[,c("idyear","sp")])
  sp=idsp$sp[match(rownames(cws),idsp$idyear)]
  #data(wrld_simpl)
  #wo=wrld_simpl
  #wo=spTransform(wo,CRS=CRS("+proj=longlat +ellps=WGS84"))
  #usca=wo[wo$NAME %in% c("United States","Canada"),]
 #plot(usca,border="gray",add=T) 
 #plot_usmap(include = c("MN", "WI"))
 #plot(z)
  #  cwsl=spTransform(SpatialPoints(na.omit(data.frame(cws[,1],cws[,2])),proj4string=CRS("+proj=utm +zone=16 +north")),CRS=CRS(proj4string(usca)))
  layout(matrix(1:4,2,2,byrow=T))
  col=rep(col1,e=nrow(cws)/2)
  par(mar=c(4,4,1,0),xpd=T)
  plot(cws,col=col,pch=ifelse(sp=="m",NA,16),
       xlab="Easting",ylab="Northing",cex=cex,axes=T)
  for(i in 1:nrow(cw1)){arrows(cw1[i,1],cw1[i,2],cs[i,1],cs[i,2],lty=ifelse(sp[i]=="m",0,1),length=.10,col=adjustcolor(col1[2],.7))}
  text(279000,5317500,"A")
  text(295000,5315000,"B")
  #segments(270000,5312500,300000,5317500)
  par(mar=c(2,2,2,2),xpd=T)
  legend("topleft",legend=c("Winter range","Summer range"),pch=16,col=col1,
         bty="n",cex=.6,pt.cex=1.5)
  rose.diag(dirs[sp=="d"], bins = 18,col=adjustcolor(col1[1],.8),border=col1[1],shrink=1)
  points(dirs[sp=="d"],col=col1[1],cex=cex)
  imageOnPlot(x=-1,y=1.2,sp="deer",scale=c(.3,.3),col=col1[1])
  
  col=rep(col2,e=nrow(cws)/2)
  par(mar=c(4,4,1,0),xpd=T)
  plot(cws[sp=="m",],col=col[sp=="m"],xlab="Easting",ylab="Northing",cex=cex,axes=T,pch=16)
  for(i in 1:nrow(cw1)){arrows(cw1[i,1],cw1[i,2],cs[i,1],cs[i,2],lty=ifelse(sp[i]=="m",1,0),length=.10,col=adjustcolor(col2[2],.7))}
  #segments(270000,5312500,300000,5317500)
  text(279000,5317500,"A")
  text(295000,5315000,"B")
  par(mar=c(2,2,2,2),xpd=T)
  legend("topleft",legend=c("Winter range","Summer range"),pch=16,col=col2,
         bty="n",cex=.6,pt.cex=1.5)
  rose.diag(dirs[sp=="m"], bins = 18,col=adjustcolor(col2[1],.8),border=col2[1],shrink=1)
  points(dirs[sp=="m"],col=col2[1],cex=cex)
  imageOnPlot(x=-1,y=1.2,sp="moose",scale=c(.3,.3),col=col2[1])
}

####
##ESTIMATING ISOPLETHS
####
kiso=function(kde,iso=.95,plot.it=F,col=1,cex=1,pch=16,disj=1,adj=.3){
volume=kde$estimate/sum(kde$estimate)
o=order(volume,decreasing=T)[cumsum(volume[order(volume,decreasing=T)])<iso]
days=sort(kde$eval.points[o])
#POLYGON
djoin=mean(diff(kde$eval.points))
cores=split(days, cumsum(c(TRUE, diff(days)!=djoin*disj)))
corepol=lapply(cores,function(x){
x=c(x[1],x,x[length(x)])
y=c(0,kde$estimate[kde$eval.points %in% x],0)
return(data.frame(x,y))
})
#PLOT
if(plot.it==T){lapply(corepol,function(x){polygon(x$x,x$y,col=adjustcolor(col,adj),border=NA)})}
##points(rep(0,length(days))~days,col=col,cex=cex,pch=pch)
return(corepol)
}


########
#BURN OUT FIRST LOCATIONS POST  CAPTURES
########
burn.out=function(data,id="id",date="datetime",cut=7){
  s=split(data,data@data[,id])
  out=lapply(s,function(x){x[x@data[,date]>(x@data[1,date]+cut*60*60*24),]}) 
  return(do.call(rbind,out))
}


####
#GETVOLUME TO BUILD ISOPLETHS
####
rastToIso=function(s){
ss <-  as.matrix(s/sum(as.matrix(s),na.rm=T))
index <- 1:length(ss)
vord <- ss[order(ss,decreasing=T)]
idord <- index[order(ss,decreasing=T)]
vsu <- cumsum(vord)
vreord <- vsu[order(idord)]
reg <- raster(matrix(vreord,dim(ss)[1],dim(ss)[2]))
extent(reg)=extent(s)
return(reg)
}


####
#DEALING WITH RANDOM COEFS  AN INTERACTIONS- GLMMTMB
####
#habs=unique(ssf1$hab)
#habs=habs[-which(habs=="mf")]

#coef=coefficients(full.fit)
#cpop=coef$cond$strata
#crand=coef$cond$id[,-1]
#crand=unlist(coef$cond$strata[1,-1])

#crand=reshape(crand, idvar = "id", ids = row.names(crand),
#              times = names(crand), timevar = "par",
#              varying = list(names(crand)), direction = "long")

#for(i in 1:length(habs))
#h=habs[i]
#lab=paste(h,c("win","dis","sum","wolf"),sep=" ")
#col=crand[,grep(h,colnames(crand))]
#col=cbind(col[,1],col[,1]+col[,2:ncol(col)])
#colnames(col)=lab
#col=reshape(col, idvar = "id", ids = row.names(col),
#              times = names(col), timevar = "par",
#              varying = list(names(col)), direction = "long")
#vioplot(col[,2]~col$par,las=2,cex.axis=.6,xlab="",col="white",ylab="Selection strength")
#abline(h=0,lty=3)
#pars=levels(factor(col$par))
#for(i in 1:length(pars)){
#  points(jitter(rep(i,length(col[col$par==pars[i],2])),.4),
#         col[col$par==pars[i],2],pch=21,col="grey",cex=.4,bg="grey")
#}


####
##SIMPLE RANDOM EFFECT
####
#habs=unique(ssf1$hab)
#habs=habs[-which(habs=="mf")]

#coef=coefficients(full.fit)
#cpop=coef$cond$strata
#crand=unlist(coef$cond$strata[1,-1])

#crand=reshape(crand, idvar = "id", ids = row.names(crand),
#              times = names(crand), timevar = "par",
#              varying = list(names(crand)), direction = "long")

#for(i in 1:length(habs))
#  h=habs[i]
#lab=paste(h,c("win","dis","sum","wolf"),sep=" ")
#col=crand[,grep(h,colnames(crand))]
##col=cbind(col[,1],col[,1]+col[,2:ncol(col)])
#colnames(col)=lab
#col=reshape(col, idvar = "id", ids = row.names(col),
#            times = names(col), timevar = "par",
#            varying = list(names(col)), direction = "long")
#vioplot(col[,2]~col$par,las=2,cex.axis=.6,xlab="",col="white",ylab="Selection strength")
#abline(h=0,lty=3)
#pars=levels(factor(col$par))
#for(i in 1:length(pars)){
#  points(jitter(rep(i,length(col[col$par==pars[i],2])),.4),
#         col[col$par==pars[i],2],pch=21,col="grey",cex=.4,bg="grey")
#}


###
#LINE TO RASTER HOLDING DISTANCES
##
require(rgeos)
#r=raster(extent(map),res=res(map))
#dd = gDistance(Sl, as(r,"SpatialPoints"), byid=TRUE)
#r[] = apply(dd,1,min)
#plot(r)

###
##GRAPHS MODES
###
plot.ssfMode=function(x,y,wp=c("actual","min")){
hab=c("mf","df","open","urb","wet","wat")
mode=c("w","d","s")
sl=c(x$frame$slayer,y$frame$slayer)
if(wp=="actual"){wol=mean(c(x$frame$wolf_pressure,y$frame$wolf_pressure),na.rm=T)}
if(wp=="min"){wol=min(c(x$frame$wolf_pressure,y$frame$wolf_pressure),na.rm=T)}
h=expand.grid(hab=hab,mode=mode,slayer=mean(sl,na.rm=T),id=NA,strata=NA,
  wolf_pressure=wol)
nd=dummy_cols(h,select_columns="hab",ignore_na=T)
nd=dummy_cols(nd,select_columns="mode",ignore_na=T)
pred=predict(x,nd,se.fit=T,type="response")
nd$p=pred$fit
nd$error=pred$se.fit
deer=nd  
pred=predict(y,nd,se.fit=T,type="response")
nd$p=pred$fit
nd$error=pred$se.fit
moose=nd
#ELEVATION FOR DEER
if(wp=="actual"){
cod=confint(x)
cod=t(cod[grep("slayer",rownames(cod)),])
ed=cbind(w=cod[,1],cod[,1]+cod[,-1])[,-4][,c(1,3,2)]
ed=exp(ed)
}
if(wp=="min"){
  cod=confint(x)
  wod=t(cod[grep("slayer:wolf",rownames(cod)),])*wol
  cod=t(cod[grep("slayer",rownames(cod)),])
  ed=cbind(w=cod[,1],cod[,1]+cod[,-1])[,-4][,c(1,3,2)]
  ed=apply(ed,2,function(x){x+wod})
  ed=exp(ed)
}

#ELEVATION FOR MOOSE
if(wp=="actual"){
cod=confint(y)
cod=t(cod[grep("slayer",rownames(cod)),])
em=cbind(w=cod[,1],cod[,1]+cod[,-1])[,-4][,c(1,3,2)]
em=exp(em)
}

if(wp=="min"){
  cod=confint(y)
  wod=t(cod[grep("slayer:wolf",rownames(cod)),])*wol
  cod=t(cod[grep("slayer",rownames(cod)),])
  em=cbind(w=cod[,1],cod[,1]+cod[,-1])[,-4][,c(1,3,2)]
  em=apply(em,2,function(x){x+wod})
  em=exp(em)
}

par(mfrow=c(3,1),mar=c(1,4,0,0))
ylim=range(c(deer$p+1.96*deer$error,deer$p-1.96*deer$error,moose$p+1.96*moose$error,moose$p-1.96*moose$error,
             range(rbind(ed,em))))
for(i in 1:length(mode)){
  plot(c(0.7,7.4),ylim,ylab="Selection strength",xlab="Habitat",xaxt="n",type="n")
  axis(1,1:7,c(hab,"Elevation"))
  abline(h=1,lty=3)
dt=deer[deer$mode==mode[i],]  
points(1:length(hab),dt$p)  
segments(1:length(hab),dt$p-1.96*dt$error,1:length(hab),dt$p+1.96*dt$error) 
points(7,ed[3,i])
segments(7,ed[1,i],7,ed[2,i])
mt=moose[moose$mode==mode[i],]  
adj=.2
points(1:length(hab)+adj,mt$p,pch=16)  
segments(1:length(hab)+adj,mt$p-1.96*mt$error,1:length(hab)+adj,mt$p+1.96*mt$error) 
points(7.2,em[3,i],pch=16)
segments(7.2,em[1,i],7.2,em[2,i])
}
return(list(deer,moose,ed,em))
}

######
##SHORT CUT SSF PLOT
######
plot.ssfMode2=function(list,par=T,col=rep(1,3),col.wolf=1,yrange=NULL,pchsp=c(1,22,2,17,16,15),cex=1.3,legend=F,wolf=T,letters=LETTERS[1:3]){
  deer=list[[1]]
  moose=list[[2]]
  ed=list[[3]][,3:1]
  em=list[[4]][,3:1]
  mode=rev(unique(deer$mode))
  hab=colnames(deer)[grep("hab_",colnames(deer))]
  hab=unlist(lapply(strsplit(hab,"_"),function(x){x[[2]]}))
  if(par){layout(matrix(c(1,2,3),3,1),heights = c(1,1,1.4))}
  
  ylim=range(c(deer$p+1.96*deer$error,deer$p-1.96*deer$error,moose$p+1.96*moose$error,moose$p-1.96*moose$error,
               range(rbind(ed,em))))
  if(is.null(yrange)){yrange=range(ylim)}
  
  col=col
  for(i in 1:length(mode)){
    if(i==1){
      par(mar=c(0,4,3,0),xpd=T)
      plot(c(0.7,7.4),ylim,ylab="",xlab="Habitat",xaxt="n",type="n",ylim=yrange)
      pchsp2=pchsp[5:6]
      #if(legend){legend("topright",legend=c("Deer","Moose"),pch=c(16,15),bty="n",pt.cex=1.3,col=col)}
      if(legend){imageOnPlot(7.2,2.1,sp="deer",scale=c(.3,.15),col="mediumpurple4")}
      if(legend){imageOnPlot(6.5,2.1,sp="moose",scale=c(.3,.2),col="steelblue4")}
      if(wolf){imageOnPlot(x=1,y=2.7,sp="wolf",scale=c(.5,.03),col=col.wolf)}
      else(imageOnPlot(x=1,y=2.7,sp="nowolf",scale=c(.5,.03),col=col.wolf))
      mtext(letters[1],3,-1.3,adj=.02)
    }
    if(i==2){
      par(mar=c(0,4,1,0))
      plot(c(0.7,7.4),ylim,ylab="Selection strength",xlab="Habitat",xaxt="n",type="n",ylim=yrange)
      pchsp2=pchsp[3:4]
      mtext(letters[2],3,-1.3,adj=.02)
    }
    if(i==3){
      par(mar=c(4,4,1,0))
      plot(c(0.7,7.4),ylim,ylab="",xlab="Habitat",xaxt="n",type="n",ylim=yrange)
      pchsp2=pchsp[1:2]
      axis(1,1:7,c(hab,"Elevation"),las=1,cex.axis=.9)
      mtext(letters[3],3,-1.3,adj=.02)
    }
    par(xpd=F)
    abline(h=1,lty=3)
    dt=deer[deer$mode==mode[i],]
    #imageOnPlot(1:length(hab),dt$p,"deer",scale=c(0.1,0.2))
    points(1:length(hab),dt$p,col=col[1],pch=pchsp2[1],cex=cex)  
    segments(1:length(hab),dt$p-1.96*dt$error,1:length(hab),dt$p+1.96*dt$error,col=col[1]) 
    points(7,ed[3,i],col=col[1],pch=pchsp2[1],cex=cex)
    segments(7,ed[1,i],7,ed[2,i],col=col[1])
    
    mt=moose[moose$mode==mode[i],]  
    #imageOnPlot(1:length(hab)+adj,mt$p,"moose",scale=c(0.1,0.2))
    adj=.2
    points(1:length(hab)+adj,mt$p,col=col[2],pch=pchsp2[2],cex=cex) 
    segments(1:length(hab)+adj,mt$p-1.96*mt$error,1:length(hab)+adj,mt$p+1.96*mt$error,col=col[2]) 
    points(7.2,em[3,i],col=col[2],pch=pchsp2[2],cex=cex)
    segments(7.2,em[1,i],7.2,em[2,i],col=col[2])
  }
  
}
#####
##PLOT FOR PLOT IMAGES ON PLOT
####
imageOnPlot=function(x=0,y=0,sp="wolf",scale=c(1,1),col=2){
  require(png)
  if(sp=="wolf"){img=readPNG("wolf.png")}
  if(sp=="deer"){img=readPNG("deer.png")}
  if(sp=="moose"){img=readPNG("moose.png")}
  if(sp=="nowolf"){img=readPNG("nowolf.png")}
  if(sp=="graydeer"){img=readPNG("graydeer.png")}
  if(sp=="graymoose"){img=readPNG("graymoose.png")}
  img[img!=1]=0
  img=as.raster(img,interpolate=T)
  if(sp=="deer"){img[img==names(sort(table(img),dec=T))[3]]=col}
  if(sp!="deer"){img[img==names(sort(table(img),dec=T))[2]]=col}
  #plot(img)
  for(i in 1:length(x)){
rasterImage(img,x[i]-scale[1],y[i]-scale[1],x[i]+scale[1],y[i]+scale[2],col=col)
  }
}


#####
#DUMMY MAPS
#####
dummy.raster=function(map,table,season="w",scale=T,predator=c("min","max","mean","actual")){
raster=map[[1]]
raster[raster==table$index[table$hab=="crop"]]=NA
cont=map[[-1]]
if(scale){cont=scale(cont)}
layer=cont[[1]]
pred=cont[["wolf_pressure"]]
if(predator=="actual"){pred=pred}
if(predator=="min"){pred[]=min(matrix(pred),na.rm=T)}
if(predator=="max"){pred[]=max(matrix(pred),na.rm=T)}
if(predator=="mean"){pred[]=mean(matrix(pred),na.rm=T)}
cont=stack(layer,pred)
osea=c("w","d","s")
osea=osea[osea!=season]
rlist=list()
lab=paste0("hab_",table$hab)  
for(i in 1:length(table$hab)){
rlist[[i]]=raster==(table$index[table$hab==table$hab[i]])
}
names(rlist)=lab
res=stack(rlist)
raster[]=1
assign(paste0("mode_",season),raster)
raster[]=0
assign(paste0("mode_",osea[1]),raster)
assign(paste0("mode_",osea[2]),raster)
lsea=list(mode_w=mode_w,mode_d=mode_d,mode_s=mode_s)
raster[]=NA
lrand=list(strata=raster,id=raster)
res=stack(cont,res,stack(lsea),stack(lrand))
return(res)
}


######################################
#PREPARE DATA FOR SPATIAL PREDICTION##
######################################
spat.pred=function(maps,season="w",predator=c("min","max","mean","actual")){
require(fastDummies)
#RASTER STACK
nd=values(maps)
nd=data.frame(nd,coordinates(maps))
nd=na.omit(nd)
hab=data.frame(index=1:8,hab=c("wat","urb","open","df","cf","mf","wet","crop"))
dumhab=dummy.raster(maps,tab=hab,season=season)
varrast=stack(maps,dumhab)
#px=predict(varrast,x,ra.rm=F)
nd$slayer=scale(nd$slayer)
if(predator=="actual"){nd$wolf_pressure=scale(nd$wolf_pressure)}
if(predator=="min"){nd$wolf_pressure=min(scale(nd$wolf_pressure))}
if(predator=="max"){nd$wolf_pressure=max(scale(nd$wolf_pressure))}
if(predator=="mean"){nd$wolf_pressure=mean(scale(nd$wolf_pressure))}
nd$hab=hab$hab[match(nd$hab,hab$index)]
nd=nd[nd$hab!="crop",]
nd$hab=as.vector(nd$hab)
nd$mode=factor(rep(season,nrow(nd)),levels=c("w","s","d"))
nd=dummy_cols(nd,select_columns="hab",ignore_na=T)
nd=dummy_cols(nd,select_columns="mode",ignore_na=T)
nd$strata=NA
nd$id=NA
return(nd)
}

##########################
##GLMMTMB IN PARALLEL#####
##########################
predict.parallel=function(nd,model,clu=100,core=5){
require(parallel)
clu=100
par=rep(1:clu,e=round(nrow(nd)/clu),clu)
nd$list=par[1:nrow(nd)]
cl <- makeCluster(core,type="FORK")
clusterEvalQ(cl, {library(glmmTMB)})
clusterExport(cl, deparse(quote(model)))
ndl=split(nd,nd$list)
t=Sys.time()
res=parLapply(cl,ndl,function(y){predict(model,y)})
Sys.time()-t
stopCluster(cl)
nd$p=do.call(c,res)
map=rasterFromXYZ(nd[,c("x","y","p")])
map=rastToIso(exp(map))
return(list(nd,map))
}

########################
##GLMMTMB IN LAPPLY#####
########################
predict.lapply=function(nd,model,clu=100){
  t=Sys.time()
  clu=100
  par=rep(1:clu,e=round(nrow(nd)/clu),clu)
  nd$list=par[1:nrow(nd)]
  ndl=split(nd,nd$list)
  res=lapply(ndl,function(x){print(x$list[1]);return(predict(model,x))})
  nd$p=do.call(c,res)
  map=rasterFromXYZ(nd[,c("x","y","p")])
  map=rastToIso(exp(map))
  print(Sys.time()-t)
  return(list(nd,map))
}

#####################
##PLOT SPATIAL MODEL#
#####################
plot.space=function(nr,level=.5,col.wolf=1,col=c(NA,"lightblue","lightgreen",2),par=T,axis=T,wolf=T,sea=T){
  if(par==T){par(mfrow=c(dim(nr)[3],1))}  
  for(i in 1:dim(nr)[3]){
    #if(colour){image(nr[[i]],col=col,ylab="Northing",xlab="Easting")}  
    #else{image(nr[[i]],col=adjustcolor(c(NA,"gray",1,2),.9),ylab="Northing",xlab="Easting")
    if(i<dim(nr)[3] & axis==T){
      xaxt="n";yaxt="n";ylab="";xlab=""
      if(i==1){par(mar=c(0,4,2,0),xpd=T)}
      if(i>1){par(mar=c(0,4,0,0),xpd=T)}
      image(nr[[i]],col=adjustcolor(col,.9),ylab=xlab,xlab=xlab,xaxt=xaxt,yaxt=yaxt)
      #if(i==1){legend(x=239000,y=5410000,legend="Core overlap",pch=15,bty="n")}
      if(i==1){legend(x=239000,y=5410000,legend=c("Moose range","Deer range","Core overlap"),
                      pch=15,bty="n",col=c(col[2],col[3],col[4]),pt.cex=1.2)}
      if(i==1 & wolf==T){imageOnPlot(x=247000,y=5430000,sp="wolf",scale=c(8000,10000),col=col.wolf)}
      if(i==1 & wolf==F){imageOnPlot(x=247000,y=5430000,sp="nowolf",scale=c(8000,10000),col=col.wolf)}
    }
    if(i<dim(nr)[3] & axis==F){
      xaxt="n";yaxt="n";ylab="";xlab=""
      if(i==1){par(mar=c(0,0,2,2),xpd=T)}
      if(i>1){par(mar=c(0,0,0,2),xpd=T)}
      image(nr[[i]],col=adjustcolor(col,.9),ylab=xlab,xlab=xlab,xaxt=xaxt,yaxt=yaxt)
      if(i==1 & wolf==T){imageOnPlot(x=247000,y=5430000,sp="wolf",scale=c(8000,10000),col=col.wolf)}
      if(i==1 & wolf==F){imageOnPlot(x=247000,y=5430000,sp="nowolf",scale=c(8000,10000),col=col.wolf)}
    }
    if(i==dim(nr)[3] & axis==T){
      par(mar=c(4,4,0,0),xpd=T)
      image(nr[[i]],col=adjustcolor(col,.9),ylab="Northing",xlab="Easting")
    }
    if(i==dim(nr)[3] & axis==F){
      par(mar=c(4,0,0,2),xpd=T)
      xaxt="n";yaxt="n";ylab="";xlab=""
      image(nr[[i]],col=adjustcolor(col,.9),ylab=ylab,xlab=xlab,xaxt=xaxt,yaxt=yaxt)
    }
    
    ##
    #PLOT MOOSE AND DEER
    ##
    if(i==1){
      imageOnPlot(x=250000,y=5340000,sp="graymoose",scale=c(9000,10000),col="steelblue4")
      imageOnPlot(x=320000,y=5320000,sp="graydeer",scale=c(10000,10000),col="mediumpurple4")
    }
    if(i!=1){
      imageOnPlot(x=250000,y=5340000,sp="graydeer",scale=c(9000,10000),col="mediumpurple4")
      imageOnPlot(x=320000,y=5320000,sp="graymoose",scale=c(10000,10000),col="steelblue4")
    }
    
    vals <- getValues(nr[[i]])
    areaOver=prod(length(subset(vals, vals == 3)),res(nr[[i]]))/prod(length(vals),res(nr[[i]]))
    if(sea==F){mtext(paste0(round(areaOver*100,1),"%"),3,-1.5,adj=0,cex=.7)}
    if(sea==T){
      s=c("Winter range","Transience time","Summer range")[i]
      #mtext(paste(s,paste0("(",round(areaOver*100,1),"%",")")),3,-1.5,adj=0,cex=.7)
    }
  }
  return(nr)
}


####
##CUT TRAJECTORY BEFORE DEATH
####
cut.death=function(spdf,got=60,time=60){
got=got*60*60*24
time=got+(time*60*60*24)
return(spdf[spdf$datetime>=spdf$death-time & spdf$datetime<spdf$death-got,]) 
}