#################################################################
#R script to do an automatic classification of flying and non flying
#behavior based on the acceleration data for Eidolon helvum 
#studies with eobs-tags
##################################################################              

library(maptools); library(move)
## download the gps and the acceleration data separately from movebank per study ##
filenameACC<-"Eidolon helvum; Ghana, Kibi; dry season 2013 [eobs]_acc.csv" #acc file name
filenameGPS<-"Eidolon helvum; Ghana, Kibi; dry season 2013 [eobs]_gps.csv" #gps file name
outACC<-"Eidolon_Ghana_Kibi_dry_2013_acc&behav.csv" #output acc file name
outGPS<-"Eidolon_Ghana_Kibi_dry_2013_gps&behav.csv" #output gps file name

accALL<-read.csv(paste0(getwd(),"/",filenameACC)) #read in the acceleration data
accALL$timestamp<-as.POSIXct(as.character(accALL$timestamp),
                             format="%Y-%m-%d %H:%M:%OS",tz="UTC")
myindACC<-unique(accALL$individual.local.identifier)
myindACC

gpsMove<-move(paste0(getwd(),"/",filenameGPS)) # read the gps file as moveStack object
gpsMoveL<-split(gpsMove) # if several individuals, convert the moveStack into a list 
myindGPS<-names(gpsMoveL)
myindGPS #check if individuals names are the same in acc file as in move object

AllInd<-list()
for(ind in myindACC){
  print(ind)
  gps <- gpsMoveL[[ind]]
  acc <- accALL[accALL$individual.local.identifier==ind,]
  
  ##creating table with the samples per acc burst
  accaxes <- strsplit(as.character(acc$eobs.accelerations.raw),                                                    		split=" ", fixed=T) 
  numberSamples<-length(unlist(accaxes[1]))
  accaxes2 <- matrix(unlist(accaxes), ncol=numberSamples, byrow=T) 
  colnames(accaxes2) <- paste0("V", 1:numberSamples)
  
  acc2 <-data.frame(acc$timestamp)
  colnames(acc2)<-"timestamp"
  acc2$timestamp<-as.POSIXct(acc2$timestamp, 
                             format="%Y-%m-%d %H:%M:%OS", tz="UTC")
  acc2$var.y <- NA
  acc2$var.z <- NA
  for(j in 1:nrow(accaxes2)){
    yax<-as.numeric(accaxes2[j,seq(2,numberSamples,3)])
    zax<-as.numeric(accaxes2[j,seq(3,numberSamples,3)])
    acc2$var.y[j] <- var(yax)
    acc2$var.z[j] <- var(zax) 
  }
  
  ## clustering via kmeans the variance of the z and y axis
  myk<-acc2[,c("var.y","var.z")]
  mykmeans <- kmeans(myk,2)
  acc2$cluster<-mykmeans$cluster
  
  ## determine which cluster is flying or not-flying
  agmn<-aggregate(acc2$var.z,by=list(acc2$cluster),FUN=min) 
  acc2$activity<-NA
  
  if(agmn$x[agmn$Group.1==1]>agmn$x[agmn$Group.1==2])
         {acc2$activity[acc2$cluster==1]<-"Flying" ; 
         acc2$activity[acc2$cluster==2]<-"NotFlying"} 
  
  if(agmn$x[agmn$Group.1==2]>agmn$x[agmn$Group.1==1])
         {acc2$activity[acc2$cluster==2]<-"Flying" ; 
        acc2$activity[acc2$cluster==1]<-"NotFlying"} 
  
  ## assigning day or night to each busts (according to sunset and sunrise)
  centroidLatLong<-centroid(gps@coords)
  acc2$long <- centroidLatLong[1]
  acc2$lat <- centroidLatLong[2]
  myindSP <- as.matrix(acc2[,c("long","lat")])
  datetime <- as.POSIXct(acc2$timestamp, 
                         format="%Y-%m-%d %H:%M:%OS", tz="UTC")
  sunrise<-sunriset(myindSP, datetime,
                    proj4string=CRS("+proj=longlat +ellps=WGS84 	
                                    +datum=WGS84 +towgs84=0,0,0"), 
                    direction=c("sunrise"), POSIXct.out=T)
  acc2$sunrise <- sunrise$time
  sunset<-sunriset(myindSP, datetime,
                   proj4string= CRS("+proj=longlat +ellps=WGS84 	
                                    +datum=WGS84 +towgs84=0,0,0"), 
                   direction=c("sunset"), POSIXct.out=T)
  acc2$sunset <- sunset$time
  acc2$DayTime <- ifelse(acc2$timestamp>acc2$sunrise & 
                           acc2$timestamp<acc2$sunset,"day","night")
  
  ## assigning each gps point with the acc activity classification nearest in time
  gps$activity<-NA
  gps$DayTime<-NA
  for (h in 1:nrow(gps@data)){ 
    tDiff<-abs(acc2$timestamp-gps@timestamps[h]) ## searches for the nearst point in time of the acc to the gps point
    minN<-which.min(tDiff) 
    gps[h,c("activity","DayTime")] <- acc2[minN,c("activity","DayTime")] 
    gps$tdiff[h]<- min(as.difftime(tDiff, units="secs")) ## time difference between the selected acc burst and gps point, in seconds
  }
  gps$activity<-factor(gps$activity,level=c("Flying","NotFlying"))
  gps$DayTime<-factor(gps$DayTime,levels=c("night","day"))
  
  AllInd[ind]<-gps ## creating a list with all inidividuals of the study
  
  ## create csv files with gps fixes for individuals of the study
  myTable <- data.frame(cbind(as.POSIXct(gps@timestamps,
                                         format="%Y-%m-%d %H:%M:%OS",tz="UTC"),gps@coords,
                              as.character(gps$activity), 
                              as.character(gps$DayTime),ind))
  colnames(myTable) <- c("timestamp","longitude","latitude",
                         "activity","DayTime","individuals")
  myTable$timestamp<-as.POSIXct(gps@timestamps, 
                                format="%Y-%m-%d %H:%M:%OS", tz="UTC")
  if(ind==myindACC[1]){GpsAllIndv<-myTable}else{GpsAllIndv<-rbind(GpsAllIndv,myTable)}
  
  ## adds to the original acc file the behavioral classification of each burst
  acc$behavioural.classification<-acc2$activity
  if(ind==myindACC[1]){AccAllInd <- acc}else{AccAllInd <-rbind(AccAllInd,acc)}
}
write.csv(AccAllInd,file=paste0(getwd(),"/",outACC)) ## save csv file of raw acc data with behavioral classification 
write.csv(GpsAllIndv,file=paste0(getwd(),"/",outGPS)) ## save csv file of gps points with behavioral classification

####################
### example plots ###
####################
AllMove <- moveStack(AllInd) #convert list to moveStack

plot(AllMove[["X2402"]],col=c("blue","red")[AllMove[["X2402"]]$activity],pch=19) #blue=flying; red=NotFlying
plot(AllMove[["X2402"]][AllMove[["X2402"]]$activity=="NotFlying"], pch=19, col="red") # plot only NotFlying points
plot(AllMove[["X2402"]][AllMove[["X2402"]]$activity=="NotFlying" & AllMove[["X2402"]]$DayTime=="day"], pch=19, col="red") #plot only NotFlying points during the day

#################################################################
### check in detail how well the classification Flying/NotFlying is doing  ###
##################################################################
accALL<-read.csv(paste0(getwd(),"/Eidolon helvum; Ghana, Kibi; wet season 2012 [eobs]_acc.csv")) #read in the acceleration data
myindACC<-unique(accALL$individual.local.identifier)
acc <- accALL[accALL$individual.local.identifier==myindACC[1],] #example for individual 2394

accaxes <- strsplit(as.character(acc$eobs.accelerations.raw), split=" ", fixed=T) #creating table with the samples per acc burst
numberSamples<-length(unlist(accaxes[1]))
accaxes2 <- matrix(unlist(accaxes), ncol=numberSamples, byrow=T) 
colnames(accaxes2) <- paste0("V", 1:numberSamples)

acc2 <-data.frame(acc$timestamp)
colnames(acc2)<-"timestamp"
acc2$timestamp<-as.POSIXct(acc2$timestamp, format="%Y-%m-%d %H:%M:%OS", tz="UTC")
acc2$var.x <- NA  
acc2$var.y <- NA
acc2$var.z <- NA

for(j in 1:nrow(accaxes2)){
  xax<-as.numeric(accaxes2[j,seq(1,numberSamples,3)])
  yax<-as.numeric(accaxes2[j,seq(2,numberSamples,3)])
  zax<-as.numeric(accaxes2[j,seq(3,numberSamples,3)])
  if(j==1){xaxB<-xax}else{xaxB<-c(xaxB,xax)}
  if(j==1){yaxB<-yax}else{yaxB<-c(yaxB,yax)}
  if(j==1){zaxB<-zax}else{zaxB<-c(zaxB,zax)}
  acc2$var.x[j] <- var(xax)
  acc2$var.y[j] <- var(yax)
  acc2$var.z[j] <- var(zax) 
}
allaxis<-data.frame(cbind(xaxB,yaxB,zaxB))

#clustering via kmeans the variance of the z and y axis
myk<-acc2[,c("var.y","var.z")]
mykmeans <- kmeans(myk,2)
acc2$cluster<-mykmeans$cluster

#determine which cluster is flying or not flying
agmn<-aggregate(acc2$var.z,by=list(acc2$cluster),FUN=min) 
acc2$activity<-NA
if(agmn$x[agmn$Group.1==1]>agmn$x[agmn$Group.1==2])	{acc2$activity[acc2$cluster==1]<-"Flying" ; 	acc2$activity[acc2$cluster==2]<-"NotFlying"} 
if(agmn$x[agmn$Group.1==2]>agmn$x[agmn$Group.1==1])	{acc2$activity[acc2$cluster==2]<-"Flying" ; 	acc2$activity[acc2$cluster==1]<-"NotFlying"} 

allaxis$cluster<- as.factor(rep(c(acc2$activity),each=264))

### plots of the acc data as they appear in the acc viewer of movebank
# subset of the 1st 3h of data
zz<-allaxis[1:47520,] 
plot(zz$zaxB,type="l",col="blue") # z axis
lines(zz$yaxB,col="green") # y axis
lines(zz$xaxB,col="red")# x axis
points(zz$zaxB,col=c("black","yellow")[zz$cluster],pch=20,cex=.6) #black is Flying, yellow is NotFlying

# subset of 60 min of data somewhere in the middle
yy<-allaxis[174240:190080,] 
plot(yy$zaxB,type="l",col="blue")
lines(yy$yaxB,col="green")
lines(yy$xaxB,col="red")
points(yy$zaxB,col=c("black","yellow")[yy$cluster],pch=20,cex=.6) #black is Flying, yellow is NotFlying

#all data
plot(allaxis$zaxB,type="l",col="blue")
lines(allaxis$yaxB,col="green")
lines(allaxis$xaxB,col="red")
points(allaxis$zaxB,col=c("black","yellow")[allaxis$cluster],pch=20,cex=.6)