#####
#HOW TO ESTIMATE THE SURFACE OF PREDATION RISK BASED ON WOLF TRAJECTORIES
####


###
#LIBRARIES
###
library(rgdal)
library(rgeos)
library(ggmap)
library(png)
library(adehabitatLT)
library(adehabitatHR)
library(animation)
library(lattice)
library(nls.multstart)
library(raster)
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=strptime(mo$datetime,format="%Y-%m-%d %H:%M:%S")$year+1900
mo$idyear=paste(mo$id,mo$year)
mo=SpatialPointsDataFrame(mo[,c("Easting","Northing")],data=mo,proj4string=CRS(proj))

#LOAD SIMULATED DEER DATA
de=read.csv("deer_track.csv")
de$id=de$GPDeer
de$year=strptime(de$datetime,format="%Y-%m-%d %H:%M:%S")$year+1900
de$idyear=paste(de$id,de$year)
de=SpatialPointsDataFrame(de[,c("Easting","Northing")],data=de,proj4string=CRS(proj))


#LOAD SIMULATED WOLF DATA
wo=read.csv("wolf_track.csv")
wo$id=wo$GPwolf
wo$year=strptime(wo$datetime,format="%Y-%m-%d %H:%M:%S")$year+1900
wo$idyear=paste(wo$id,wo$year)
wo=SpatialPointsDataFrame(wo[,c("Easting","Northing")],data=wo,proj4string=CRS(proj))


###
#DEFINE STUDY AREA
###
join=rbind(de[,"id"],mo[,"id"])
mcp=mcp(join[,-1],percent=100)
mcpB=gBuffer(mcp,width=0)
plot(wo,axes=T,cex=.1)
plot(mcpB,axes=T,add=T,border=2)

###
#HOME RANGE AND OVERLAP AND PACK CLUSTER NETWORKING
###
k=kernelUD(wo[,"GPwolf"],same4all=T,grid=1000)
pol=getverticeshr(k,95)
plot(pol,axes=T)

over=kerneloverlaphr(k,method="BA")
over

cov=unique(wo[,c("GPwolf","Pack")]@data)
cov$GPwolf
cov$strata=paste(cov$Pack,1:length(cov$Pack))
wo$strata=cov$strata[match(paste(wo$GPwolf,wo$Pack),paste(cov$GPwolf,cov$Pack))] 
rownames(over)=colnames(over)=cov$strata
pol@data$idd=as.vector(unique(wo$strata))


#NETWORK
library(igraph)
neto=graph_from_adjacency_matrix(over,mode="undirected",diag=F,weighted=T)
#V(neto)$size=pol$wei/10
E(neto)$width=sqrt((E(neto)$weight*30))
plot(neto)

par(mar=rep(0,4))
plot(neto,vertex.color="gray",edge.curved=0.1,
     edge.color=1, vertex.label.cex=.6)

#CUTTING WEAK CONNECTIONS
neto2=delete_edges(neto,E(neto)[weight<quantile(E(neto)$weight,.50)])
e=ends(neto2, es=E(neto2), names=T)
#org=data.frame(vhf=V(neto2)$name,sex=pol$sex,size=pol$wei)
#a=org[match(e[,1],org$vhf),"sex"]
#b=org[match(e[,2],org$vhf),"sex"]
#pair=paste0(a,b)
#col.edge=ifelse(pair=="FF","red",ifelse(pair=="MM","black",ifelse(pair=="MF"|pair=="FM","yellow","gray")))
par(mar=c(0,0,0,0))
#plot(neto2,vertex.color=ifelse(pol$sex=="M","gray50","tomato"),edge.curved=.1,
#     edge.color=col.edge, vertex.label.cex=.6)

plot(neto2,vertex.color="gray50",edge.curved=.1,
        edge.color=1, vertex.label.cex=.6)
     

#PREDICTING GROUPS
neto2=delete_edges(neto,E(neto)[weight<quantile(E(neto)$weight,.80)])
cebo=cluster_label_prop(neto2)
cebo=cluster_fast_greedy(neto2)
#jpeg("sociobiologyHR.jpeg",quality=100,units="cm",h=20,w=30,res=500)
#layout(mat=matrix(c(1,1,2,3,4,5,6,7,8), ncol=3),widths=c(1,.7))
par(mar=c(0,0,0,0),mfrow=c(1,2))
cols=adjustcolor(c("orange","lightblue","lightgreen","purple","tomato","grey"),1)
plot(cebo,neto2,col="gray50",mark.col=cols
     ,mark.border=NA,vertex.label.cex=.6,edge.curved=0)
ids=unique(wo$GPwolf)
#plot(mcp(wo,50),border=F,axes=T,xaxt="n",yaxt="n")
wcenter=spTransform(gCentroid(wo),CRS=CRS("+proj=longlat"))
goomap(wcenter,zoom=10,crs=CRS(proj4string(wo)),adjlat=0.1,adjlon=0.1,nosignal=2*10^3)
for(i in 1:length(cebo)){
  par(xpd=F)
  #plot(bbox[,1],bbox[,2],type="n",cex.axis=0.8,cex.main=0.8,xaxt="n",yaxt="n",xlab="",ylab="")
  #rasterImage(n,bbox[1,1],bbox[1,2],bbox[2,1],bbox[2,2])
  plot(wo[wo@data$strata %in% cebo[[i]],],col=cols[i],axes=T,cex=.2,add=T,pch=16,xpd=F)
  #plot(pol[pol$idd%in% cebo[[i]],],border=cols[i],add=T,lwd=3)
}

#####
#TERRITORY UTILIZATION DENSITY
####
k=kernelUD(wo[,"strata"],same4all=T,grid=2000,ex=.1)
lapply(k,function(x){res(raster(x))})
sk=lapply(k,function(x){raster(x)/sum(as.matrix(x))})
r=stack(sk)
lr=list()
for(i in 1:length(cebo)){
  tempr=r[[which(gsub("\\."," ",names(r)) %in% cebo[[i]])]]
  if(dim(tempr)[3]>1){lr[[i]]=sum(tempr)}
  else{lr[[i]]=tempr}
}
slr=lapply(lr,function(x){x/sum(matrix(x))})
ter=stack(slr)
fear=sum(ter)
sfear=rastToIso(fear)
sfear=1-sfear
plot(sfear,col=rev(heat.colors(100)))

#SAVE FINAL RASTER TO USE IT IN STEP SELECTION FUNCTION
writeRaster(sfear,"wolf_pressure.tif",overwrite=TRUE)