##Figure 2 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)
})

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

#Bioclimatic-zone specific temporal responses
species_temp = c('Moth', 'Butterfly', 'Bird', 'Plant', 'Phytopl.')
data_temp = data.frame(Species = factor(rep(species_temp, each = 5*4), levels = c('Plant', 'Phytopl.', 'Butterfly', 'Moth', 'Bird')),
                       Response = factor(rep(rep(c('CTI', 'CTIdiv' ,'CTBI', 'CTBIdiv'), each = 5), 5), levels = c('CTI', 'CTIdiv' ,'CTBI', 'CTBIdiv')),
                       Area = factor(rep(rep(c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal'), 4), 5),
                                     levels = c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal')),
                       Mean = NA, Min = NA, Max = NA)

for (i in 1:length(species_temp)) {
  #CTI, thermal range and thermal diversity
  # Mean
  resp_temp = c(models[[i]][[1]]$summary.lincomb.derived[1:5,2],
                models[[i]][[2]]$summary.lincomb.derived[1:5,2],
                models[[i]][[3]]$summary.lincomb.derived[1:5,2],
                models[[i]][[4]]$summary.lincomb.derived[1:5,2])
  data_temp$Mean[data_temp$Species==species_temp[i]] = resp_temp
  
  # Min
  resp_temp = c(models[[i]][[1]]$summary.lincomb.derived[1:5,4],
                models[[i]][[2]]$summary.lincomb.derived[1:5,4],
                models[[i]][[3]]$summary.lincomb.derived[1:5,4],
                models[[i]][[4]]$summary.lincomb.derived[1:5,4])
  data_temp$Min[data_temp$Species==species_temp[i]] = resp_temp
  
  # Max
  resp_temp = c(models[[i]][[1]]$summary.lincomb.derived[1:5,6],
                models[[i]][[2]]$summary.lincomb.derived[1:5,6],
                models[[i]][[3]]$summary.lincomb.derived[1:5,6],
                models[[i]][[4]]$summary.lincomb.derived[1:5,6])
  data_temp$Max[data_temp$Species==species_temp[i]] = resp_temp
  
}

#panel letters
graphLabels = data.frame(Response = factor(c("CTI", "CTIdiv", "CTBI", "CTBIdiv"), levels = c("CTI", "CTIdiv", "CTBI", "CTBIdiv")),
                         Species = .8, Mean = c(0.0175, 0.0009, 0.00275, 0.00125),
                          Letter = c("A)", "B)", "C)", "D)"))

g1 = ggplot(data = data_temp, aes(x = Species, y = Mean, fill = Area)) +
  geom_bar(stat = 'identity', position=position_dodge(.8), width=.8) +
  geom_point(position=position_dodge(.8), size = .6) +
  geom_errorbar(aes(ymin = Min, ymax = Max), position=position_dodge(.8), width = .5) +
  geom_text(data = graphLabels, aes(x = Species, y = Mean, label = Letter), inherit.aes = FALSE) +
  theme_classic() +
  geom_hline(yintercept = 0) +
  facet_wrap(~Response, scales = "free") +
  ylab(expression(paste(Delta, ' °C / year'))) +
  xlab('') +
  scale_fill_moma_d("VanGogh", name = '') +
  theme(legend.position = 'bottom', legend.text=element_text(size=6),
        axis.text.x = element_text(angle = 45, vjust = .5, hjust=.5))

#coordinate system EPSG: 3067
#get a background map
world = ne_countries(scale = "medium", returnclass = "sf")
world_plane = st_transform(world, 3067)

# crs_lea = "+proj=laea +lat_0=60 +lon_0=25 +x_0=-359705.51 +y_0=6615192.09 +datum=WGS84 +units=m +no_defs"
world_lea = st_transform(world, 3035)

#bioclimatic zones
bioclim = st_read('data/bioclimatic_zones/BioclimaticZonesFinland.shp')
bioclim$PaajakoNro = factor(c('Middle boreal', 'Northern boreal', 'Hemiboreal', 'Southern boreal'),
                            levels = c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal'))


g2 = ggplot() +
  geom_sf(data = world_lea, fill = "antiquewhite3", alpha = .5) +
  coord_sf(xlim = c(1e6, 7e6), ylim = c(0, 6e6), expand = FALSE) +
  annotate("rect", xmin = 4.6e6, xmax = 5.5e6, ymin = 4e6, ymax = 5.5e6, alpha = 0, color = "black", lwd = 1) +
  theme(legend.position = 'none', axis.text = element_text(size = 6),
        plot.title = element_text(size = 8)) +
  theme(panel.background = element_rect(fill = "azure"),
        panel.border = element_rect(fill = NA),
        panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", linewidth = 0.5)) +
  ggtitle('Study area')

g3 = ggplot() +
  geom_sf(data = world_plane, fill = "antiquewhite3", alpha = .5) +
  geom_sf(data = bioclim, aes(fill = PaajakoNro), alpha = .4) +
  coord_sf(xlim = c(0, 8e5), ylim = c(6.55e6, 7.8e6), expand = FALSE) +
  scale_fill_manual(values = moma.colors("VanGogh", n = 5, type = 'discrete')[2:5]) +
  ggtitle('Bioclimatic zones') +
  theme(legend.position = 'none', axis.text = element_text(size = 6),
        plot.title = element_text(size = 8)) +
  theme(panel.background = element_rect(fill = "azure"),
        panel.border = element_rect(fill = NA),
        panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", linewidth = 0.5))

cm2in = 0.393700787

png(paste0('Figures/Fig2.png'), width = (210-2*2)*.8, height = (297-2.54*2)/2, res = 600, units = 'mm')
g_layout = rbind(c(1,2),c(1,3))
grid.arrange(g1, g2, g3, layout_matrix = g_layout, widths = c(.8,.2))
dev.off()

tiff(paste0('Figures/Fig2.tif'), width = (210-2*2)*.8, height = (297-2.54*2)/2, res = 200, units = 'mm')
g_layout = rbind(c(1,2),c(1,3))
grid.arrange(g1, g2, g3, layout_matrix = g_layout, widths = c(.8,.2))
dev.off()

pdf(paste0('Figures/Fig2.pdf'), width = 16*cm2in, height = 14*cm2in)
g_layout = rbind(c(1,2),c(1,3))
grid.arrange(g1, g2, g3, layout_matrix = g_layout, widths = c(.8,.2))
dev.off()
