rm(list=ls())
library(tidyverse)
library(ggpubr)

# species vectors for plotting:
cols <- c("AI" = "#FDE725FF", "KO" = "#5DC863FF","AC" = "#21908CFF",  "AV" = "#3B528BFF", "SA" = "#440154FF")
speciesnames <- c("AI" = "Acanthus ilicifolius", "KO" = "Kandelia obovata", "AC" = "Aegiceras corniculatum", "AV" = "Avicennia marina", "SA" = "Sonneratia apetala")

# load data
df.size <- data.frame(read.csv("dfbranchdimensions.csv"))
df.size$species <- factor(df.size$species, levels = c("AI", "KO", "AC", "AV","SA"))
df.archi <- data.frame(read.csv("dfarchitecture.csv"))
df.archi$species <- factor(df.archi$species, levels = c("AI", "KO", "AC", "AV","SA"))

df.archi$treeID <- paste(df.archi$site, df.archi$species, df.archi$tree, sep= "")

mean.archi <- df.archi %>%
  group_by(species, branchorder) %>%
  summarise(meandiam = mean(diam_cm, na.rm = TRUE),
            sddiam = sd(diam_cm, na.rm = TRUE))

p.branchorder_diameter <-
  ggplot(df.archi,
         aes(x = as.factor(branchorder), y = diam_cm)) + 
    geom_boxplot(aes(colour = species), lwd= 0.8,  width = 0.8) +
  scale_colour_manual(values = cols, labels = speciesnames) +
  labs(y = "Branch diameter ø (cm)",
       x = "Branch hierarchy (stem to twig)", 
       colour = "Species") +
  theme_classic() +
  scale_y_log10(breaks = seq(from = 0, to = 30, by = 2))

p.branchorder_diameter
  
### --- branch length per diameter -----
# make regression line
mod.len <- lm(length_m ~ diameter_cm, df.size)
diams <- seq(min(df.size$diameter_cm), max(df.size$diameter_cm), 0.1)
pred.len <- predict(mod.len, list(diameter_cm = diams))
df.predlen <- data.frame(diameter_cm = diams, length_m = pred.len)
df.predlen$species <- NA

# make plot
p.length <- 
  ggplot(df.size, 
         aes(x = diameter_cm, 
             y = as.numeric(length_m))) +
  geom_point(aes(colour = species), size = 2) +
  scale_colour_manual(values = cols) +
  scale_fill_manual(values = cols) +
  theme_classic() +
  geom_line(data = df.predlen,  aes(x = diameter_cm, y = length_m)) +
  labs(x = "Diameter ø (cm)", y = expression(paste("Branch length L (m)")), colour = "Species",
       title = "(a) Branch length per base diameter") +
  scale_x_continuous(breaks = seq(0, 5, 1), limits = c(0, 4.5)) 

### --- branch surface area, without leaves -----
# make regression line
mod.area_nl <- lm(formula = area_m2 ~ I(diameter_cm^2), data = filter(df.size, leaves == "no"))
diams <- seq(min(df.size[df.size$leaves == "no",]$diameter_cm), max(df.size[df.size$leaves == "no",]$diameter_cm), 0.1)
pred.area_nl <- predict(mod.area_nl, list(diameter_cm = diams))
df.pred_area_nl <- data.frame(diameter_cm = diams, area_m2 = pred.area_nl)

# make plot
p.area_nl <- 
  ggplot(filter(df.size, leaves == "no"), 
         aes(x = diameter_cm, y = area_m2)) +
  geom_point(aes(colour = species), size = 2) + 
  geom_line(data = df.pred_area_nl, 
            aes(x = diameter_cm, y = area_m2), colour = "black") +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  theme_classic() +
  labs(x = expression(paste("Branch diameter ", symbol("\306"), " (cm)")),
       y = expression(paste("Branch projected surface area ", A[proj], " (", m^2, ")")),
       colour = "Species", title = "(c) Branches without leaves") +
  scale_x_continuous(breaks = seq(0, 4, 1), limits = c(0, 4.5)) +
  scale_y_continuous(breaks = seq(0, 1, 0.1), limits = c(0, 1)) 


### --- branch surface area, with leaves -----
# prepare regression line
mod.area <- lm(area_m2 ~ I(diameter_cm^2), data = filter(df.size, leaves == "yes"))
diams <- seq(min(df.size$diameter_cm), max(df.size$diameter_cm), 0.1)
pred.area <- predict(mod.area, list(diameter_cm = diams))
df.predarea <- data.frame(diameter_cm = diams, area_m2 = pred.area)
df.predarea$species <- NA

# make plot
p.area_wl <- 
  ggplot(filter(df.size, leaves == "yes"), 
         aes(x = diameter_cm, y = area_m2)) +
  geom_point(aes(colour = species), size = 2) +
  theme_light() +
  geom_line(data = df.predarea, 
            aes(x = diameter_cm, y = area_m2)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  scale_fill_manual(values = cols) +
  theme_classic() +
  labs(x = expression(paste("Branch diameter ", symbol("\306"), " (cm)")),
       y = expression(paste("Projected surface area ", A[proj], " (", m^2, ")")),
       colour = "Species", title = "(b) Branches with leaves") +
  scale_x_continuous(breaks = seq(0, 5, 1), limits = c(0, 4.5)) +
  scale_y_continuous(breaks = seq(0, 1, 0.1), limits = c(0, 1)) 

# plot data
ggarrange(p.length, p.area_wl,
          p.area_nl2, 
          ncol = 2, nrow = 2,
          common.legend = TRUE, legend = "right")

### ---- statistics ----
summary(mod.len)
summary(mod.area)
summary(mod.area_nl)