#####################################################
# Marco FW Gauger
# 2021/08/23
#
# This is the script to obtain fig 5 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("geosphere" %in% rownames(installed.packages()) == FALSE, 
       install.packages("geosphere", dependencies = T),
       ("geosphere is already installed"))
ifelse("dplyr" %in% rownames(installed.packages()) == FALSE, 
       install.packages("dplyr", dependencies = T),
       ("dplyr is already installed"))
ifelse("factoextra" %in% rownames(installed.packages()) == FALSE, 
       install.packages("factoextra", dependencies = T),
       ("factoextra is already installed"))
ifelse("FactoMineR" %in% rownames(installed.packages()) == FALSE, 
       install.packages("FactoMineR", dependencies = T),
       ("FactoMineR is already installed"))
ifelse("ggsci" %in% rownames(installed.packages()) == FALSE, 
       install.packages("ggsci", dependencies = T),
       ("ggsci is already installed"))
ifelse("cowplot" %in% rownames(installed.packages()) == FALSE, 
       install.packages("cowplot", dependencies = T),
       ("cowplot is already installed"))
ifelse("vegan" %in% rownames(installed.packages()) == FALSE, 
       install.packages("vegan", dependencies = T),
       ("vegan is already installed"))
ifelse("RColorBrewer" %in% rownames(installed.packages()) == FALSE, 
       install.packages("RColorBrewer", dependencies = T),
       ("RColorBrewer is already installed"))
ifelse("viridis" %in% rownames(installed.packages()) == FALSE, 
       install.packages("viridis", dependencies = T),
       ("viridis is already installed"))
ifelse("gplots" %in% rownames(installed.packages()) == FALSE, 
       install.packages("gplots", dependencies = T),
       ("gplots is already installed"))
ifelse("parameters" %in% rownames(installed.packages()) == FALSE, 
       install.packages(c("parameters","bayestestR")),
       ("parameters is already installed"))



require(gplots)
require(RColorBrewer)
require(vegan)
require(cowplot)
require(ggsci)
require(FactoMineR)
require(factoextra)
require(geosphere)
require(dplyr)



analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
names(analysed_data_hour)
analysed_data_hour$month_campaign[analysed_data_hour$month_campaign==c("11_13")] <- "11_12"
j <- aggregate(analysed_data_hour$dp10m_ELAP1, list(analysed_data_hour$month_campaign), sum)
names(j) <- c("campaign","dph")
j$d10pm <- 
  aggregate(analysed_data_hour$dp10m_ELAP1, list(analysed_data_hour$month_campaign), sum)$x
j$moon <- 
  aggregate(analysed_data_hour$lunar_phase , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$SST <- 
  aggregate(analysed_data_hour$SST , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$SST_BLAP <- 
  aggregate(analysed_data_hour$SST_BLAP , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$CHL <- 
  aggregate(analysed_data_hour$CHL_BLAP , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$depth <- 
  aggregate(analysed_data_hour$depth , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$sediment <- 
  aggregate(analysed_data_hour$sediment , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$sediment <- as.numeric(j$sediment)
j$distance_mangrove <- 
  aggregate(analysed_data_hour$distance_mangrove , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$distance_coast <- 
  aggregate(analysed_data_hour$distance_coast , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$distance_BLAP <- 
  aggregate(analysed_data_hour$distance_BLAP , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$delta_SST <- j$SST - j$SST_BLAP
j$daylength <- geosphere::daylength(24, as.Date(aggregate(analysed_data_hour$datehour , list(analysed_data_hour$month_campaign),  dplyr::first)$x, format="%d.%m.%Y"))


normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}

pca_j <- j
for(k in 2:length(j)){pca_j[,k] <- normalize(j[,k])}
names(pca_j)


parameters::check_kmo(pca_j[,c(5,7,14)]) # SST, CHL & daylength

rownames(pca_j) <- j$campaign

res.pca_j2 <- PCA(pca_j[,c(5,7,14)], ncp = 2, graph = FALSE)
# Compute hierarchical clustering on principal components
res.hcpc_j2 <- HCPC(res.pca_j2, graph = FALSE)

map1 <- fviz_dend(res.hcpc_j2, 
                  cex = 0.7,                     
                  # Label size
                  palette = c("#1B9E77","grey25","magenta"),               
                  # Color palette see ?ggpubr::ggpar
                  rect = TRUE, rect_fill = TRUE, 
                  # Add rectangle around groups
                  rect_border = c("#1B9E77","grey25","magenta"),           
                  # Rectangle color
                  labels_track_height = 0.95    ,
                  main = ""
                  # Augment the room for labels
)
map1 + scale_fill_manual(values=c("#1B9E77","grey25","magenta"))
(map2 <- fviz_cluster(res.hcpc_j2,
                      repel = TRUE,            # Avoid label overlapping
                      show.clust.cent = TRUE, # Show cluster centers
                      palette =c("#1B9E77","grey25","magenta"),         # Color palette see ?ggpubr                     ggtheme = theme_minimal(),
                      main = "",
                      labelsize = 9)+
    theme_classic() )
map2 + guides(color = FALSE )+ 
  scale_color_manual(values=c("#1B9E77","grey25","magenta"))


map2 <- map2 + guides(color = FALSE )



plot_grid(map1+ylim(-.4,1.25), map2, labels = c('A', 'B'), label_size = 9)


ggsave("Supplemental Figure S3.png",plot_grid(map1+ylim(-.4,1.25), map2, labels = c('A', 'B'), label_size = 12), dpi=300, units="cm", width=15, height=10)