knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE, error = FALSE, eval = TRUE, cache = FALSE, tidy.opts = list(width.cutoff = 60), tidy = TRUE)

#code to set the document active folder (the one in which this .RMD file is located) as the working directory, if it doesn't work just set it manually.
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

comp <- complexity %>% 
  mutate(mini = paste(.$Code, .$Point, sep = "_")) %>% 
  mutate(site = recode(Code,"NBK" = "Rade de Brest - Keraliou", "MBE" = "Belle-ile", "MCM" = "Camaret", "MGL" = "Glenan","MME"= "Meaban","MMO" = "Molene", "MMX" = "Morlaix", "MRZ" = "Rade de Brest - Rozegat", "MPP" = "Baie de Saint-Brieuc", "MTR" = "Trevignon"))%>% 
  mutate(site = factor(site, levels = c("Baie de Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Rade de Brest - Keraliou", "Rade de Brest - Rozegat", "Trevignon"))) %>%
  mutate(Maerl_density = recode(site, "Rade de Brest - Keraliou" =  2288, "Belle-ile" = 1311, "Camaret" = 1322, "Glenan" = 2155, "Meaban" = 2222, "Molene" = 977, "Morlaix" = 3000, "Baie de Saint-Brieuc" = 277, "Rade de Brest - Rozegat" = 3377, "Trevignon" = 988)) %>%  ## add maerl density data at site level
  mutate(Sphericity = ((S^2)/(L*I))^(1/3)) %>% 
  mutate(DR1 = S/L) %>% 
  mutate(DR2 = ((L-I)/(L-S))) %>% 
  mutate(DR3 = I/L) %>% 
  mutate(point_name = paste(.$site, .$Point, sep = " ")) %>% 
  mutate(point_name = gsub("Rade de Brest - ", "", point_name)) %>%
  mutate(point_name = gsub("Baie de ", "", point_name)) %>%
  mutate(point_name = gsub("è", "e", point_name)) %>% 
  mutate(point_name = gsub("é", "e", point_name)) %>% 
  select(mini, site, point_name, everything(), -Code) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.factor)

# Allow us to see intra and intersite variability 

compnum <- bccomp %>% 
  select(Branching_density:DR3)
complexpca <- rda(compnum, scale = TRUE)

#summary(complexpca)
screeplot(complexpca,bstick=TRUE) #only look at pc1 and pc2, maybe pc5

## save as PDF
# pdf("03_Figures/pca_complex_wb.pdf", width = 20, height = 8.75)

pca_victor(pca = complexpca, metadata = bccomp, main.group = "site", scale.fill = pal, scale.colour = pal, goodness.thresh = 0.0, add.centroids = TRUE, nudge.x = c(0, 0.09, -0.055, .14, -0.05, 0, 0.06, -0.14, 0, 0), nudge.y = c(-.05, .06, -.05, .005, 0.055, -.05, -.05, .005, .05, -.05), stat1 = "ellipse", conf.level = .8, ysites = c(-0.65, 0.65), xsites = c(-0.65, 0.65), ysp = c(-2.1, 1.5), xsp = c(-2.1, 1.5), font.size = 20/.pt, axis.size = 22, axis.text = 20, c.size = 4)

# dev.off()


pca_victor(axes = c(1,5),pca = complexpca, metadata = comp, main.group = "site", scale.fill = pal, scale.colour = pal, stat1 = "ellipse" ,ysites = c(-.8, .8), xsites = c(-.8, .8), ysp= c(-2.1,1.6), xsp = c(-2.1,1.6))

#Extract PC1 and PC2 scores for exploration against environmental variables 
bccomp <- bccomp %>% 
  mutate(PC1_score =  scores(complexpca,choices = 1, scaling = 1, display = "sites")) %>% 
  mutate(PC2_score =  scores(complexpca,choices = 2, scaling = 1, display = "sites")) %>% 
  mutate(PC5_score =  scores(complexpca,choices = 5, scaling = 1, display = "sites"))

#get the coordinates
triplot_pos <- comp %>% 
  select(mini,Point, point_name, Sample, site, DR1, DR2, DR3) %>% 
  mutate(x = (1-DR3 + (0.5 * DR1))) %>% # x coordinates based on Graham and Midgley (2000), “Graphical Representation of Particle Shape Using Triangular Diagrams.”
  mutate(y =  ((DR1*.866) * 1.15)) %>% # y coordinates
  mutate(Point = as.factor(Point)) %>% 
  mutate(site = gsub("Rade de Brest - ", "", site)) %>% 
  mutate(site = gsub("Baie de ", "", site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon")))



## Here's an option to create a plot with Sneed & Folk original classes  
# triplot_SF <- ggplot()+ #Sneed and Folk classes
#   geom_polygon(aes(x = c(0, 1, .5, 0), y = c(0, 0, 1, 0)), fill = "grey92")+ #ggplot's defaut bg colour
#   geom_segment(aes(x = .35, y = .7, xend = .65, yend = .7), size = 1, colour = "white") +
#   geom_segment(aes(x = .25, y = .5, xend = .75, yend = .5), size = 1, colour = "white") +
#   geom_segment(aes(x = .15, y = .3, xend = .85, yend = .3), size = 1, colour = "white") +
#   geom_segment(aes(x = .33, y = 0, xend = .45, yend = .7), size = 1, colour = "white") +
#   geom_segment(aes(x = .33, y = 0, xend = .45, yend = .7), size = 1, colour = "white") +
#   geom_segment(aes(x = .67, y = 0, xend = .55, yend = .7), size = 1, colour = "white") +
#   xlim(-0.25, 1.25) +
#   theme(line = element_blank(), panel.background = element_rect(fill = NA))
  
triplot <- ggplot()+ #Here's the option used in figure 3, based on Bosence's Classes
  geom_polygon(aes(x = c(0, 1, .5, 0), y = c(0, 0, 1, 0)), fill = "grey92")+ #ggplot's defaut bg colour
  geom_segment(aes(x = .35, y = .7, xend = .65, yend = .7), size = 1, colour = "white") +
  geom_segment(aes(x = .25, y = .5, xend = .5, yend = .35), size = 1, colour = "white") +
  geom_segment(aes(x = .5, y = .35, xend = .75, yend = .5), size = 1, colour = "white") +
  geom_segment(aes(x = .5, y = 0, xend = .5, yend = .35), size = 1, colour = "white") +
  geom_segment(aes(x = .7, y = 0, xend = .85, yend = .3), size = 1, colour = "white") +
  geom_segment(aes(x = .15, y = 0.3, xend = .3, yend = 0), size = 1, colour = "white") +
  xlim(-0.25, 1.25) +
  theme(line = element_blank(), panel.background = element_rect(fill = NA))+
  coord_fixed(ratio = 1, clip = "off")

# create one plot for each site with a facet grid theme
plist <- list()
pal.tab <- data.frame(site = levels(triplot_pos$site), pal = pal)
triplot_pos <- triplot_pos %>% 
  left_join(pal.tab)
for(i in unique(triplot_pos$site)){
  plist[[i]] <- triplot + geom_point(data  = triplot_pos[triplot_pos$site == i, ], aes(x= x, y = y, col = site, shape = Point), size = 2.5)
  plist[[i]] <- plist[[i]] + facet_wrap(~site)
  plist[[i]] <- plist[[i]] + scale_colour_manual(values = pal.tab[pal.tab$site == i, "pal"])
  plist[[i]] <- plist[[i]] + scale_fill_manual(values = pal.tab[pal.tab$site == i, "pal"])
  plist[[i]] <- plist[[i]] + labs(x = "", y = "") 
  plist[[i]] <- plist[[i]] + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(), legend.position = "NULL", strip.text = element_text(size = 30, face = "bold"))
  plist[[i]] <- plist[[i]] + coord_fixed(ratio = 1, clip = "off")
  plist[[i]]
  }
### Get the means so we have the centroids
centroids <-  triplot_pos %>%
  group_by(site) %>%
  summarise_at(vars(c(x, y)),mean)  %>%
  ungroup() %>% 
  rename(x_cent = x, y_cent = y)
  
triplot_pos <-  triplot_pos %>% 
  mutate(Title_all = as.factor("All Sites")) 

classes <- data.frame(classes = c("SS", "SD", "SE"), x = c(.50, .35, .65), y = c(.55, .24, .24)) %>%
  mutate(Title_sneed = as.factor("Bosence's Classes"))

triplot_pos <- triplot_pos %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon")))
p0 <- triplot
p0 <- p0 + geom_text(data = classes, aes(x = x, y = y, label = classes), size = 18/.pt, fontface = "bold")
p0 <- p0 + labs(x = "L-I/L-S", y = "")
p0 <- p0 + annotate(x = .15, y = .55, geom = "text", label = "S/L", size = 20/.pt, angle = 60)
p0 <- p0 + annotate(x = .85, y = .55, geom = "text", label = "I/L", size = 20/.pt, angle = -60)
p0 <- p0 + annotate(x = -.1, y = -.05, geom = "text", label = "Discoidal", fontface="bold", size = 25/.pt)
p0 <- p0 + annotate(x = 1.1, y = -.05, geom = "text", label = "Ellipsoidal", fontface="bold", size = 25/.pt)
p0 <- p0 + annotate(x = .5, y = 1.05, geom = "text", label = "Spheroidal", fontface="bold", size = 25/.pt)
p0 <- p0 + facet_wrap(~Title_sneed)
p0 <- p0 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(), legend.position = "NULL", axis.title.x = element_text(size = 20, vjust = -.05), strip.text = element_text(size = 30, face = "bold"))
p0 <- p0 + coord_fixed(ratio = 1, clip = "off")

p11 <- triplot + stat_ellipse(data  = triplot_pos, aes(x= x, y = y, colour = site, fill = site), alpha = 0.3, geom = "polygon", level = .75, linetype = 2)
  p11 <- p11 + geom_point(data = centroids, aes(x = x_cent, y = y_cent, fill = site), size = 3.5, shape = 23, alpha = .9, colour = "NA")
  p11 <- p11 + facet_wrap(~Title_all)
  p11 <- p11 + scale_colour_manual(values = pal)
  p11 <- p11 + scale_fill_manual(values = pal)
  #trick so that p11 and p0 have the same size
  p11 <- p11 + labs(x = " ", y = " ") 
  p11 <- p11 + labs(x = "L-I/L-S", y = "")
  p11 <- p11 + annotate(x = .15, y = .55, geom = "text", label = "S/L", size = 20/.pt, angle = 60, col = "NA")
  p11 <- p11 + annotate(x = .85, y = .55, geom = "text", label = "I/L", size = 20/.pt, angle = -60, col = "NA")
  p11 <- p11 + annotate(x = -.1, y = -.05, geom = "text", label = "Discoidal", fontface="bold", size = 25/.pt, col = "NA")
  p11 <- p11 + annotate(x = 1.1, y = -.05, geom = "text", label = "Ellipsoidal", fontface="bold", size = 25/.pt, col = "NA")
  p11 <- p11 + annotate(x = .5, y = 1.05, geom = "text", label = "Spheroidal", fontface="bold", size = 25/.pt, col =  "NA")
  #again here in theme we add colour = "NA"so that the text occupies the same space as in p0 but is transparent...
  p11 <- p11 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.title.x = element_text(size = 20, vjust = -.05, colour = "NA"), axis.ticks = element_blank(), axis.text = element_blank(), legend.position = "NULL", strip.text = element_text(size = 30, face = "bold"))
  p11 <- p11 + coord_fixed(ratio = 1, clip = "off")


pind <- ggarrange(plist$`Saint-Brieuc`, plist$`Belle-ile`,plist$Camaret, plist$Glenan, plist$Meaban, plist$Molene, plist$Morlaix, plist$Keraliou, ncol = 2, nrow = 4)
ptry <- ggarrange(ggarrange(p11, p0, ncol = 1, nrow = 2), ggarrange(plist$Rozegat, plist$Trevignon, align = "hv"), ncol = 1, heights = c(3,1))
(pf <- ggarrange(pind, ptry , widths = c(1,1), heights = c(1,1), align = "hv"))  

# pdf("03_Figures/shapetriplots.pdf", height = 16, width = 20)
# pf
# dev.off()

#---- Permanova with all variables between sites ----

compnum_sc <- as.data.frame(scale(compnum)) ## scale the transformed variables
comp_dis <- vegdist(compnum_sc, method = "eucl") ## create eucliden distance matrix
comp_bd <- betadisper(comp_dis,comp$site) ## test for homogeneity of multivariate dispersions

permutest(comp_bd)# Some sites are different from eachother (p < 0.001)

tukey <- TukeyHSD(comp_bd) 
tukey <- do.call(rbind.data.frame, tukey) 

tukey %>% 
  filter(`p adj` < 0.05) #Problematic sites are Molene Saint-Brieuc

#
permcomp <- adonis2(formula= compnum_sc ~ comp$site/comp$point_name, method = "eucl") #fixed effects with points nested in sites
permcomp

pairwiseAdonis::pairwise.adonis(comp_dis, comp$site, p.adjust.m = "holm") #all sites are significantly different.

## Sites Saint-Brieuc and Molene are the ones that pose most of the problems

sansSB <- compnum_sc %>% 
  filter(!bccomp$site %in% c("Saint-Brieuc", "Molene"))
compsansSB <- bccomp %>% 
  filter(!site %in% c("Saint-Brieuc", "Molene"))
compsansSB <- droplevels(compsansSB)
sansSB_dis <- vegdist(sansSB, method = "eucl") ## create eucliden distance matrix
sansSB_bd <- betadisper(sansSB_dis, compsansSB$site) ## test for homogeneity of multivariate dispersions

permutest(sansSB_bd)# MET!

#(tukeySB <- TukeyHSD(sansSB_bd)) 


permcompSB <- adonis2(formula= sansSB_dis ~ compsansSB$site/compsansSB$point_name, method = "eucl") #fixed effects with point nested in sites
permcompSB

## almost nothing changes from the first result! we can proceed with the interpretation assuming the differences in dispersion are not interefering as much.

granulo <-  point_imp %>% 
  filter(site!= "Rhuys") %>%     #remove rows on Rhuys site, not analysed in this project, only spring data
  filter(annee %in% c(2007:2018)) %>% 
  mutate(point_name = gsub("Belle-Ile", "Belle-ile", point_name)) %>% 
  mutate(point_name = gsub("Rade de Brest - ", "", point_name)) %>%
  mutate(point_name = gsub("Rade de Brest-", "", point_name)) %>%
  mutate(point_name = gsub("è", "e", point_name)) %>% 
  mutate(point_name = gsub("é", "e", point_name)) %>% 
   mutate(site = gsub("Rade de Brest - ", "", .$site)) %>% 
  mutate(site = gsub("Baie de ", "", .$site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  select(mini, everything()) %>%
  as.data.frame() %>%
  set_rownames(paste( .$mini, .$annee)) %>% 
  rename(OM = pourcentage_MO) %>% 
  mutate_if(is.character, as.factor)

granulo <- droplevels(granulo)

label_pos <- granulo %>% 
  select(Gravel:mud) %>% 
  aggregate(., by = list(granulo$site), FUN = mean) %>% 
  rename(., site = Group.1) %>% 
  mutate(site = gsub("Rade de Brest - ", "", site)) %>% 
  mutate(site = gsub("Baie de ", "", site)) %>% 
  mutate(total = Sand + Gravel + mud) %>% 
  mutate_if(is.character, as.factor) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon")))
  
#all samples
p <- ggtern(data = granulo, mapping = aes(x = Sand, y= mud, z = Gravel))+
  geom_point(aes(colour = site))+
  stat_mean_ellipse(aes(colour = site, fill = site), geom = "polygon", alpha = 0.3)+
  geom_label_viewport(x = .75, y = .27, label = "Trevignon", fill = pal[10], fontface="bold", size = 20/.pt)+
  geom_label_viewport(x = .65, y = .22, label = "Rozegat", fill = pal[9], fontface="bold", size = 20/.pt)+
  geom_label_viewport(x = .62, y = .42, label = "Keraliou", fill = pal[8], fontface="bold", size = 20/.pt)+
  geom_label_viewport(x = .34, y = .19, label = "Camaret", fill = pal[3], fontface="bold", size = 20/.pt)+
  scale_fill_manual(values = pal) +
  scale_colour_manual(values = pal) +
  labs(x = "Sand", y = "Mud", z = "Gravel", )+
  theme(legend.position = "null", axis.title = element_text(face = "bold", size = 22), axis.text = element_text(size = 20))+
  theme_arrowlength()
  

#pdf(file = "03_Figures/triplot_granulo.pdf", width = 20)
p

#dev.off()

  
# Just an idea of what it would look like only with mean values

ggtern(data = granulo_mean, mapping = aes(x = Sand, y= mud, z = Gravel))+
  geom_point(aes(colour = site))+
  stat_mean_ellipse(aes(colour = site, fill = site), geom = "polygon", alpha = 0.3)+
  scale_fill_manual(values = pal) +
  scale_colour_manual(values = pal) +
  labs(x = "Sand", y = "Mud", z = "Gravel", )+
  theme_arrowlength()

hydro <- MA_hydro_int_six_month %>% 
  select(-c(theme, method_echant))%>% #remove useless info
  filter(site!= "Rhuys",saison == "Printemps") %>%   #remove rows on Rhuys site, not analysed in this project, only spring data
  filter(annee %in% c(2007:2018)) %>% 
  mutate(mini = recode(point_name, "Keraliou 1"= "NBK_1","Keraliou 2"= "NBK_2","Keraliou 3"= "NBK_3","Belle-Ile 1"= "MBE_1", "Belle-Ile 2"= "MBE_2", "Belle-Ile 3"= "MBE_3", "Camaret 1"= "MCM_1", "Camaret 2"= "MCM_2", "Camaret 3"= "MCM_3", "Glénan 1"= "MGL_1", "Glénan 2"= "MGL_2", "Glénan 3"= "MGL_3", "Meaban 1"= "MME_1", "Meaban 2"= "MME_2", "Meaban 3"= "MME_3", "Molène 1"= "MMO_1", "Molène 2"= "MMO_2", "Molène 3"= "MMO_3", "Morlaix 1"= "MMX_1", "Morlaix 2"= "MMX_2", "Morlaix 3"= "MMX_3", "Rade de Brest-Rozegat 1" = "MRZ_1", "Rade de Brest-Rozegat 2" = "MRZ_2", "Rade de Brest-Rozegat 3" = "MRZ_3", "Saint-Brieuc 1" = "MPP_1", "Saint-Brieuc 2" = "MPP_2", "Saint-Brieuc 3" = "MPP_3", "Trévignon 1" = "MTR_1", "Trévignon 2" = "MTR_2", "Trévignon 3" = "MTR_3" )) %>%
  mutate(point_name = gsub("Belle-Ile", "Belle-ile", point_name)) %>% 
  mutate(point_name = gsub("Rade de Brest - ", "", point_name)) %>%
  mutate(point_name = gsub("Rade de Brest-", "", point_name)) %>%
  mutate(point_name = gsub("è", "e", point_name)) %>% 
  mutate(point_name = gsub("é", "e", point_name)) %>% #remove special ch.
  mutate(site = gsub("Rade de Brest - ", "", .$site)) %>% 
  mutate(site = gsub("Baie de ", "", .$site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  mutate(annee = as.numeric(levels(.$annee)[.$annee])) %>% 
  select(mini, everything()) %>%
  as.data.frame() %>%
  set_rownames(paste( .$mini, .$annee)) %>% 
  mutate_if(is.character, as.factor)

hydro <- droplevels(hydro)

# ggplot(data = hydro, aes(x = point_name, y = bottomT_mean, colour = site))+
#   geom_boxplot()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))
# 
# ggplot(data = hydro, aes(x = point_name, y = current_mean, colour = site))+
#   geom_boxplot()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))


#let's get the mean for now

hydro_mean <- hydro %>% 
  select(bottomT_min:current_sd) %>% 
  aggregate(., by = list(hydro$mini, hydro$site, hydro$point_name), FUN = median) %>% 
  rename(., mini = Group.1, site = Group.2, point_name = Group.3) %>% 
  set_rownames(paste( .$mini))

# ggplot(data = hydro_mean, aes(x = point_name, y = current_mean, colour = site))+
#   geom_point()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))
# ggplot(data = hydro_mean, aes(x = point_name, y = bottomT_mean, colour = site))+
#   geom_point()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))
# ggplot(data =hydro_mean, aes(x = current_max))+
#   stat_density(col = "red")

#let's get the sd 
hydro_sd <- hydro %>% 
  select(bottomT_min:current_sd) %>% 
  aggregate(., by = list(hydro$mini, hydro$site, hydro$point_name), FUN = sd) 
colnames(hydro_sd) <- paste(colnames(hydro_sd), "sd", sep = "_")
hydro_sd <- hydro_sd %>% 
  rename(., mini = Group.1_sd, site = Group.2_sd, point_name = Group.3_sd) %>% 
  mutate_if(is.character, as.numeric)

# ggplot(data = hydro_sd, aes(x = point_name, y = current_mean_sd, colour = site))+
#   geom_point()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))
# ggplot(data = hydro_sd, aes(x = point_name, y = bottomT_mean_sd, colour = site))+
#   geom_point()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))

bathy <- depth_point_sub %>% 
  filter(theme == "Bancs de Maërl") %>% 
  select(-c(theme, method_echant))%>% #remove useless info
  filter(site!= "Rhuys") %>%   #remove rows on Rhuys site, not analysed in this project
  mutate(mini = recode(point_name, "Keraliou 1"= "NBK_1","Keraliou 2"= "NBK_2","Keraliou 3"= "NBK_3","Belle-Ile 1"= "MBE_1", "Belle-Ile 2"= "MBE_2", "Belle-Ile 3"= "MBE_3", "Camaret 1"= "MCM_1", "Camaret 2"= "MCM_2", "Camaret 3"= "MCM_3", "Glénan 1"= "MGL_1", "Glénan 2"= "MGL_2", "Glénan 3"= "MGL_3", "Meaban 1"= "MME_1", "Meaban 2"= "MME_2", "Meaban 3"= "MME_3", "Molène 1"= "MMO_1", "Molène 2"= "MMO_2", "Molène 3"= "MMO_3", "Morlaix 1"= "MMX_1", "Morlaix 2"= "MMX_2", "Morlaix 3"= "MMX_3", "Rade de Brest-Rozegat 1" = "MRZ_1", "Rade de Brest-Rozegat 2" = "MRZ_2", "Rade de Brest-Rozegat 3" = "MRZ_3", "Saint-Brieuc 1" = "MPP_1", "Saint-Brieuc 2" = "MPP_2", "Saint-Brieuc 3" = "MPP_3", "Trévignon 1" = "MTR_1", "Trévignon 2" = "MTR_2", "Trévignon 3" = "MTR_3" )) %>%
  mutate(point_name = gsub("Belle-Ile", "Belle-ile", point_name)) %>% 
  mutate(point_name = gsub("Rade de Brest - ", "", point_name)) %>%
  mutate(point_name = gsub("Rade de Brest-", "", point_name)) %>%
  mutate(point_name = gsub("è", "e", point_name)) %>% 
  mutate(point_name = gsub("é", "e", point_name)) %>% #remove special ch.
  mutate(site = gsub("Rade de Brest - ", "", .$site)) %>% 
  mutate(site = gsub("Baie de ", "", .$site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  mutate(depth = -depth) %>% 
  select(mini, everything()) %>%
  as.data.frame() %>%
  set_rownames(paste( .$mini)) %>% 
  mutate_if(is.character, as.factor)

bathy <- droplevels(bathy)

# ggplot(data = bathy, aes(x = point_name, y = depth, colour = site))+
#   geom_point()+
#   scale_color_manual(values = pal)+
#   theme(axis.text.x = element_text(angle = 90))

fetch <- fetch.dfDCE %>% 
  mutate(site = recode(site, "Trevignon " = "Trevignon")) %>%
  mutate(point_name = paste(site, point, sep = " ")) %>% 
  mutate(mini = recode(point_name, "Rade de Brest - Keraliou 1"= "NBK_1","Rade de Brest - Keraliou 2"= "NBK_2","Rade de Brest - Keraliou 3"= "NBK_3","Belle-ile 1"= "MBE_1", "Belle-ile 2"= "MBE_2", "Belle-ile 3"= "MBE_3", "Camaret 1"= "MCM_1", "Camaret 2"= "MCM_2", "Camaret 3"= "MCM_3", "Glenan 1"= "MGL_1", "Glenan 2"= "MGL_2", "Glenan 3"= "MGL_3", "Meaban 1"= "MME_1", "Meaban 2"= "MME_2", "Meaban 3"= "MME_3", "Molene 1"= "MMO_1", "Molene 2"= "MMO_2", "Molene 3"= "MMO_3", "Morlaix 1"= "MMX_1", "Morlaix 2"= "MMX_2", "Morlaix 3"= "MMX_3", "Rade de Brest - Rozegat 1" = "MRZ_1", "Rade de Brest - Rozegat 2" = "MRZ_2", "Rade de Brest - Rozegat 3" = "MRZ_3", "Baie de Saint-Brieuc 1" = "MPP_1", "Baie de Saint-Brieuc 2" = "MPP_2", "Baie de Saint-Brieuc 3" = "MPP_3", "Trevignon 1" = "MTR_1", "Trevignon 2" = "MTR_2", "Trevignon 3" = "MTR_3" )) %>%
  mutate(point_name = gsub("Rade de Brest - ", "", point_name)) %>%
  mutate(point_name = gsub("Rade de Brest-", "", point_name)) %>%
  mutate(point_name = gsub("Baie de ", "", point_name)) %>%
    mutate(site = gsub("Rade de Brest - ", "", site)) %>% 
  mutate(site = gsub("Baie de ", "", site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  mutate(site = gsub("Rade de Brest - ", "", .$site)) %>% 
  mutate(site = gsub("Baie de ", "", .$site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  select(mini, site, point_name, everything()) %>% 
  set_rownames(paste( .$mini)) %>% 
  mutate_if(is.character, as.factor)

fetch.melt <- melt(fetch.dfDCE,id.vars=c("site","point","longitude","latitude"))

# Retrieve a map layer for France
france_data <- map_data("france")
names(france_data) <- c("X","Y","PID","POS","region","subregion")

# Zoom on Brittany
xlim <- c(-7.5,-1.8)
ylim <- c(46,49)
francemap <- clipPolys(france_data, xlim=xlim,ylim=ylim, keepExtra=TRUE)


p <- ggplot(aes(x=longitude,y=latitude,size=value, col = site), data=fetch.melt, fill = "black") + coord_equal()
p <- p + geom_polygon(data=francemap,aes(X,Y,group=PID),fill="gray35",inherit.aes=F,show.legend=F)
p <- p + geom_point(aes(text=paste(site,point,value)), alpha=0.3)
p <- p + facet_wrap(~variable,ncol=3) 
p <- p + scale_fill_manual(values = pal)
p <- p + scale_colour_manual(values = pal)
p <- p + theme(legend.position="none",panel.background = element_rect(fill = "#b6ced6"), panel.grid = element_blank())
p <- p +  ylim(c(min(francemap$Y), max(francemap$Y)))
p <- p + ylim(c(min(francemap$Y), max(francemap$Y)))
p <- p + ylab("")+ xlab("")
p <- p +scale_y_continuous(breaks=c(46.5,47,47.5,48,48.5,49,49.5,50), labels = c("46.5°N","47ºN","47.5°N", "48ºN","48.5°N", "49ºN","49.5°N","50°N")) 
p <- p + scale_x_continuous(breaks=c(-6,-5,-4,-3,-2,-1,0), labels = c("6°W","5°W","4°W","3°W","2°W","1°W","0°")) 
p <- p + coord_cartesian(xlim = xlim, ylim = ylim, expand = FALSE)

p

# get the fetch max 
fetch_max <- fetch %>% 
  gather(North:West, key = "Direction", value = Fetch) %>% 
  group_by(mini, point_name, latitude, longitude) %>% 
  summarise(Fetch_max = max(Fetch)) %>% 
  left_join(fetch) %>% 
  select(mini, site, point_name, -point, latitude, longitude, North:Average, Fetch_max) %>% 
  set_rownames(paste(.$mini))

p <- ggplot() 
p <- p + geom_polygon(data=francemap,aes(X,Y,group=PID),fill="gray90",inherit.aes=F,show.legend=F)
p <- p + geom_point(data = fetch_max, aes(x=longitude,y=latitude,size=Fetch_max, fill=site), shape = 21, alpha  = 0.6) + coord_equal() + scale_size_continuous(range = c(2,15), name = "Fetch", breaks = c(50, 150, 300))
p <- p + scale_fill_manual(values = pal)
p <- p + guides(fill = FALSE, size = guide_legend(title.position = "top"))
p <- p + theme(legend.key = element_rect(fill = NA), legend.background = element_blank(), panel.background = element_rect(fill = "#b6ced6"), panel.grid = element_blank(), axis.text = element_text(size = 14))
p <- p + ylab("")+ xlab("")
p <- p +scale_y_continuous(breaks=c(46.5,47,47.5,48,48.5,49,49.5,50), labels = c("46.5°N","47ºN","47.5°N", "48ºN","48.5°N", "49ºN","49.5°N","50°N")) 
p <- p + scale_x_continuous(breaks=c(-7, -6,-5,-4,-3,-2,-1,0), labels = c( "7°W","6°W","5°W","4°W","3°W","2°W","1°W","0°")) 
p <- p + coord_cartesian(xlim = xlim, ylim = ylim, expand = FALSE)

p

# Color Brittany (Bretagne)
fill_bret <- rep("gray90",nrow(france_data))
fill_bret[france_data$region=="Finistere"] <- "gray35"
fill_bret[france_data$region=="Morbihan"] <- "gray35"
fill_bret[france_data$region=="Ille-et-Vilaine"] <- "gray35"
fill_bret[france_data$region=="Cotes-Darmor"] <- "gray35"


france_map <- ggplot() + coord_map()
france_map <- france_map + scale_x_continuous(expand=c(0,0))
france_map <- france_map + theme(axis.ticks=element_blank(), axis.title=element_blank(), axis.text=element_blank())
france_map <- france_map + geom_polygon(data=france_data, mapping=aes(x=X, y=Y,group=PID), fill=fill_bret,col=fill_bret) 
france_map <- france_map + theme_void() 

### world

world <- map_data("world")

# Color France
fill_world <- rep("gray90",nrow(world))
fill_world[world$region=="France"] <- "gray35"
col_world <- rep("gray75", nrow(world) )
col_world[world$region=="France"] <- "gray35"
worldmap <- ggplot(world, aes(x=long, y=lat, group=group)) 
worldmap <- worldmap + geom_path(size = .2)
worldmap <- worldmap + scale_y_continuous(breaks=(-2:2) * 30) 
worldmap <- worldmap + scale_x_continuous(breaks=(-4:4) * 45)+ theme(panel.grid = element_line(colour = "gray35", size = 0.2))
worldmap <- worldmap + coord_map("ortho", orientation=c(41, -5, 0)) + ylab("") + xlab("")
worldmap <- worldmap + geom_polygon(data=world, mapping=aes(x=long, y=lat,group=group), fill=fill_world,col=col_world,size = 0.2) 
worldmap <- worldmap+ theme(axis.text=element_blank(),axis.ticks=element_blank(),panel.border=element_blank(),panel.background = element_blank(), rect = element_blank())


fetch_max <- fetch_max %>%
  as.data.frame() %>% 
  mutate(rename(., long = longitude, lat = latitude))


(map_fetch <- p +  
    geom_rect(aes(xmin = -6.3, xmax = -4.56, ymin = 46.1, ymax = 47.5), fill = alpha("gray35", .15), color = "gray35", lwd = 0.75)+
    annotation_custom(grob = ggplotGrob(worldmap), xmin = -7.4, xmax = -6, ymin = 47.5, ymax = 48.7)+
    geom_rect(aes(xmin = -6.7, xmax = -6.52, ymin = 48.1, ymax = 48.25), fill = NA, color = "gray35", lwd = 0.75)+
    geom_line(aes(x = c(-6.3, -6.7), y = c(47.5, 48.25)), color = "gray35", lwd = 0.65, linetype = 2) +
    geom_line(aes(x = c(-6.52, -4.56), y = c(48.25, 47.5)), color = "gray35", lwd = 0.65, linetype = 2) +
    annotation_custom(grob = ggplotGrob(france_map), xmin = -6.3, xmax = -4.6, ymin = 46.1, ymax = 47.5 )+
    scalebar(data=fetch_max, dist=45, dist_unit="km",  height = 0.02, transform=TRUE, model="WGS84", box.fill = c("gray35", "white"), location = "bottomright", st.color = "gray35", st.size = 8,  anchor = c(x = -2.38, y =  46.1))+
  annotation_north_arrow(location = "br", which_north = "true", height = unit(1.5, "cm"), width = unit(1.5, "cm"), style = north_arrow_fancy_orienteering(fill = c("white", "gray35"), text_face = "bold", text_col = "gray35", text_size = 16, line_col = "gray35", line_width = 2))+
  geom_text(aes(label = "France", x = -5.42, y = 46.85), fontface = "bold", colour = "gray35", size = 12)+
  geom_text(aes(label = "Brittany", x = -2.9, y = 48.20), fontface = "bold", colour = "gray35", size = 16) +
  geom_text(aes(label = "Saint-Brieuc", x = -2.35, y = 48.65), fontface = "bold", colour = "gray35", size = 10)+
  geom_text(aes(label = "Rozegat", x = -4, y = 48.315), fontface = "bold", colour = "gray35", size = 10) +
  geom_text(aes(label = "Keraliou", x = -4.06, y = 48.415), fontface = "bold", colour = "gray35", size = 10)+
  geom_text(aes(label = "Morlaix", x = -3.5, y = 48.7), fontface = "bold", colour = "gray35", size = 10) +
  geom_text(aes(label = "Molène", x = -5.3, y = 48.4), fontface = "bold", colour = "gray35", size = 10) +
  geom_text(aes(label = "Camaret", x = -5.05, y = 48.25), fontface = "bold", colour = "gray35", size = 10) +
  geom_text(aes(label = "Trévignon", x = -3.35, y = 47.82), fontface = "bold", colour = "gray35", size = 10)+
  geom_text(aes(label = "Glénan", x = -4.4, y = 47.7), fontface = "bold", colour = "gray35", size = 10)+ 
  geom_text(aes(label = "Méaban", x = -2.5, y = 47.53), fontface = "bold", colour = "gray35", size = 10)+ 
  geom_text(aes(label = "Belle-île", x = -3.2, y = 47.22), fontface = "bold", colour = "gray35", size = 10)+
  theme(legend.position = c(.81, .1), legend.direction = "horizontal", legend.title = element_text(face = "bold", colour = "gray35", size = 22), legend.text = element_text(size = 20, colour = "gray35"), axis.text = element_text(size = 20))
  )

# pdf("03_Figures/fetch_map.pdf", width= 16, height=11.25)
# map_fetch
# dev.off()

env <- granulo %>% 
  left_join(hydro) %>% 
  left_join(bathy) %>% 
  left_join(fetch_max) %>% 
  mutate(site = gsub("Rade de Brest - ", "", site)) %>% 
  mutate(site = gsub("Baie de ", "", site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  mutate(site = gsub("Rade de Brest - ", "", .$site)) %>% 
  mutate(site = gsub("Baie de ", "", .$site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  rename("Grain size" = "Mean.fw.um", "Temperature" = "bottomT_mean", "Organic Matter" = "OM", "Current" = "current_mean", "Depth" = "depth", "Fetch" = "Fetch_max", "Trask" = "Trask(So)", "Mud" = "mud")

  
envnum <- env %>% 
  select(Temperature, Current, Trask, `Organic Matter`, Sand, Gravel, Mud, Fetch, Depth, `Grain size`)

boxcoxenv <-   envnum %>% 
  select(-`Grain size`) %>% 
  mutate(Mud = boxcoxTransform(Mud, boxcox(.$Mud, optimize = TRUE )$lambda)) %>%
  mutate(Sand = boxcoxTransform(Sand, boxcox(.$Sand, optimize = TRUE )$lambda)) %>%
  mutate(Gravel = boxcoxTransform(Gravel, boxcox(.$Gravel, optimize = TRUE )$lambda)) %>%
  mutate(Current = boxcoxTransform(Current, boxcox(.$Current, optimize = TRUE )$lambda)) %>%
  mutate(Temperature = boxcoxTransform(Temperature, boxcox(.$Temperature, optimize = TRUE )$lambda)) %>%
  mutate(`Organic Matter` = boxcoxTransform(`Organic Matter`, boxcox(.$`Organic Matter`, optimize = TRUE )$lambda)) %>%
  mutate(Trask = boxcoxTransform(Trask, boxcox(.$Trask, optimize = TRUE )$lambda)) %>%
  mutate(Depth = boxcoxTransform(Depth, boxcox(.$Depth, optimize = TRUE )$lambda)) %>%
  mutate(Fetch = boxcoxTransform(Fetch, boxcox(.$Fetch, optimize = TRUE )$lambda))

bcenvpca <- rda(boxcoxenv, scale = TRUE)

#summary(bcenvpca)
screeplot(bcenvpca,bstick=TRUE) #only look at pc1 and pc2, maybe pc5

# pdf("03_Figures/pca_env_wb.pdf", width = 20, height = 8.75)

pca_victor(pca = bcenvpca, metadata = env, main.group = "site", scale.fill = pal, scale.colour = pal, scale.shape = c(15:17), nudge.x = c(-.11, .17, 0, -0.15, .08, -.08, 0.04, 0, 0, 0 ), nudge.y = c(-.05, .01, .065, .01, .065, .065, -.06, .065, -.05, -.05), stat1 = "ellipse", conf.level = .9, ysites = c(-0.65, 0.65), xsites = c(-0.65, 0.65), ysp = c(-2.3, 2.3), xsp = c(-2.3, 2.3), font.size = 20/.pt, axis.size = 22, axis.text = 20, c.size = 4, add.centroids = TRUE, goodness.thresh = 0)

# dev.off()

#let's check the axis 5 just in case
pca_victor(pca = bcenvpca, metadata = env, main.group = "site", scale.fill = pal, scale.colour = pal, scale.shape = c(15:17), stat1 = "ellipse", conf.level = .9, ysites = c(-0.65, 0.65), xsites = c(-0.65, 0.65), ysp = c(-1.8, 1.8), xsp = c(-1.8, 1.8), font.size = 20/.pt, axis.size = 22, axis.text = 20, c.size = 4, add.centroids = TRUE,axes = c(1,5))

modtab <- bccomp %>% 
  left_join(granulo_mean) %>% 
  left_join(granulo_sd) %>% 
  left_join(hydro_mean) %>% 
  left_join(hydro_sd) %>% 
  left_join(bathy) %>% 
  left_join(fetch_max) %>% 
  mutate(site = gsub("Rade de Brest - ", "", site)) %>% 
  mutate(site = gsub("Baie de ", "", site)) %>% 
  mutate(site = factor(site, levels = c("Saint-Brieuc", "Belle-ile", "Camaret", "Glenan", "Meaban", "Molene", "Morlaix", "Keraliou", "Rozegat", "Trevignon"))) %>% 
  set_rownames(paste(.$mini, .$Sample, sep = "_"))

GGally::ggpairs(modtab %>% select(bottomT_mean, bottomT_max, bottomT_sd)) #Max temperature and variation in temperature are very correlated, but not mean temperature

GGally::ggpairs(modtab %>% select(current_mean, current_sd, current_max)) # all very correlated

varparbc <- varpart(Y = comprda, X = ~ Mud + Trask + Sand, ~ Fetch + Current, ~ Depth,  data = boxcoxenvrda, permutations = 9999)

varparbc

# pdf(file = "03_Figures/varpar.pdf", width = 20, height = 12)
plot(varparbc, digits = 2, bg = c("#02998B","#380142", "#e8bb5a"), lwd = 1.5, cex = 2, cutoff = 0.001)

# dev.off()

## Check it with sites
(bccompenvirdaS <- rda(comprda ~ Mud + Fetch + Sand + Depth + Current + Trask + site, data = boxcoxenvrda, scale = TRUE))

anova(bccompenvirdaS,permutations = how(nperm = 9999))

(R2a.allS <- RsquareAdj(bccompenvirdaS)$adj.r.squared) #added 10% of explained variation in comparison to the model without sites, but could be biased due to influence of currents that are basically the same in points of each site and maerl densities which are at the site level.

# pdf(file = "03_Figures/rda.pdf", width= 16, height=11.25)

autoplot.rda.victor(bccompenvirdaS, metadata = modtab, scale.fill = pal, scale.colour = pal, thresh = 0, stat = "ellipse", lvl = .9, legend.position = c(.4,.9), ylim = c(-2,2.5), xlim = c(-2,2.5))+
  theme(legend.title = element_text(face = "bold", colour = "gray35", size = 14), legend.text = element_text(size = 14, colour = "gray35"), legend.direction = "horizontal", legend.background = element_blank(), legend.box.background = element_blank(), legend.key = element_blank(), ) + guides(fill = guide_legend(title.position = "top"))

# dev.off()

(varparbcS <- varpart(Y = comprda, X = ~ Mud + Trask + Sand, ~ Fetch + Current, ~ Depth, ~site,  data = boxcoxenvrda, permutations = 9999))

#pdf(file = "03_Figures/varpar.pdf", width = 20)
plot(varparbcS, digits = 2, bg = c("#02998B","#380142", "#e8bb5a", "black"), cutoff = 0.001)

#dev.off()