## Fig_S1 in supporting information

#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('CTImean', 'CTIdiv', 'STBImean', 'STBIdiv')
species_title = c('Moth', 'Butterfly', 'Bird', 'Plant', 'Phytoplankton')


##Distribution of sampling sites
#loop over species
g1 = list()
i_sp = 0

for (species_group in c('moth', 'butterfly', 'bird', 'plant', 'phytoplankton')) {
  i_sp = i_sp + 1
  data_fig = data_all[[i_sp]]
  coord = data.frame(unique(cbind(data_fig$Longitude, data_fig$Latitude)))
  colnames(coord) = c('X', 'Y')
  
  g1[[i_sp]] = ggplot() +
    geom_sf(data = world_plane, fill = "antiquewhite3", alpha = .5) +
    geom_point(data=coord, aes(x = X, y = Y), size = .5) +
    coord_sf(xlim = c(0, 8e5), ylim = c(6.55e6, 7.8e6), expand = FALSE) +
    theme_minimal() +
    xlab('') +
    ylab('') +
    ggtitle(species_title[i_sp]) +
    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))
  
  
  if (i_sp != 4) {
    g1[[i_sp]] = g1[[i_sp]] +
      theme(axis.text.x = element_blank(), axis.text.y = element_blank())
  }
}

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


#Scatter plots
for (s in 1:length(species_title)) {
  ind = 1:(length(models[[s]][[1]]$summary.random$SiteID_tr_2$mean)/2)
  for (i in 1:(length(response_all)-1)) {
    for (j in (i+1):length(response_all)) {
      ind_plot = ind_plot+1
      if (i == 1 & j == 2) {
        temp = data.frame(x = models[[s]][[i]]$summary.random$SiteID_tr_2$mean[ind],
                          y = models[[s]][[j]]$summary.random$SiteID_tr_2$mean[ind],
                          xlab = response_all[i],
                          ylab = response_all[j])
      } else {
        temp = rbind(temp,
                     data.frame(x = models[[s]][[i]]$summary.random$SiteID_tr_2$mean[ind],
                          y = models[[s]][[j]]$summary.random$SiteID_tr_2$mean[ind],
                          xlab = response_all[i],
                          ylab = response_all[j]))
        
      }
    }
  }
  
  temp$xlab = factor(temp$xlab, levels = response_all)
  temp$ylab = factor(temp$ylab, levels = response_all)
  
  g2 = ggplot(data = temp) +
    geom_point(aes(x=x, y=y)) +
    geom_smooth(aes(x=x, y=y), method = 'lm') +
    theme_bw() +
    facet_grid(ylab~xlab, scales = 'free') +
    ggtitle(species_title[s]) +
    xlab('') +
    ylab('') +
    theme(axis.text.x = element_text(angle = 60, vjust = .5, hjust=.5))
    
  
  # png(paste0('Figures/Manuscript/Scatter_plots_', species_title[s], '.png'), width = (210-2*2)*.7, height = (297-2.54*2)/2, res = 600, units = 'mm')
  # g2
  # dev.off()

  ggsave(plot = g2, filename = paste0('Figures/Manuscript/Scatter_plots_', species_title[s], '.png'), width = (210-2*2)*.7, height = (297-2.54*2)/2, dpi = 600, units = 'mm')
  
}


