#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain fig S 1 of the supplementary material of the manuscript submitted to PeerJ Gauger et al 2021 Diel influences on bottlenose dolphin acoustic detection in a coastal lagoon in the southwestern Gulf of California
#####################################################
ifelse("dplyr" %in% rownames(installed.packages()) == FALSE, 
       install.packages("dplyr", dependencies = T),
       ("dplyr is already installed"))
ifelse("suncalc" %in% rownames(installed.packages()) == FALSE, 
       install.packages("suncalc", dependencies = T),
       ("suncalc is already installed"))
ifelse("scales" %in% rownames(installed.packages()) == FALSE, 
       install.packages("scales", dependencies = T),
       ("scales is already installed"))
require(dplyr)
require(suncalc)
require(scales)

########
Sys.setlocale("LC_ALL","English")
Sys.setenv(TZ='UTC')
###################

analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
analysed_data_hour$date <-as.Date(analysed_data_hour$date,format="%d.%m.%Y")
analysed_data_hour$datehour_lt <- as.POSIXlt(paste0(analysed_data_hour$datehour,":00"), format="%d.%m.%Y %H:%M:%S", tz="MST")
analysed_data_hour$datehour <- as.POSIXct(paste0(analysed_data_hour$datehour,":00"), format="%d.%m.%Y %H:%M:%S", tz="GMT")

table(analysed_data_hour$lunar_phase_2)
table(analysed_data_hour$cluster_hcpc,analysed_data_hour$lunar_phase_2)
table(analysed_data_hour$month, analysed_data_hour$cluster_hcpc)

marea01 <- read.csv("Supplemental Data S2.csv", sep=";")
marea01$datehour <- as.POSIXct(marea01$datehour, format="%d.%m.%Y %H:%M:%S", tz="GMT")
marea01$datehour_lt <- as.POSIXlt(marea01$datehour, format="%d.%m.%Y %H:%M:%S", tz="MST")
marea01$Date_lt <-as.Date(substr(marea01$datehour_lt,1,10))
marea01$Date <-as.Date(marea01$datehour)
MAX = max(marea01$datehour_lt)

SS <- suncalc::getSunlightTimes(date = as.Date(unique(marea01$Date)), lat = 24, lon = -110, data = NULL,keep = c("sunrise", "sunset"), tz = "UTC")

K = aggregate(analysed_data_hour$cluster_hcpc, list(analysed_data_hour$cluster_hcpc, analysed_data_hour$month_campaign, analysed_data_hour$date),dplyr::first)
K2=aggregate(K$Group.3, list(K$Group.2), first)  
K2
# dates between deployments
sort(diff.Date(sort(K2$x)))

##### ######

  unique(analysed_data_hour$month_campaign)
#x11()  
png("Figure S1.png", width = 150, height = 150, units = "mm", pointsize = 12, res=300)
    par(mfcol=c(8,3), oma=c(4,4,3,1), mar=c(0,1,0,1))
    ##### #####
    K = aggregate(analysed_data_hour$cluster_hcpc, list(analysed_data_hour$cluster_hcpc, analysed_data_hour$month_campaign, analysed_data_hour$date),dplyr::first)
    names(K) <- c("cluster_hcpc","month_campaign","date","x")
    K$Station <- as.numeric(aggregate(analysed_data_hour$Station, list(analysed_data_hour$cluster_hcpc, analysed_data_hour$month_campaign, analysed_data_hour$date),dplyr::first)$x)
    K$lunar_phase <- (aggregate(analysed_data_hour$lunar_phase_2, list(analysed_data_hour$cluster_hcpc, analysed_data_hour$month_campaign, analysed_data_hour$date),dplyr::first)$x)
    K$month_campaign <- (aggregate(analysed_data_hour$month_campaign, list(analysed_data_hour$cluster_hcpc, analysed_data_hour$month_campaign, analysed_data_hour$date),dplyr::first)$x)
    K$date<-as.Date(K$date)
    K <- K[order(K$cluster_hcpc, K$date),]
    ll=1
    C = sort(unique(K$month_campaign[K$cluster_hcpc == ll]))
    L = K$date[K$cluster_hcpc == ll]
    LL= aggregate(K$date[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLp= aggregate(K$Station[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLm= aggregate(K$lunar_phase[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLc= aggregate(K$month_campaign[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    
    for(k in 1:6){
      LL= K$date[K$month_campaign == C[k]]
      
      plot(marea01$datehour_lt[marea01$Date_lt == L[k]], marea01$tide[marea01$Date_lt == L[k]] , col="white", ylab="", xlab="", yaxt="n", xaxt="n", ylim=c(-1,1),cex=.01)
      if(k%in% c(1)) mtext(side=3, "Aug-Oct", line=.5)
      axis(side=2, las=2, labels=F)
      if(k==8) axis.POSIXct(side=1,at=seq(marea01$datehour_lt[marea01$Date_lt == LL[1]][1],marea01$datehour_lt[marea01$Date_lt == LL[1]][139],by="hours"))
      for(l in 1:length(LL)){
        par(new=T)
        plot(marea01$datehour_lt[marea01$Date_lt == LL[l]], marea01$tide[marea01$Date_lt == LL[l]], ylab="", xlab="", type="l", xaxt="n", yaxt="n", ylim=c(-1,1), col=alpha("white",.01))
        rect(0,-1,SS$sunrise[SS$date== LL[l]],1, border=alpha("grey",1/length(LL)),col=alpha("grey",1/length(LL)))
        rect(SS$sunset[SS$date== LL[l]],-1,MAX,1, border=alpha("grey",1/length(LL)), col=alpha("grey",1/length(LL)))
        lines(marea01$datehour_lt[marea01$datehour%in% analysed_data_hour$datehour & marea01$Date_lt == LL[l]],  marea01$tide[marea01$datehour%in% analysed_data_hour$datehour & marea01$Date_lt == LL[l]])
        #points(analysed_data_hour$datehour_lt[analysed_data_hour$dp10m_ELAP1 > 0], analysed_data_hour$tide[analysed_data_hour$dp10m_ELAP1 > 0], cex= analysed_data_hour$dp10m_ELAP1[analysed_data_hour$dp10m_ELAP1 > 0]/10+1, pch=20, bg="white")
      }
      mtext(side=1, LL[1], line=-1, adj = .995, cex=.6)
      if(k<5) mtext(side=3, "*", line=-2.5, adj = .995, cex=1.6)
      mtext(side=3, paste(LLc[k]), line=-1, adj = .005, cex=.6)
      mtext(side=1, paste("station",LLp[k]), line=-2, adj = .995, cex=.6)
      mtext(side=3, line=-1.0,paste(length(analysed_data_hour$tide[analysed_data_hour$month_campaign %in% C[k]]),"hs"),adj = .995, cex=.6)
      par(new=T)
      plot(c(0,1),c(0,1), xaxt="n", yaxt="n", ylab="", xlab="", col=alpha("white",.01))
      points(0.05,0.1, pch=21, bg=1, cex=3)
      if(LLm[k] == "Waning")  points(0.06,0.1, pch=19, col="grey", cex=2.75)
      if(LLm[k] == "Waxing")  points(0.04,0.1, pch=19, col="grey", cex=2.75)
      if(LLm[k] == "New")   points(0.05,0.1, pch=21, bg="white", cex=3)
    }  
    par(new=T)
    plot(marea01$datehour_lt[marea01$Date_lt == L[k]], marea01$tide[marea01$Date_lt == L[k]], col=alpha("white",.01), ylab="", xlab="", yaxt="n", xaxt="n", ylim=c(-1,1))
    axis.POSIXct(side=1,at=seq(marea01$datehour_lt[marea01$Date_lt == L[k]][1],marea01$datehour_lt[marea01$Date_lt == L[k]][139],by="hours"))
    plot(c(0,1),c(0,1), axes=F, xlab="", ylab="", col="white")
    plot(c(0,1),c(0,1), axes=F, xlab="", ylab="", col="white")
    text(.2,.75,"moon phase", adj=0)  points(0.4,0.35, pch=19, col=1, cex=2)
    
    points(0.1,0.65, pch=19, col=1, cex=2)
    points(0.1,0.45, pch=19, col=1, cex=2)
    points(0.09,0.45, pch=19, col="white", cex=1.75)
    points(0.1,0.25, pch=21, bg="white", cex=2)
    points(0.1,0.05, pch=19, col=1, cex=2)
    points(0.11,0.05, pch=19, col="white", cex=1.75)
    
    text(.2,.65,"New ", adj=0)  
    text(.2,.45,"Waxing", adj=0)  
    text(.2,.25,"Full ", adj=0)  
    text(.2,.05,"Waning", adj=0)  
    
    mtext(side=2, outer=T, "tide [m]", line=2.5)
    
    ll=2
    C = sort(unique(K$month_campaign[K$cluster_hcpc == ll]))
    L = K$date[K$cluster_hcpc == ll]
    LL= aggregate(K$date[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLp= aggregate(K$Station[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLm= aggregate(K$lunar_phase[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLc= aggregate(K$month_campaign[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    for(k in c(5:8,1:4)){
      LL= K$date[K$month_campaign == C[k]]
      
      plot(marea01$datehour_lt[marea01$Date_lt == L[k]], marea01$tide[marea01$Date_lt == L[k]] , col="white", ylab="", xlab="", yaxt="n", xaxt="n", ylim=c(-1,1),cex=.01)
      if(k%in% c(5)) mtext(side=3, "Oct-Mar", line=.5)
      axis(side=2, las=2, labels=F)
      if(k==8) axis.POSIXct(side=1,at=seq(marea01$datehour_lt[marea01$Date_lt == LL[1]][1],marea01$datehour_lt[marea01$Date_lt == LL[1]][139],by="hours"))
      for(l in 1:length(LL)){
        par(new=T)
        plot(marea01$datehour_lt[marea01$Date_lt == LL[l]], marea01$tide[marea01$Date_lt == LL[l]], ylab="", xlab="", type="l", xaxt="n", yaxt="n", ylim=c(-1,1), col=alpha("white",.01))
        rect(0,-1,SS$sunrise[SS$date== LL[l]],1, border=alpha("grey",1/length(LL)),col=alpha("grey",1/length(LL)))
        rect(SS$sunset[SS$date== LL[l]],-1,MAX,1, border=alpha("grey",1/length(LL)), col=alpha("grey",1/length(LL)))
        lines(marea01$datehour_lt[marea01$datehour%in% analysed_data_hour$datehour & marea01$Date_lt == LL[l]],  marea01$tide[marea01$datehour%in% analysed_data_hour$datehour & marea01$Date_lt == LL[l]])
        #points(analysed_data_hour$datehour_lt[analysed_data_hour$dp10m_ELAP1 > 0], analysed_data_hour$tide[analysed_data_hour$dp10m_ELAP1 > 0], cex= analysed_data_hour$dp10m_ELAP1[analysed_data_hour$dp10m_ELAP1 > 0]/10+1, pch=20, bg="white")
      }
      mtext(side=1, LL[1], line=-1, adj = .995, cex=.6)
      if(k<5) mtext(side=3, "*", line=-2.5, adj = .995, cex=1.6)
      mtext(side=3, paste(LLc[k]), line=-1, adj = .005, cex=.6)
      mtext(side=1, paste("station",LLp[k]), line=-2, adj = .995, cex=.6)
      mtext(side=3, line=-1.0,paste(length(analysed_data_hour$tide[analysed_data_hour$month_campaign %in% C[k]]),"hs"),adj = .995, cex=.6)
      par(new=T)
      plot(c(0,1),c(0,1), xaxt="n", yaxt="n", ylab="", xlab="", col=alpha("white",.01))
      points(0.05,0.1, pch=21, bg=1, cex=3)
      if(LLm[k] == "Waning")  points(0.06,0.1, pch=19, col="grey", cex=2.75)
      if(LLm[k] == "Waxing")  points(0.04,0.1, pch=19, col="grey", cex=2.75)
      if(LLm[k] == "New")   points(0.05,0.1, pch=21, bg="white", cex=3)
    }  
    par(new=T)
    plot(marea01$datehour_lt[marea01$Date_lt == L[k]], marea01$tide[marea01$Date_lt == L[k]], col=alpha("white",.01), ylab="", xlab="", yaxt="n", xaxt="n", ylim=c(-1,1))
    axis.POSIXct(side=1,at=seq(marea01$datehour_lt[marea01$Date_lt == L[k]][1],marea01$datehour_lt[marea01$Date_lt == L[k]][139],by="hours"))
    mtext(side=1, line=2.5, outer=T, "hour")
   
    
    ll=3
    C = sort(unique(K$month_campaign[K$cluster_hcpc == ll]))
    L = K$date[K$cluster_hcpc == ll]
    LL= aggregate(K$date[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLp= aggregate(K$Station[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLm= aggregate(K$lunar_phase[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    LLc= aggregate(K$month_campaign[K$month_campaign%in%C], list(K$month_campaign[K$month_campaign%in%C]), dplyr::first)$x
    for(k in 1:7){
      LL= K$date[K$month_campaign == C[k]]
      
      plot(marea01$datehour_lt[marea01$Date_lt == L[k]], marea01$tide[marea01$Date_lt == L[k]], col="white", ylab="", xlab="", yaxt="n", xaxt="n", ylim=c(-1,1),cex=.01)
      if(k%in% c(1)) mtext(side=3, "May-Jul", line=.5)
      axis(side=2, las=2, labels=F)
      for(l in 1:length(LL)){
        par(new=T)
        plot(marea01$datehour_lt[marea01$Date_lt == LL[l]], marea01$tide[marea01$Date_lt == LL[l]], ylab="", xlab="", type="l", xaxt="n", yaxt="n", ylim=c(-1,1), col=alpha("black",1))
        rect(0,-1,SS$sunrise[SS$date== LL[l]],1, border=alpha("grey",1/length(LL)),col=alpha("grey",1/length(LL)))
        rect(SS$sunset[SS$date== LL[l]],-1,MAX,1, border=alpha("grey",1/length(LL)), col=alpha("grey",1/length(LL)))
        lines(marea01$datehour_lt[marea01$datehour%in% analysed_data_hour$datehour & marea01$Date_lt == LL[l]],  marea01$tide[marea01$datehour%in% analysed_data_hour$datehour & marea01$Date_lt == LL[l]])
        #points(analysed_data_hour$datehour_lt[analysed_data_hour$dp10m_ELAP1 > 0], analysed_data_hour$tide[analysed_data_hour$dp10m_ELAP1 > 0], cex= analysed_data_hour$dp10m_ELAP1[analysed_data_hour$dp10m_ELAP1 > 0]/10+1, pch=20, bg="white")
      }
      mtext(side=1, LL[1], line=-1, adj = .995, cex=.6)
      mtext(side=3, line=-1.0,paste(length(analysed_data_hour$tide[analysed_data_hour$month_campaign %in% C[k]]),"hs"),adj = .995, cex=.6)
      mtext(side=3, paste(LLc[k]), line=-1, adj = .005, cex=.6)
      mtext(side=1, paste("station",LLp[k]), line=-2, adj = .995, cex=.6)
      par(new=T)
      plot(c(0,1),c(0,1), xaxt="n", yaxt="n", ylab="", xlab="", col=alpha("white",.01))
      points(0.05,0.1, pch=21, bg=1, cex=3)
      if(LLm[k] == "Waning")  points(0.06,0.1, pch=19, col="grey", cex=2.75)
      if(LLm[k] == "Waxing")  points(0.04,0.1, pch=19, col="grey", cex=2.75)
      if(LLm[k] == "New")   points(0.05,0.1, pch=21, bg="white", cex=3)
    }  
    par(new=T)
    plot(marea01$datehour_lt[marea01$Date_lt == L[k]], marea01$tide[marea01$Date_lt == L[k]], col=alpha("white",.01), ylab="", xlab="", yaxt="n", xaxt="n", ylim=c(-1,1))
    axis.POSIXct(side=1,at=seq(marea01$datehour_lt[marea01$Date_lt == L[k]][1],marea01$datehour_lt[marea01$Date_lt == L[k]][139],by="hours"))
    
    
    
    
    
    ##### #####
    
  dev.off()