#Diminution temps 
# analyses valvo TAOS M2 Lucie
# library
library(dplyr) # tidyverse beuuarggg
library(tidyr) 
library(data.table) # quick
library(parallel)
library(ggplot2) # graph
library(arrow) # speed up filtering and reading
library(lubridate)
library(scattermore) # plot high density data
library(RcppRoll) # quick rolling function
library(stringr)
library(ggpubr)



# Localisation dossier parquet avec données normées
#parquetinsitu<-open_dataset("C:/Users/lbourdon/Desktop/Données valvo in situ/parquet/insitunorm")
parquetinsitu<-open_dataset("//datawork/datawork-taos-s/intranet/valvo/Data/parquet/insitunorm")

# Mise en place d'un objet avec données base
metadat<-parquetinsitu%>%
  dplyr::select(TagID,spp,idvalvo,id)%>%
  distinct()%>%
  collect()

## 2) Boucle calcul indicateurs statistiques 10 minutes-----------------------------
indicateurs<-data.frame()# création objet vide dans lequel stock toutes les données de la boucle

#Boucle 
for(i in 1:nrow(metadat)){ 
  # i=14
  print(metadat$TagID[i])
  op <- options(digits.secs=3)
  dat<-parquetinsitu%>%
    dplyr::filter(TagID==metadat$TagID[i],id==metadat$id[i],is.finite(magX))%>%
    dplyr::select(TagID,id,Timestamp,magX,norm,mean,max,min)%>%collect()%>%
    dplyr::mutate(time=dmy_hms(Timestamp))
  options(op)
  dat <- dat[-1,] # on sup la premiere ligne qui contient NA suite lissage pour pouvoir ensuite cumsum
  dat<-dat%>%
    transmute(TagID,id,DateTime=time,norm,max,min)
  # Mise en place seuil fermeture => quantile 85% (85% des valeurs mesurées sont en dessous et 5 au dessus)
  s1 <- quantile(dat$norm,0.85,na.rm=T)
  
  #ggplot(dat, aes(x= DateTime, y= norm))+geom_point()+ggtitle(print(metadat$TagID[i]))+geom_scattermore()
  dat <- dat%>%
    mutate_at(.funs = list(prec = lag),
              .vars = vars(norm))
  # TEMPS OUV/FERM
  dat <- dat%>%
    mutate(close=ifelse(norm>s1,1,0)) # Création colonne close = si ind considéré comme fermé alors = 1, sinon =0
  dat <- dat%>%
    mutate(open=ifelse(norm<s1,1,0)) # Création colonne open = si ind considéré comme ouvert alors = 1, sinon =v0
  
  # NB OUV/FERM
  dat <- dat%>%
    mutate(close2=ifelse(norm>s1 & prec<s1,1,0)) # Création colonne close = si ind considéré comme fermé alors = 1, sinon =0
  dat <- dat%>%
    mutate(open2=ifelse(norm<s1 & prec>s1,1,0)) 
  
  #AMPLITUDE
  #dat <- dat%>%
  #mutate(ampli=(dat$max-dat$min)) # Création colonne amplitude exprimée en %
  
  #Indicateurs stats sur 10 min
  indicateurstmp <- dat%>%
    mutate(heure=substr(DateTime,1,15))%>% # Réduit le temps au 15 ème caractère = prend uniquement la date et la dizaine des minutes
    group_by(id,heure)%>% 
    summarise(closenumb=sum(close2,na.rm=T),# closenumb = Nombre de fermeture par 10 min
              opennumb = sum(open2,na.rm=T),# opennumb = Nombre d'ouverture par 10 min
              closetime=sum(close,na.rm=T)/2, # closetime = temps(en sec) ind fermé par 10 min
              opentime = 600-closetime, # opentime = temps(en sec) ind ouvert par 10 min
              isma=sum(abs(diff(norm)),na.rm=T),# indicateur isma = somme cumulée norm par 10 min
              #ampli =mean(ampli,na.rm=T),
              #ampli2=ampli*100/max(dat$ampli),
              mean2= mean(DateTime)) %>% 
    ungroup()
  indicateurs <- rbind(indicateurs, indicateurstmp) #Stock les données dans objet vide
}

#convert heure in time
pipo <- indicateurs %>% filter(nchar(heure) <= 10) %>% mutate(heure= paste0(heure, " 00:05:00"))
popo <- indicateurs %>% filter(nchar(heure) >= 11) %>% mutate(heure= paste0(heure, "5:00"))
ind2 <- rbind(pipo, popo)

# indicateurs <- indicateurs%>%select(-mean2)
# summary(indicateurs)
ind2$heure <- ymd_hms(ind2$heure)
ind2$jour <- substr(as.character(ind2$heure),1,10)
# ici, on vire les jours au cous desquels les rotation de valvo ont été effectuées
ind2<-ind2%>%filter(jour!="2023-02-22",jour!="2023-03-21",jour!="2023-04-19",jour!="2023-05-04")

# Graphique
graphics.off()
X11()
p1 <- ggplot(ind2,aes(x=heure,y=id,fill=log(opennumb+1)))+
  geom_raster()+
  scale_fill_gradientn(colors = hcl.colors(10, palette = "YlGnBu", alpha = NULL, rev = FALSE, fixup = TRUE))+ 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b")+
  ggtitle("Valvogramme du Nb. d'ouvertures par tranche de 10 minutes", subtitle = "Du 22/02/2023 au 17/05/2023") +
  xlab("Date") + ylab("Individu")
p2 <- ggplot(ind2,aes(x=heure,y=id,fill=closetime))+
  geom_raster()+
  scale_fill_gradientn(colors = hcl.colors(10, palette = "YlGnBu", alpha = NULL, rev = FALSE, fixup = TRUE))+ 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b")+
  ggtitle("Valvogramme du Temps de fermeture (secondes) par tranche de 10 minutes", subtitle = "Du 22/02/2023 au 17/05/2023") +
  xlab("Date") + ylab("Individu")
p3 <- ggplot(ind2,aes(x=heure,y=id,fill=log(isma+1)))+
  geom_raster()+
  scale_fill_gradientn(colors = hcl.colors(10, palette = "YlGnBu", alpha = NULL, rev = FALSE, fixup = TRUE))+ 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b")+
  ggtitle("Valvogramme de la Somme des écarts cumulés par tranche de 10 minutes", subtitle = "Du 22/02/2023 au 17/05/2023") +
  xlab("Date") + ylab("Individu")

ggsave(p1, file="//datawork/datawork-taos-s/intranet/valvo/outputs/valvogrammes_230726/p1.png",
       width = 8, height = 4, dpi = 300, units = "in", device='png')
ggsave(p2, file="//datawork/datawork-taos-s/intranet/valvo/outputs/valvogrammes_230726/p2.png",
       width = 8, height = 4, dpi = 300, units = "in", device='png')
ggsave(p3, file="//datawork/datawork-taos-s/intranet/valvo/outputs/valvogrammes_230726/p3.png",
       width = 8, height = 4, dpi = 300, units = "in", device='png')

write.csv2(ind2, file="//datawork/datawork-taos-s/intranet/valvo/Data/Data-indicateurs/indicateurs10min_230706.csv")

## 3) Boucle calcul indicateurs statistiques 1 heure-----------------------------
indicateurs<-data.frame()# création objet vide dans lequel stock toutes les données de la boucle
#Boucle 
for(i in 1:nrow(metadat)){ 
  # i=14
  print(metadat$TagID[i])
  op <- options(digits.secs=3)
  dat<-parquetinsitu%>%
    dplyr::filter(TagID==metadat$TagID[i],id==metadat$id[i],is.finite(magX))%>%
    dplyr::select(TagID,id,Timestamp,magX,norm,mean,max,min)%>%collect()%>%
    dplyr::mutate(time=dmy_hms(Timestamp))
  options(op)
  dat <- dat[-1,] # on sup la premiere ligne qui contient NA suite lissage pour pouvoir ensuite cumsum
  dat<-dat%>%
    transmute(TagID,id,DateTime=time,norm,max,min)
  # Mise en place seuil fermeture => quantile 85% (85% des valeurs mesurées sont en dessous et 5 au dessus)
  s1 <- quantile(dat$norm,0.85,na.rm=T)
  
  #ggplot(dat, aes(x= DateTime, y= norm))+geom_point()+ggtitle(print(metadat$TagID[i]))+geom_scattermore()
  dat <- dat%>%
    mutate_at(.funs = list(prec = lag),
              .vars = vars(norm))
  # TEMPS OUV/FERM
  dat <- dat%>%
    mutate(close=ifelse(norm>s1,1,0)) # Création colonne close = si ind considéré comme fermé alors = 1, sinon =0
  dat <- dat%>%
    mutate(open=ifelse(norm<s1,1,0)) # Création colonne open = si ind considéré comme ouvert alors = 1, sinon =v0
  
  # NB OUV/FERM
  dat <- dat%>%
    mutate(close2=ifelse(norm>s1 & prec<s1,1,0)) # Création colonne close = si ind considéré comme fermé alors = 1, sinon =0
  dat <- dat%>%
    mutate(open2=ifelse(norm<s1 & prec>s1,1,0)) 
  
  #AMPLITUDE
  #dat <- dat%>%
  #mutate(ampli=(dat$max-dat$min)) # Création colonne amplitude exprimée en %
  #Indicateurs stats sur 1 heure
  indicateurstmp <- dat%>%
    mutate(heure=substr(DateTime,1,13))%>% # Réduit le temps au 15 ème caractère = prend uniquement la date et la dizaine des minutes
    group_by(id,heure)%>% 
    summarise(closenumb=sum(close2,na.rm=T),# closenumb = Nombre de fermeture par heure
              opennumb = sum(open2,na.rm=T),# opennumb = Nombre d'ouverture par heure
              closetime=sum(close,na.rm=T)/2, # closetime = temps(en sec) ind fermé par heure
              opentime = 3600-closetime, # opentime = temps(en sec) ind ouvert par heure
              isma=sum(abs(diff(norm)),na.rm=T),# indicateur isma = somme cumulée norm par heure
              #ampli =mean(ampli,na.rm=T),
              #ampli2=ampli*100/max(dat$ampli),
              mean2= mean(DateTime)) %>% 
    ungroup()
  indicateurs <- rbind(indicateurs, indicateurstmp) #Stock les données dans objet vide
}


#convert heure in time
pipo <- indicateurs %>% filter(nchar(heure) <= 10) %>% mutate(heure= paste0(heure, " 00:00:00"))
popo <- indicateurs %>% filter(nchar(heure) >= 11) %>% mutate(heure= paste0(heure, ":00:00"))
ind3 <- rbind(pipo, popo)

# indicateurs <- indicateurs%>%select(-mean2)
# summary(indicateurs)
ind3$heure <- ymd_hms(ind3$heure)
ind3 <- ind3 %>% arrange(id, heure)
ind3$jour <- substr(as.character(ind3$heure),1,10)
# ici, on vire les jours au cous desquels les rotation de valvo ont été effectuées
ind3<-ind3%>%filter(jour!="2023-02-22",jour!="2023-03-21",jour!="2023-04-19",jour!="2023-05-04")

# Graphique
graphics.off()
X11()
p4 <- ggplot(ind3,aes(x=heure,y=id,fill=log(opennumb+1)))+
  geom_raster()+
  scale_fill_gradientn(colors = hcl.colors(10, palette = "YlGnBu", alpha = NULL, rev = FALSE, fixup = TRUE))+ 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b")+
  ggtitle("Valvogramme du Nb. d'ouvertures par tranche d'une heure", subtitle = "Du 22/02/2023 au 17/05/2023") +
  xlab("Date") + ylab("Individu")
p5 <- ggplot(ind3,aes(x=heure,y=id,fill=closetime))+
  geom_raster()+
  scale_fill_gradientn(colors = hcl.colors(10, palette = "YlGnBu", alpha = NULL, rev = FALSE, fixup = TRUE))+ 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b")+
  ggtitle("Valvogramme du Temps de fermeture (secondes) par tranche d'une heure", subtitle = "Du 22/02/2023 au 17/05/2023") +
  xlab("Date") + ylab("Individu")
p6 <- ggplot(ind3,aes(x=heure,y=id,fill=log(isma+1)))+
  geom_raster()+
  scale_fill_gradientn(colors = hcl.colors(10, palette = "YlGnBu", alpha = NULL, rev = FALSE, fixup = TRUE))+ 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b")+
  ggtitle("Valvogramme de la Somme des écarts cumulés par tranche d'une heure", subtitle = "Du 22/02/2023 au 17/05/2023") +
  xlab("Date") + ylab("Individu")

ggsave(p4, file="//datawork/datawork-taos-s/intranet/valvo/outputs/valvogrammes_230726/p4.png",
       width = 8, height = 4, dpi = 300, units = "in", device='png')
ggsave(p5, file="//datawork/datawork-taos-s/intranet/valvo/outputs/valvogrammes_230726/p5.png",
       width = 8, height = 4, dpi = 300, units = "in", device='png')
ggsave(p6, file="//datawork/datawork-taos-s/intranet/valvo/outputs/valvogrammes_230726/p6.png",
       width = 8, height = 4, dpi = 300, units = "in", device='png')

write.csv2(ind3, file="//datawork/datawork-taos-s/intranet/valvo/Data/Data-indicateurs/indicateurs1heure_230706.csv")