# This script allows the creation of trend figures for the 4 birds species retained in the paper
# These figures are slightly modified compared to the ones obtained in the routine
# So that they are more readable (e.g, increased writing size)
# NB : the script works if the working directory is the root of the indicator routine

rm(list = ls()) 

# Result path
result_path = "C:/Users/mvimont01/Downloads/maj_appendix_birds/"

# Load libraries
library(tidyverse)
library(ggalluvial)

# FUNCTIONS ----

# Load all functions required from the routine
source("R/formatTrendEstimates.R")
source("R/formatYearlyVariations.R")
source("R/identifyConvIssue.R")
source("R/makeSlopesFromVar.R")
source("R/supportFunctions.R")
source("R/fillAbsence.R")

# Updated version of the function 'percEvolutionToText' 
percEvolutionToText_bis <- function(data, distribution, language){
  
  # Extract percentage of change (and confidence interval)
  dataChange = apply(data[,c("estimate", "infIC", "supIC")], 2, function(x){
    
    # Calculate value for first year
    start = x * 0
    
    # Calculate Value for relative last year
    stop = x * (data$maxYear - data$minYear)
    
    # Turn to exponential if counting or occurrence data
    if(distribution != "gaussian"){
      start = exp(start)
      stop = exp(stop)
    }
    
    # Measure the percentage of difference
    perc = round(100 * (stop- start), 1)
    
    # Return the percentage of difference
    return(perc)
  })
  
  # Turn to text
  changeText = paste0(ifelse(dataChange["estimate"] > 0, "+", ""),
                      dataChange["estimate"], "% (",
                      
                      ifelse(dataChange["infIC"] > 0, "+", ""),
                      dataChange["infIC"], "% ; ",
                      
                      ifelse(dataChange["supIC"] > 0, "+", ""),
                      dataChange["supIC"], "%)")
  
  # Add a 
  if(language == "french"){ 
    preText = "Tendance linéaire "}
  else{preText = "Linear trend over "}
  
  changeText = paste0(preText, data$min, " - ", data$max, " : ", changeText)
  
  return(changeText)
}

# Updated version of the function 'plotLinearTrends' 
plotLinearTrends_bis <- function(speciesList, dataLongTerm, dataYearlyVariations, dataShortTerm = NULL,
                                 language = "french", dataName = NULL, distribution, plotGamm = TRUE, plotST = TRUE,
                                 uncertainty = TRUE, N = 100, weight = TRUE, save = TRUE, path = ""){
  
  plots <- sapply(speciesList, function(sp){ 
    # Filter tables for considered species
    dataLT_sp = dataLongTerm[dataLongTerm$species == sp,]
    dataVar_sp = dataYearlyVariations[dataYearlyVariations$species == sp,]
    
    if(nrow(dataLT_sp)==0|nrow(dataLT_sp)==0) {
      message('No plot for species ', sp, ': no model was found for either ',
              'long-term trend or yearly variations\n')
      plot = NULL
      
    }else if(any(is.na(c(dataLT_sp$estimate, dataVar_sp$estimate)))){
      message('No plot for species ', sp, ' due to no-convergence of either',
              " long-term trend or yearly variations\n")
      plot = NULL
      
    }else{
      # Extract values of the temporal series
      yearValues = dataLT_sp$minYear:dataLT_sp$maxYear
      
      ###################
      # GAMM VARIATIONS #
      ###################
      pathToGamm = here::here("outputs", repo, "models", "gammVariations", paste0(sp, ".rdata"))
      
      if(plotGamm & file.exists(pathToGamm)){
        # Load the corresponding file
        load(here::here("outputs", repo, "models", "gammVariations", paste0(sp, ".rdata")))
        
        if(is.null(gammVariations$warnings) & is.null(gammVariations$error)){
          # Extract values of gam
          valuesToPlotGAM = plot(gammVariations$value$gam, se = TRUE, n = 50 * length(yearValues))
          
          # Make data.frame
          dataGammVariations_sp = data.frame(year = valuesToPlotGAM[[1]]$x, 
                                             estimate = valuesToPlotGAM[[1]]$fit, 
                                             infIC = valuesToPlotGAM[[1]]$fit - 1.96 * valuesToPlotGAM[[1]]$se,
                                             supIC = valuesToPlotGAM[[1]]$fit + 1.96 * valuesToPlotGAM[[1]]$se,
                                             group = "GAM",
                                             rep = 1)
          
          # Centrer sur 0
          dataGammVariations_sp[,c("estimate", "infIC", "supIC")] = dataGammVariations_sp[,c("estimate", "infIC", "supIC")] - mean(dataGammVariations_sp$estimate)
          
          # Passer à l'exponentielle si données d'occurrence ou de comptage
          if(distribution != "gaussien"){
            dataGammVariations_sp[,c("estimate", "infIC", "supIC")] = exp(dataGammVariations_sp[,c("estimate", "infIC", "supIC")])
          }
        }else{
          plotGamm = FALSE
        }
        
      }else{
        plotGamm = FALSE
      }
      
      ####################
      # SHORT-TERM TREND #
      ####################    
      if(plotST & !is.null(dataShortTerm)){
        dataST_sp = dataShortTerm[dataShortTerm$species == sp,]
        if(!is.na(dataST_sp$estimate)){
          plotST = TRUE
          
          # Extract values range of year
          yearValuesST = dataST_sp$minYear : dataST_sp$maxYear
          
          if(uncertainty){
            # Make slopes from yearly variations
            STSlopes = makeSlopesFromVar(dataVar_sp, N = N, limits = yearValuesST, weight = weight, group =  "ST")
            
            # Aesthetics
            lwdST = .8
            alphaST = .1
          }else{
            # Values of abundance for the whole years range
            estST = dataST_sp$estimate * (yearValuesST - min(yearValuesST))
            
            # Center on the middle value of the yearly variations
            midVal = mean(dataVar_sp$estimate[dataVar_sp$year >= min(yearValuesST)])
            estST = estST - mean(estST) + midVal
            
            # Turn to exponential if required
            estST = sapply(estST, function(x) ifelse(distribution == "gaussian", x, exp(x)))
            # Make dataframe
            STSlopes = data.frame(year = yearValuesST,
                                  estimate = estST,
                                  infIC = NA,
                                  supIC = NA,
                                  group = "ST",
                                  rep = 1)
            
            # Aesthetics
            lwdST = 1.1
            alphaST = 1
          }
          
        }else{
          plotST = FALSE
        }
      }else{
        plotST = FALSE
      }
      
      
      ###################
      # LONG-TERM TREND #
      ###################
      if(uncertainty){
        LTSlopes = makeSlopesFromVar(dataVar_sp, N = N, limits = yearValues, weight = weight, group =  "LT")
        # Aesthetics
        lwdLT = .8
        alphaLT = .1
        
      }else{
        # Values of abundance for the whole years range
        estLT = dataLT_sp$estimate * (yearValues - min(yearValues))
        
        # Center on middle value
        estLT = estLT - mean(estLT)
        
        # Turn to exponential if required
        estLT = sapply(estLT, function(x) ifelse(distribution == "gaussian", x, exp(x)))
        
        # Make the data.frame
        LTSlopes = data.frame(year = yearValues,
                              estimate = estLT,
                              infIC = NA,
                              supIC = NA,
                              group = "LT",
                              rep = 1)
        
        # Aesthetics
        lwdLT = 1.1
        alphaLT = 1
        
      }
      
      #####################
      # YEARLY VARIATIONS #
      #####################
      # Centrer sur le milieu de la courbe
      dataVar_sp[,c("estimate", "infIC", "supIC")] = dataVar_sp[,c("estimate", "infIC", "supIC")] - mean(dataVar_sp$estimate)
      
      # Passer à l'exponentiel, si données d'occurrence ou de comptage
      if(distribution != "gaussian"){
        dataVar_sp[,c("estimate", "infIC", "supIC")] = exp(dataVar_sp[,c("estimate", "infIC", "supIC")])
      }
      
      # Garder les colonnes d'intérêt
      dataVar_sp = dataVar_sp[,c("year","estimate", "infIC", "supIC")]
      
      # Ajouter les
      dataVar_sp$group = "VAR"
      dataVar_sp$rep = 1
      
      # FINAL DATAFRAME #
      dataPlot = rbind(LTSlopes, dataVar_sp)
      
      if(plotST){
        dataPlot = rbind(dataPlot, STSlopes) 
      }
      if(plotGamm){
        dataPlot = rbind(dataPlot, dataGammVariations_sp)
      }
      
      ###################################
      #   TITLE, SUBTITLES & AXIS NAME  #
      ###################################
      
      if(is.null(dataName)){
        title <- sp
      }else{
        # Filter french names for considered species
        dataName_sp = dataName[dataName$species == sp,]
        sc_name = paste0("(",dataName_sp$scientific_name, ")")
        
        if(nrow(dataName_sp) >0){
          fr_name = dataName_sp$french_name
          title <- bquote(.(fr_name)~italic(.(sc_name)))
          if(language == "english"){
            en_name = dataName_sp$english_name
            title <- bquote(.(en_name)~italic(.(sc_name)))
          }
        }else{
          title <- sp
        }
        
      }
      
      # Subtitle
      legendLT = percEvolutionToText_bis(data = dataLT_sp, distribution = distribution, language = language)
      
      if(plotST){
        legendST = percEvolutionToText_bis(data = dataST_sp, distribution = distribution, language = language)
      }
      
      legendGAM = ifelse(language == "french", "Tendance non linéaire", "Non-linear trend")
      legendVar = ifelse(language == "french", "Estimations annuelles", "Annual estimates")
      
      # Name of y axis
      yName = ifelse(language == "french", "Abondance relative", "Relative abundance")
      
      ######################
      #   REFERENCE LINE   # 
      ######################
      
      # Identify reference intercept
      yRef = ifelse(distribution == "gaussian", 0, 1)
      
      ####################
      #   SIGNIFICANCE   #
      ####################
      
      # Turn sinificance of the trend to text ----
      signLT <- significanceToText(dataLT_sp$pval)
      
      # x-position of the significance
      xSignLT <- dataLT_sp$maxYear + 1 ; 
      
      # y-position of the significance
      ySignLT <- mean(LTSlopes$estimate[LTSlopes$year == dataLT_sp$maxYear]) 
      
      # If ST trend exists, do the same ----
      if(plotST){
        xSignST <- xSignLT
        signST <- significanceToText(dataST_sp$pval)
        ySignST <- mean(STSlopes$estimate[STSlopes$year == dataST_sp$maxYear ]) 
      }
      
      ###################
      #   GROWTH RATE   #
      ###################
      
      # Extract growth rate
      grLT = ifelse(distribution != "gaussian", exp(dataLT_sp$estimate), dataLT_sp$estimate)
      
      # Round it
      grLT = round(grLT, 3)
      if(grLT == 1){grLT = paste0(1, ".0")}
      
      # Extract x & y coordinate
      xGrLT = dataLT_sp$minYear
      yGrLT = mean(LTSlopes$estimate[LTSlopes$year == min(LTSlopes$year)]) + 0.15
      
      if(plotST){
        # Extract growth rate
        grST = ifelse(distribution != "gaussian", exp(dataST_sp$estimate), dataST_sp$estimate)
        # Round it
        grST = round(grST, 3)
        # Extract x & y coordinate
        xGrST = dataST_sp$minYear - 1.25
        yGrST = mean(STSlopes$estimate[STSlopes$year == min(STSlopes$year)])
        
        if(grST == 1){
          grST = paste0(grST, ".0")
          yGrST = yGrST - 0.15
          xGrST = xGrST + 1.25
          
            }
        
      }
      
      ##############
      #   COLORS   #
      ##############
      if (plotGamm & plotST){
        # Color values
        colValues = c( "#F8C457","#429EBD","#FF4359", "#F8C457")
        colValues_guide = c("#429EBD","#FF4359", "#F8C457", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendST, legendGAM, legendVar)
        
        # Other characteristics
        linetype = c(1, 1, 1, 0) ; shape = c(NA, NA, NA, 18)
        size = c(1.1, 1.1, 1.1, 3) ; alpha = c(1, 1, 1, 1)
        
      }else if(plotGamm & !plotST){
        # Color values
        colValues = c( "#F8C457","#429EBD", "#F8C457") 
        colValues_guide = c("#429EBD", "#F8C457", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendGAM, legendVar)
        
        # Other characteristics
        linetype = c(1, 1, 0) ; shape = c(NA, NA, 18)
        size = c(1.1, 1.1, 3) ; alpha = c(1, 1, 1)
        
      }else if(!plotGamm & plotST){
        # Color values
        colValues = colValues_guide = c("#429EBD","#FF4359", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendST, legendVar) 
        
        # Other characteristics
        linetype = c(1, 1, 0) ; shape = c(NA, NA, 18)
        size = c(1.1, 1.1, 3) ; alpha = c(1, 1, 1)
        
      }else{
        # Color values
        colValues = colValues_guide = c("#429EBD", "#F8C457")
        
        # Labels for each color
        colLabels = c(legendLT, legendVar) 
        
        # Other characteristics
        linetype = c(1, 0) ; shape = c(NA, 18)
        size = c(1.1, 3) ; alpha = c(1, 1)
      }
      
      ############
      #   PLOT   #
      ############
      if(all(dataPlot$estimate < 10^10)){
        
        # Initialiser le plot 
        plot <- ggplot2::ggplot()
        
        # ADD GAMM INFORMATION
        if(plotGamm){
          plot <- plot +  
            # Add the confidence interval
            ggplot2:: geom_ribbon(data = dataPlot[dataPlot$group == "GAM",], ggplot2::aes(x = year, ymin = infIC, ymax = supIC), fill = "#F8C457", alpha = .2) +
            
            # Add the curve associated with the GAM
            ggplot2::geom_line(data = dataPlot[dataPlot$group == "GAM",], ggplot2::aes(x = year, y = estimate, col = group), linewidth = 1.05)
          
        }
        
        # ADD YEARLY VARIATIONS & LONG-TERM TREND
        plot = plot +
          
          # Add yearlyVariations
          ggplot2::geom_point(data = dataPlot[dataPlot$group == "VAR",], ggplot2::aes(x = year, y = estimate, col = group), shape = 18, size = 2.5) + 
          
          # Add the long-term trend
          ggplot2::geom_line(data = dataPlot[dataPlot$group == "LT",], ggplot2::aes(x = year, y = estimate, group = rep, col = group), 
                             alpha = alphaLT, linewidth = lwdLT) + 
          
          # Add LT significance information
          ggplot2::annotate(geom = "text", label = signLT, color = "#429EBD", size = 8, x = xSignLT, y = ySignLT) +
          
          # Add LT growth rate
          ggplot2::annotate(geom = "text", label = grLT, color = "#429EBD", size = 8, x = xGrLT, y = yGrLT)
        
        # ADD SHORT-TERM TREND
        if(plotST){
          plot = plot +
            
            # Add the short-term trend
            ggplot2::geom_line(data = dataPlot[dataPlot$group == "ST",], ggplot2::aes(x = year, y = estimate, group = rep, col = group), 
                               alpha = alphaST, linewidth = lwdST) + 
            
            # Add LT significance information
            ggplot2::annotate(geom = "text", label = signST, color = "#FF4359", size = 8, x = xSignST, y = ySignST) +
            
            # Add LT growth rate
            ggplot2::annotate(geom = "text", label = grST, color = "#FF4359", size = 8, x = xGrST, y = yGrST)
          
          
        }
        
        # FORMAT THE PLOT : AXIS, TITLES, ...
        plot <-  plot + 
          
          # Add the reference line
          ggplot2::geom_hline(yintercept = yRef, col = "gray", linetype = "dashed", linewidth = 1.01) +
          
          # Title
          ggplot2::ggtitle(title) +
          
          # Format x and y-axis
          ggplot2::scale_x_continuous(ifelse(language == "french", "Année", "Year"), 
                                      breaks = yearValues, limits = c(dataLT_sp$minYear - 1, dataLT_sp$maxYear +1)) + 
          
          # Rename y-axis                            
          ggplot2::ylab(yName) +
          
          # Change theme
          ggplot2::theme_bw() +
          
          # Adapt y-axis range
          ggplot2::coord_cartesian(ylim =  c(0.5, 2)) +  
          
          # Format theme
          ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, vjust = 1, hjust = 1), 
                         legend.text = ggplot2::element_text(size = 17),
                         axis.text = ggplot2::element_text(size = 17),
                         axis.title = ggplot2::element_text(size = 18),
                         title = ggplot2::element_text(size = 21), 
                         panel.grid.minor = ggplot2::element_blank(),
                         legend.position = "top", 
                         legend.justification = "left") +
          
          # Format type of trend colour legend
          ggplot2::scale_colour_manual(name = "", 
                                       values = colValues, 
                                       labels = colLabels, 
                                       guide = ggplot2::guide_legend(nrow = 2, 
                                                                     by.row = FALSE,
                                                                     override.aes = list(linetype = linetype,
                                                                                         color = colValues_guide,
                                                                                         shape = shape,
                                                                                         size = size,
                                                                                         alpha = alpha))) 
      } else{
        plot = ggplot2::ggplot() + ggplot2::theme_bw()
        cat("No plot for species ", sp," due to extreme values in some coefficients\n")
      }
      
    }
    
    # Save results
    if(save){
      ggplot2::ggsave(plot = plot, device = "png", filename = paste0(sp, ".png"), 
                      path = path, width = 25, height = 18, units = "cm")
    }
    
    return(plot) 
  }
  )
  
  return(plots)
}

# DATA IMPORT & FORMATTING ----
# Load data
data = read.csv("data/STOC_2023/countingData.csv")

# Extract the list of species
speciesList = unique(data$species)

# Create a unique ID for year/site
data$ID = paste0(data$year, data$site)

# Add 0's to the data.frame
data = fillAbsence(data[,c("ID","species", "year", "site", "abondance")], interestVar = "abondance", speciesList, "once")

# Make a dataframe with correspondency of species ID, scientific and english name
dataName = data.frame(species = c("PHYSIB", "CARCAR", "CORCOR", "CETCET"),
                      scientific_name = c("Phylloscopus sibilatrix", "Carduelis carduelis",
                                          "Corvus corone", "Cettia cetti"),
                      english_name = c("Wood warbler", "European goldfinch",
                                       "Carrion crow", "Cetti's warbler"))

# Format trends
speciesList = c("PHYSIB", "CARCAR", "CORCOR", "CETCET")
dataLongTermTrend = formatTrendEstimates(data, speciesList, repo = "STOC_2023", "longTermTrend")
dataYearlyVariations = formatYearlyVariations(speciesList, repo = "STOC_2023", "yearlyVariations", contr = "mean")
dataShortTermTrend = formatTrendEstimates(data[data$year >= max(data$year)-10,], speciesList, repo = "STOC_2023", "shortTermTrend")

# Make figures 
repo = "STOC_2023"
plotLinearTrends_bis(speciesList, dataLongTermTrend, dataYearlyVariations, dataShortTermTrend,
                     language = "english", dataName = dataName, distribution = "nbinom2", plotGamm = TRUE, plotST = TRUE,
                     uncertainty = TRUE, N = 100, weight = TRUE, save = TRUE, path = result_path)
