##Figure 4 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 models and data
load('data/Model_fits.RData')

#coordinate system EPSG: 3067
#get a background map
world = ne_countries(scale = "medium", returnclass = "sf")
world_plane = st_transform(world, 3067)
#bioclimatic zones
bioclim = st_read('data/bioclimatic_zones/BioclimaticZonesFinland.shp')

#response variables
response_all = c('CTI', 'CTIdiv', 'STBImean', 'STBIdiv')
species_title = c('Moth', 'Butterfly', 'Bird', 'Forest plant', 'Phytoplankton')

g1 = list()
test_temp = list()
ind_plot = 0

for (s in 1:length(species_title)) {
  
  ind = 1:(length(models[[s]][[1]]$summary.random$SiteID_tr_2$mean)/2)
  ind_plot = ind_plot+1
      
  temp = data.frame(x = models[[s]][[1]]$summary.random$SiteID_tr_2$mean[ind],
                    y = models[[s]][[2]]$summary.random$SiteID_tr_2$mean[ind])
      
  g1[[ind_plot]] = ggplot(data = temp) +
    geom_point(aes(x=x, y=y)) +
    geom_smooth(aes(x=x, y=y), method = 'lm') +
    theme_classic() +
    ggtitle(species_title[s]) +
    scale_x_continuous(breaks = pretty(temp$x, n = 2)) +
    scale_y_continuous(breaks = pretty(temp$y, n = 2)) +
    theme(axis.text.x = element_text(angle = 45, vjust = .5, hjust=.5, size = 6),
          axis.text.y = element_text(angle = 65, vjust = .5, hjust=.5, size = 6),
          axis.title.x = element_text(size = 8),
          axis.title.y = element_text(size = 8),
          plot.title = element_text(size = 8)) +
    xlab('') +
    ylab('')
  
  test_temp[[s]] = cor.test(temp$x, temp$y)
  
}

g1[[4]] = g1[[4]] +
  xlab(expression(paste(Delta, ' CTI/year'))) +
  ylab(expression(paste(Delta, ' CTIdiv/year')))
  
cm2in = 0.393700787

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

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

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

