library(rgdal)
library(raster)
library(sf)
library(ggplot2)
library(sp)
library(stringr)
library(ggpubr)
library(colorspace)
library(tidyr)
library(viridis)
library(scales)
library(patchwork)
library(ggh4x)
library(ggridges)
library(cowplot)
library(gridExtra)

# # overall summary plot -  could be done with Hp or WL saved because we only want the All one
# # Normal ...........................................
# Fin data .............................................
setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/4.Predictions/3.CommunityMeasures")
load('PAtrendWL_FIN_2808.RDATA')
PArich_FINi = PArich[, 4:6] ## we look at the global all species only Richall
PArich_FINi$X = coef_trend_FIN$X # same site order
PArich_FINi$Y = coef_trend_FIN$Y # same site order
coef_trend_FINi = coef_trend_FIN[, 4:7]
Resid2_site_FINi = Resid2_site_FIN[, 4:7]

load('ABtrendWL_FIN_2808.RDATA')
AB_FINi = ABeve[, 4:6]
AB_FINi$X = coef_trend_FIN$X # same site order
AB_FINi$Y = coef_trend_FIN$Y # same site order
coef_trend_FINbi = coef_trend_FIN[, 4:7]
Resid2_site_FINbi = Resid2_site_FIN[, 4:7]

######## 
# 1 # Mean across years
# Prepare mean  .................................................................
colnames(PArich_FINi)[1] = c('All')
colnames(AB_FINi)[1] = c('All')

PArichmean_FINi = PArich_FINi %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

PArichmean_FINi$spatial = 'FINLAND'

ABmean_FINi = AB_FINi %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

ABmean_FINi$spatial = 'FINLAND'


#### # 2 # trend 


#### # 3 # Resid 


##.....................  
## Mean / trend / Variance


MaprichAM2_i <- ggplot(data = PArichmean_FINi, aes(x=X_mean, y=Y_mean, color=All_mean)) + geom_point(size=0.6) + 
  ggtitle("Average") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


MaprichAT2_i <- ggplot(data = coef_trend_FINi, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.6) + 
  ggtitle("Trend") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = colorspace::diverging_hcl(11, h = c(180, 50), c = 80, l = c(20, 95), power = c(0.7, 1.3)),
                        limits=c(-1, 1)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank())+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


##


MaprichAU2_i <- ggplot(data = Resid2_site_FINi, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.6) + 
  ggtitle("Residuals variance") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


## - Evenness ..................................................................
##-...................................................................................................

## Mean / trend / Variance


MapABAM2_i <- ggplot(data = ABmean_FINi, aes(x=X_mean, y=Y_mean, color=All_mean)) + geom_point(size=0.6) + 
  ggtitle("Average") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 



MapABAT2_i <- ggplot(data = coef_trend_FINbi, aes(x=X, y=Y, color=eveall)) + geom_point(size=0.6) + 
  ggtitle("Trend") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = colorspace::diverging_hcl(11, h = c(180, 50), c = 80, l = c(20, 95), power = c(0.7, 1.3)),
                        limits=c(-0.02, 0.02)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank())+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 

##

MapABAU2_i <- ggplot(data = Resid2_site_FINbi, aes(x=X, y=Y, color=eveall)) + geom_point(size=0.6) + 
  ggtitle("Residuals variance") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 



# 
# # Habsc...........................................
# Fin data .............................................
load('PAtrendWL_FIN_2808hab.RDATA')
PArich_FINhab = PArich[, 4:6]
PArich_FINhab$X = coef_trend_FIN$X # same site order
PArich_FINhab$Y = coef_trend_FIN$Y # same site order
coef_trend_FINhab = coef_trend_FIN[, 4:7]
Resid2_site_FINhab = Resid2_site_FIN[, 4:7]

load('ABtrendWL_FIN_2808hab.RDATA')
AB_FINhab = ABeve[, 4:6]
AB_FINhab$X = coef_trend_FIN$X # same site order
AB_FINhab$Y = coef_trend_FIN$Y # same site order
coef_trend_FINbhab = coef_trend_FIN[, 4:7]
Resid2_site_FINbhab = Resid2_site_FIN[, 4:7]

######## 
# 1 # Mean across years
# Prepare mean  .................................................................
colnames(PArich_FINhab)[1] = c('All')
colnames(AB_FINhab)[1] = c('All')

PArichmean_FINhab = PArich_FINhab %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

PArichmean_FINhab$spatial = 'FINLAND'

ABmean_FINhab = AB_FINhab %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

ABmean_FINhab$spatial = 'FINLAND'


#### # 2 # trend 
#### # 3 # Resid 

##.....................  
## Mean / trend / Variance

MaprichAM2_hab <- ggplot(data = PArichmean_FINhab, aes(x=X_mean, y=Y_mean, color=All_mean)) + geom_point(size=0.6) + 
  ggtitle("Average") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


MaprichAT2_hab <- ggplot(data = coef_trend_FINhab, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.6) + 
  ggtitle("Trend") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = colorspace::diverging_hcl(11, h = c(180, 50), c = 80, l = c(20, 95), power = c(0.7, 1.3)),
                        limits=c(-1, 1)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank())+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


##


MaprichAU2_hab <- ggplot(data = Resid2_site_FINhab, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.6) + 
  ggtitle("Residuals variance") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


## - Evenness ..................................................................
##-...................................................................................................

## Mean / trend / Variance


MapABAM2_hab <- ggplot(data = ABmean_FINhab, aes(x=X_mean, y=Y_mean, color=All_mean)) + geom_point(size=0.6) + 
  ggtitle("Average") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 



MapABAT2_hab <- ggplot(data = coef_trend_FINbhab, aes(x=X, y=Y, color=eveall)) + geom_point(size=0.6) + 
  ggtitle("Trend") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = colorspace::diverging_hcl(11, h = c(180, 50), c = 80, l = c(20, 95), power = c(0.7, 1.3)),
                        limits=c(-0.02, 0.02)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank())+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 

##

MapABAU2_hab <- ggplot(data = Resid2_site_FINbhab, aes(x=X, y=Y, color=eveall)) + geom_point(size=0.6) + 
  ggtitle("Residuals variance") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


# 
# # Climsc...........................................
load('PAtrendWL_FIN_2808clim.RDATA')
PArich_FINclim = PArich[, 4:6]
PArich_FINclim$X = coef_trend_FIN$X # same site order
PArich_FINclim$Y = coef_trend_FIN$Y # same site order
coef_trend_FINclim = coef_trend_FIN[, 4:7]
Resid2_site_FINclim = Resid2_site_FIN[, 4:7]

load('ABtrendWL_FIN_2808clim.RDATA')
AB_FINclim = ABeve[, 4:6]
AB_FINclim$X = coef_trend_FIN$X # same site order
AB_FINclim$Y = coef_trend_FIN$Y # same site order
coef_trend_FINbclim = coef_trend_FIN[, 4:7]
Resid2_site_FINbclim = Resid2_site_FIN[, 4:7]

######## 
# 1 # Mean across years
# Prepare mean  .................................................................
colnames(PArich_FINclim)[1] = c('All')
colnames(AB_FINclim)[1] = c('All')

PArichmean_FINclim = PArich_FINclim %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

PArichmean_FINclim$spatial = 'FINLAND'

ABmean_FINclim = AB_FINclim %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

ABmean_FINclim$spatial = 'FINLAND'


#### # 2 # trend 
#### # 3 # Resid 

##.....................  
## Mean / trend / Variance

MaprichAM2_clim <- ggplot(data = PArichmean_FINclim, aes(x=X_mean, y=Y_mean, color=All_mean)) + geom_point(size=0.6) + 
  ggtitle("Average") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


MaprichAT2_clim <- ggplot(data = coef_trend_FINclim, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.6) + 
  ggtitle("Trend") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = colorspace::diverging_hcl(11, h = c(180, 50), c = 80, l = c(20, 95), power = c(0.7, 1.3)),
                        limits=c(-1, 1)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank())+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


##


MaprichAU2_clim <- ggplot(data = Resid2_site_FINclim, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.6) + 
  ggtitle("Residuals variance") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


## - Evenness ..................................................................
##-...................................................................................................

## Mean / trend / Variance


MapABAM2_clim <- ggplot(data = ABmean_FINclim, aes(x=X_mean, y=Y_mean, color=All_mean)) + geom_point(size=0.6) + 
  ggtitle("Average") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 



MapABAT2_clim <- ggplot(data = coef_trend_FINbclim, aes(x=X, y=Y, color=eveall)) + geom_point(size=0.6) + 
  ggtitle("Trend") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = colorspace::diverging_hcl(11, h = c(180, 50), c = 80, l = c(20, 95), power = c(0.7, 1.3)),
                        limits=c(-0.02, 0.02)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank())+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 

##

MapABAU2_clim <- ggplot(data = Resid2_site_FINbclim, aes(x=X, y=Y, color=eveall)) + geom_point(size=0.6) + 
  ggtitle("Residuals variance") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "bottom",plot.title = element_text(size=9), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 8),
        legend.spacing.x = unit(0.5, 'cm'), legend.spacing.y = unit(0.3, 'cm'),
        legend.title = element_text(size = 0), axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks = element_blank(), 
        strip.text = element_text(size=9), strip.background = element_blank()) +
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 7,
                                barheight = 0.5)) 


## prepare plots
#..........................................................................

# Arranging the plot using cowplot
HistPart = plot_grid(MaprichAM2_i, MaprichAT2_i, MaprichAU2_i, NULL,
                    MapABAM2_i, MapABAT2_i, MapABAU2_i, 
                    ncol = 7, nrow=1, align = "hv", labels=c('a', 'b', 'c', '', 'd', 'e', 'f'),
                    rel_widths = c(1,1,1, 0.25, 1,1,1))

ClimPart = plot_grid(MaprichAM2_clim, MaprichAT2_clim, MaprichAU2_clim, NULL,
                     MapABAM2_clim, MapABAT2_clim, MapABAU2_clim, 
                     ncol = 7, nrow=1, align = "hv", labels=c('g', 'h', 'i', '', 'j', 'k', 'l'),
                     rel_widths = c(1,1,1, 0.25, 1,1,1))

HabPart = plot_grid(MaprichAM2_hab, MaprichAT2_hab, MaprichAU2_hab, NULL,
                     MapABAM2_hab, MapABAT2_hab, MapABAU2_hab, 
                     ncol = 7, nrow=1, align = "hv", labels=c('m', 'n', 'o', '', 'p', 'q', 'r'),
                    rel_widths = c(1,1,1, 0.25, 1,1,1))

# grid.arrange(arrangeGrob(HistPart, left = "Historical"), 
#              NULL,
#              arrangeGrob(ClimPart, left = "Climate scenario"), 
#              arrangeGrob(HabPart, left = "Habitat scenario"), 
#              ncol = 1, heights = c(1, 0.15, 1,1), top=' Richness                                                                                                      
#              Evenness')

library(grid)

blank <- grid.rect(gp=gpar(col="white"))

grid.arrange(arrangeGrob(HistPart, left = "Climate & Habitat"), 
             blank, ## NULL does not work anymore
             arrangeGrob(ClimPart, left = "Climate scenario"), 
             arrangeGrob(HabPart, left = "Habitat scenario"), 
             ncol = 1, heights = c(1, 0.1, 1,1), top='   Richness                                                                                                                  Evenness') + 
  annotation_custom(grid.polygon(c(0, 0.5, 1, 0.51, 0.51, 0.51),
                 c(0.65, 0.65, 0.65,0, 0.5, 1), 
                 id = c(1,1,1,2,2,2), 
                 gp = gpar(lwd = 2.5)))

HistPart <- annotate_figure(HistPart, left = text_grob("Climate & habitat", face = "bold", rot=90, size=12))
ClimPart <- annotate_figure(ClimPart, left = text_grob("Climate change only scenario", face = "bold", rot=90, size=12))
HabPart <- annotate_figure(HabPart, left = text_grob("Habitat change only scenario", face = "bold", rot=90, size=12))


grid.arrange(arrangeGrob(HistPart), 
             blank,
             arrangeGrob(ClimPart), 
             arrangeGrob(HabPart), 
             ncol = 1, heights = c(1, 0.1, 1,1), top='   Richness                                                                                                                  Evenness') + 
  annotation_custom(grid.polygon(c(0, 0.5, 1, 0.51, 0.51, 0.51),
                                 c(0.65, 0.65, 0.65,0, 0.5, 1), 
                                 id = c(1,1,1,2,2,2), 
                                 gp = gpar(lwd = 2.5)))

####
# dissaemble
MaprichAM2_i2 <- ggplot_build(MaprichAM2_i)

# changes
MaprichAM2_i2$plot$labels$title = 'Average' 
MaprichAM2_i2$plot

