##Figure 5 of the manuscript

#set path
try(setwd(''))

suppressPackageStartupMessages({
  library(terra)
  library(tidyverse)
  library(vegan)
  library(reshape)
  library(MASS)
  library(sf)
  library(rnaturalearth)
  library(rnaturalearthdata)
  library(paletteer)
  library(ggrepel)
  library(gridExtra)
  library(RColorBrewer)
  library(INLA)
  library(pracma)
  library(viridis)
  library(ggsci)
  library(wesanderson)
  library(betapart)
  library(parallel)
  library(stars)
  library(adegenet)
  library(spdep)
  library(MoMAColors)
  library(tidyterra)
  library(png)
  library(grid)
})

##read in model structures and data
load('data/Model_fits.RData')

##silhouettes
img = readPNG("Figures/Silhouettes/Bird.png")
bird = rasterGrob(img, interpolate=TRUE)
img = readPNG("Figures/Silhouettes/Moth.png")
moth = rasterGrob(img, interpolate=TRUE)
img = readPNG("Figures/Silhouettes/Butterfly.png")
butterfly = rasterGrob(img, interpolate=TRUE)
img = readPNG("Figures/Silhouettes/Plant.png")
plant = rasterGrob(img, interpolate=TRUE)
img = readPNG("Figures/Silhouettes/Phytoplankton.png")
phytoplankton = rasterGrob(img, interpolate=TRUE)

g_all_4 = list()

i_sp = 0
species_title = c('Moth', 'Butterfly', 'Bird', 'Forest plant', 'Phytoplankton')

for (species_group in c('moth', 'butterfly', 'bird', 'plant', 'phytoplankton')) {
  i_sp = i_sp + 1
  
  data_temp = data_all[[i_sp]]
  
  temp_rast = rast(paste0('data/compiled_data/Temp_change_', species_group, '.tif'))
  coord = data.frame(unique(cbind(data_temp$Longitude, data_temp$Latitude)))
  colnames(coord) = c('X', 'Y')
  temp_change = terra::extract(temp_rast, coord)
  
  data_g1 = data.frame(Year_orig = unique(data_temp$Year_orig),
                       mean = models[[i_sp]][[1]]$summary.fixed$mean[1] + unique(data_temp$Year) * models[[i_sp]][[1]]$summary.lincomb.derived[1,2],
                       temp = unique(data_temp$Year) * mean(temp_change$mean, na.rm = T))
  
  data_g1$mean = data_g1$mean-mean(data_g1$mean)
  data_g1$temp = data_g1$temp-mean(data_g1$temp)
  data_g1$min1 = data_g1$mean - 2*(models[[i_sp]][[2]]$summary.fixed$mean[1] + unique(data_temp$Year)*models[[i_sp]][[2]]$summary.lincomb.derived[1,2])
  data_g1$max1 = data_g1$mean + 2*(models[[i_sp]][[2]]$summary.fixed$mean[1] + unique(data_temp$Year)*models[[i_sp]][[2]]$summary.lincomb.derived[1,2])
  
  g_all_4[[i_sp]] = ggplot(data = data_g1) +
    geom_line(aes(x=Year_orig, y=mean, col = 'CTI'), linewidth = 1) +
    geom_ribbon(aes(x=Year_orig, y=mean, ymin = min1, ymax = max1), alpha=0.1, col = MoMAPalettes$Rattner[[1]][4], show.legend = FALSE) + 
    geom_line(aes(x=Year_orig, y=temp, col = 'Temperature'), linewidth = 1) +
    scale_color_manual(values = c(MoMAPalettes$Rattner[[1]][4], MoMAPalettes$Rattner[[1]][5])) +
    ggtitle(species_group) +
    coord_cartesian(ylim = c(-1,1)) +
    ggtitle(species_title[i_sp]) +
    theme_classic() +
    theme(legend.position = 'none', 
          axis.text.x = element_text(angle = 60, vjust = .5, hjust=.5, size = 6),
          axis.text.y = element_text(size = 6),
          axis.title.x = element_text(size = 8),
          axis.title.y = element_text(size = 8),
          plot.title = element_text(size = 8)) +
    ylab('Air temperature,\nCTI anomaly (°C)') +
    xlab('')
    
}

g_all_4[[1]] = g_all_4[[1]] +
  theme(axis.text.y = element_blank()) +
  ylab('') +
  annotation_custom(moth, xmin=1993, xmax=1998, ymin=.65, ymax=1.1)

g_all_4[[2]] = g_all_4[[2]] +
  theme(axis.text.y = element_blank()) +
  ylab('') +
  annotation_custom(butterfly, xmin=2000, xmax=2005, ymin=.65, ymax=1.1)

g_all_4[[3]] = g_all_4[[3]] +
  theme(axis.text.y = element_blank()) +
  ylab('') +
  annotation_custom(bird, xmin=1980, xmax=1988, ymin=.65, ymax=1.1)

g_all_4[[4]] = g_all_4[[4]] +
  scale_x_continuous(labels=c('1985', '', '1995', '', '2005')) +
  annotation_custom(plant, xmin=1985, xmax=1990, ymin=.65, ymax=1.1)

g_all_4[[5]] = g_all_4[[5]] +
  theme(axis.text.y = element_blank()) +
  ylab('') +
  annotation_custom(phytoplankton, xmin=1978, xmax=1983, ymin=.65, ymax=1.1)

# png(paste0('Figures/Fig4.png'), width = (210-2*2)*.7, height = (297-2.54*2)/2.5, res = 600, units = 'mm')
# grid.arrange(g_all_4[[3]], g_all_4[[2]], g_all_4[[1]], g_all_4[[5]], g_all_4[[4]], nrow = 2)
# dev.off()

cm2in = 0.393700787

png(paste0('Figures/Fig5.png'), width = (210-2*2)*.8, height = (297-2.54*2)/6, res = 600, units = 'mm')
grid.arrange(g_all_4[[4]], g_all_4[[5]], g_all_4[[2]], g_all_4[[1]], g_all_4[[3]], nrow = 1, widths = c(1.3,1,1,1,1))
dev.off()

tiff(paste0('Figures/Fig5.tif'), width = (210-2*2)*.8, height = (297-2.54*2)/6, res = 200, units = 'mm')
grid.arrange(g_all_4[[4]], g_all_4[[5]], g_all_4[[2]], g_all_4[[1]], g_all_4[[3]], nrow = 1, widths = c(1.3,1,1,1,1))
dev.off()

pdf(paste0('Figures/Fig5.pdf'), width = 16*cm2in, height = 5*cm2in)
grid.arrange(g_all_4[[4]], g_all_4[[5]], g_all_4[[2]], g_all_4[[1]], g_all_4[[3]], nrow = 1, widths = c(1.3,1,1,1,1))
dev.off()
