##########################################################################################################
#
#		Plot map community diversity - empirical - Hostplant forms
#
##########################################################################################################


# spatial prep
library(rgdal)
library(raster)
library(sf)
library(sp)

# play with data
library(stringr)
library(tidyr)
library(scales)


## PLotting
library(ggplot2)
library(ggpubr)
library(colorspace)
library(viridis)
library(cowplot)
library(gridExtra)
library(grid)
library(lattice)

# Site data ...............................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions/5.CommunityMeasures_Sites")

load('ABtrendHP_SITE_2808.RDATA') ## 
AB_SITE = ABeve
coef_trend_SITEb = coef_trend_SITE
Resid_site_SITEb = Resid_site_SITE
Resid2_site_SITEb = Resid2_site_SITE


load('PAtrendHP_SITE_2808.RDATA')
PArich_SITE = PArich


setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/2.ObsInfo")
load('MothclusterFIN.RDATA')
dim(Moth_knn)  # 109    13
head(Moth_knn)

Moth_knn2 = Moth_knn[, 7:10]
colnames(Moth_knn2)[2] = 'sample'
head(Moth_knn2)

AB_SITE$X = Moth_knn2$XEUREF# same site order
AB_SITE$Y = Moth_knn2$YEUREF
coef_trend_SITEb$X = Moth_knn2$XEUREF# same site order
coef_trend_SITEb$Y = Moth_knn2$YEUREF

coef_trend_SITE$X = coef_trend_SITEb$X
coef_trend_SITE$Y = coef_trend_SITEb$Y

PArich_SITE$X = coef_trend_SITEb$X
PArich_SITE$Y = coef_trend_SITEb$Y

Resid_site_SITEb$X = coef_trend_SITEb$X
Resid_site_SITEb$Y = coef_trend_SITEb$Y
Resid2_site_SITEb$X = coef_trend_SITEb$X
Resid2_site_SITEb$Y = coef_trend_SITEb$Y

Resid_site_SITE$X = coef_trend_SITEb$X
Resid_site_SITE$Y = coef_trend_SITEb$Y
Resid2_site_SITE$X = coef_trend_SITEb$X
Resid2_site_SITE$Y = coef_trend_SITEb$Y

# Fin data .............................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions/1.CommunityMeasures")

load('ABtrendHP_FIN_2808.RDATA')
AB_FIN = ABeve
AB_FIN$X = coef_trend_FIN$X # same site order
AB_FIN$Y = coef_trend_FIN$Y # same site order

coef_trend_FINb = coef_trend_FIN
Resid_site_FINb = Resid_site_FIN
Resid2_site_FINb = Resid2_site_FIN



load('PAtrendHP_FIN_2808.RDATA')
PArich_FIN = PArich
PArich_FIN$X = coef_trend_FIN$X # same site order
PArich_FIN$Y = coef_trend_FIN$Y # same site order


PArich_FIN_test = PArich_FIN[,c(7:11)] # p

# recalculate RichAll
PArich_FIN_test$Richall = PArich_FIN$RichH1 + PArich_FIN$RichH2 + PArich_FIN$RichH4 +
  PArich_FIN$RichH31 + PArich_FIN$RichH32 + PArich_FIN$RichH33


#PArich_FIN_test = PArich_FIN_test[PArich_FIN_test$Year == 1,]

# to check.......................................
ggplot(data = PArich_FIN_test, aes(x=X, y=Y, color=Richall)) + geom_point(size=0.8) + 
  facet_wrap(~ Year, scales = "free") +
  #ggtitle("Richness") +
  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=12), axis.text=element_text(size=6),
        axis.title = element_text(size = 9), legend.text = element_text(size = 12),
        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()) +
  #geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick')+
  guides(color = guide_colorbar(title = "",
                                label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 15,
                                barheight = 1)) 

PArich_FIN$Richall = PArich_FIN$RichH1 + PArich_FIN$RichH2 + PArich_FIN$RichH4 +
  PArich_FIN$RichH31 + PArich_FIN$RichH32 + PArich_FIN$RichH33



# Prepare mean  .................................................................
colnames(PArich_FIN)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
colnames(AB_FIN)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')

PArichmean_FIN = PArich_FIN %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

PArichmean_FIN$spatial = 'FINLAND'

colnames(PArich_SITE)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
colnames(AB_SITE)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
PArichmean_SITE = PArich_SITE %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

PArichmean_SITE$spatial = 'SITE'

ABmean_FIN = AB_FIN %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

ABmean_FIN$spatial = 'FINLAND'


ABmean_SITE = AB_SITE %>%
  group_by(SiteID) %>%
  dplyr::summarise(across(everything(), list(mean = mean)))

ABmean_SITE$spatial = 'SITE'

Combined_PA = rbind(PArichmean_FIN, PArichmean_SITE)
Combined_AB = rbind(ABmean_FIN, ABmean_SITE)


###-----------------------------------------------------------------------------------------------------------------------------
########  Plot and change trends
# 1 # Mean across years

PArichmean_FIN$All = PArichmean_FIN$herbaceous_mean + PArichmean_FIN$graminoids_mean + PArichmean_FIN$others_mean +
  PArichmean_FIN$forbs_mean + PArichmean_FIN$deciduous_mean + PArichmean_FIN$coniferous_mean

PArichmean_SITE$All = PArichmean_SITE$herbaceous_mean + PArichmean_SITE$graminoids_mean + PArichmean_SITE$others_mean +
  PArichmean_SITE$forbs_mean + PArichmean_SITE$deciduous_mean + PArichmean_SITE$coniferous_mean


PArichmean_FIN = PArichmean_FIN[, c(1:7, 13, 9:12)]
PArichmean_SITE = PArichmean_SITE[, c(1:7, 13, 9:12)]

PArichmean_FINlong = gather(PArichmean_FIN, HP, S, herbaceous_mean:All, factor_key=TRUE)
ABmean_FINlong = gather(ABmean_FIN, HP, EVE, herbaceous_mean:eveall_mean, factor_key=TRUE)

PArichmean_SITElong = gather(PArichmean_SITE, HP, S, herbaceous_mean:All, factor_key=TRUE)
ABmean_SITElong = gather(ABmean_SITE, HP, EVE, herbaceous_mean:eveall_mean, factor_key=TRUE)

PArichmean_FINlong$Slog = log(PArichmean_FINlong$S)
PArichmean_SITElong$Slog = log(PArichmean_SITElong$S)
ABmean_FINlong$EVElog = log(ABmean_FINlong$EVE)
ABmean_SITElong$EVElog = log(ABmean_SITElong$EVE)

Combined_PA = rbind(PArichmean_FINlong, PArichmean_SITElong)
Combined_AB = rbind(ABmean_FINlong, ABmean_SITElong)


#####
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions")
load('Predsmall2408.RDATA')

load('CurrentSmallFIN_1408.RDATA')
tail(mapData2)
tail(mapData2hab)

mapDataALL = cbind(mapData2[, c(1:2, 6, 15, 10:11)], mapData2hab[, c(2:10)])
head(mapDataALL)
dim(mapDataALL)
mapcoord1 = mapDataALL[mapDataALL$Year == 1999, 14:15]
mapDataALL$SiteID = paste(mapDataALL$X, '', mapDataALL$Y)
SiteFIn_un = unique(mapDataALL$SiteID)
MapdataALL.list = split(mapDataALL, f = mapDataALL$Year)
Fin_coord = MapdataALL.list[[1]][, 14:15]


## prepare trend all again with new calculation
PArich_FIN$Year = as.integer(PArich_FIN$Year)

Trend_site_FIN = Resid_site_FIN = Resid2_site_FIN = coef_trend_FIN = matrix(NA, nrow= 12340, ncol=11)
for (i in 1:length(SiteFIn_un)) {
  MapdataALL_i = PArich_FIN[PArich_FIN$SiteID == SiteFIn_un[i],]

  for (j in colnames(MapdataALL_i)[1:7]) {
    index_j = match(j, colnames(MapdataALL_i)[1:7])
    Fit_site_FIN = lm(reformulate('Year', j), data = MapdataALL_i)

    Resid_site_FIN[i, index_j] = var(resid(Fit_site_FIN))
    Resid2_site_FIN[i, index_j] = var(rstandard(Fit_site_FIN), na.rm=TRUE)
    coef_trend_FIN[i, index_j] = Fit_site_FIN$coefficients[2]
  }
}

Resid_site_FIN = as.data.frame(Resid_site_FIN)
colnames(Resid_site_FIN)[1:7] = colnames(MapdataALL_i)[1:7]
Resid_site_FIN$X = Fin_coord$X
Resid_site_FIN$Y = Fin_coord$Y
Resid_site_FIN$spatial = 'FINLAND'

Resid2_site_FIN = as.data.frame(Resid2_site_FIN)
colnames(Resid2_site_FIN)[1:7] = colnames(MapdataALL_i)[1:7]
Resid2_site_FIN$X = Fin_coord$X
Resid2_site_FIN$Y = Fin_coord$Y
Resid2_site_FIN$spatial = 'FINLAND'

coef_trend_FIN = as.data.frame(coef_trend_FIN)
colnames(coef_trend_FIN)[1:7] = colnames(MapdataALL_i)[1:7]
coef_trend_FIN$X = Fin_coord$X
coef_trend_FIN$Y = Fin_coord$Y
coef_trend_FIN$spatial = 'FINLAND'


#### # 2 # trend 

PAtrend_FIN = coef_trend_FIN[, c(1:7, 12:14)]
colnames(PAtrend_FIN)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
PAtrend_SITE = coef_trend_SITE[, c(1:7, 9:10, 8)]
colnames(PAtrend_SITE)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')

PAtrend_FINlong = gather(PAtrend_FIN, HP, S, herbaceous:Richall, factor_key=TRUE)
PAtrend_SITElong = gather(PAtrend_SITE, HP, S, herbaceous:Richall, factor_key=TRUE)


PAtrend_FINlongA = PAtrend_FINlong[PAtrend_FINlong$HP == 'Richall',]
PAtrend_FINlongA$HP = 'All'




#### # 3 # Resid 

PAResid_FIN = Resid2_site_FIN[, c(1:7, 12:14)]
PAResid_SITE = Resid2_site_SITE[, c(1:7, 9:10, 8)]
colnames(PAResid_FIN)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
colnames(PAResid_SITE)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')


PAResid_FINlong = gather(PAResid_FIN, HP, S, herbaceous:Richall, factor_key=TRUE)
PAResid_SITElong = gather(PAResid_SITE, HP, S, herbaceous:Richall, factor_key=TRUE)


PAResid_FINlongA = PAResid_FINlong[PAResid_FINlong$HP == 'Richall',]
PAResid_FINlongA$HP = 'All'




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


PArichmean_FINlongA = PArichmean_FINlong[PArichmean_FINlong$HP == 'All',]

MaprichAM2 <- ggplot(data = PArichmean_FINlongA, aes(x=X_mean, y=Y_mean, color=S)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Richness mean") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 


PArichmean_FINlong2 = PArichmean_FINlong[!PArichmean_FINlong$HP == 'All',]
levels(PArichmean_FINlong2$HP) = c("herbaceous", "graminoids", "others", "forbs", "deciduous", "coniferous", "All")
levels(PArichmean_FINlong2$HP)[4] = 'forbs'
MaprichM2 <- ggplot(data = PArichmean_FINlong2, aes(x=X_mean, y=Y_mean, color=S)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Richness mean") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 



Combined_PA = rbind(PArichmean_FINlong, PArichmean_SITElong)
Combined_PAA = Combined_PA[Combined_PA$HP == 'All',]

RichAM2.gg = ggplot(Combined_PAA, aes(x=spatial, y=S, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3)) +# ylim(-0.5, 0.5) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

Combined_PA2 = Combined_PA[!Combined_PA$HP == 'All',]
levels(Combined_PA2$HP) = c("herbaceous", "graminoids", "others", "forbs", "deciduous", "coniferous", "All")
levels(Combined_PA2$HP)[4] = 'forbs'
RichM2.gg = ggplot(Combined_PA2, aes(x=spatial, y=S, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3)) +# ylim(-0.5, 0.5) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())




PAtrend_FINlongA = PAtrend_FINlong[PAtrend_FINlong$HP == 'Richall',]
PAtrend_FINlongA$HP = 'All'

MaprichAT2 <- ggplot(data = PAtrend_FINlongA, aes(x=X, y=Y, color=S)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Richness 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.5, 0.5), oob=squish) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank())+
  geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 


PAtrend_FINlong2 = PAtrend_FINlong[!PAtrend_FINlong$HP == 'Richall',]
levels(PAtrend_FINlong2$HP)[4] = 'forbs'
MaprichT2 <- ggplot(data = PAtrend_FINlong2, aes(x=X, y=Y, color=S)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Richness 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.5, 0.5)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 



Combined_PA = rbind(PAtrend_FINlong, PAtrend_SITElong)
Combined_PAA = Combined_PA[Combined_PA$HP == 'Richall',]
Combined_PAA$HP = 'All'

RichAT2.gg = ggplot(Combined_PAA, aes(x=spatial, y=S, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits=c(-0.5, 0.5)) + 
  theme_classic() + xlab('') + ylab('') +
  geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 0.6) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

Combined_PA2 = Combined_PA[!Combined_PA$HP == 'Richall',]
levels(Combined_PA2$HP)[4] = 'forbs'
RichT2.gg = ggplot(Combined_PA2, aes(x=spatial, y=S, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits=c(-0.5, 0.5)) + 
  theme_classic() + xlab('') + ylab('') +
  geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 0.6) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())


##

PAResid_FINlongA = PAResid_FINlong[PAResid_FINlong$HP == 'Richall',]
PAResid_FINlongA$HP = 'All'

MaprichAU2 <- ggplot(data = PAResid_FINlongA, aes(x=X, y=Y, color=S)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Richness residuals") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 


PAResid_FINlong2 = PAResid_FINlong[!PAResid_FINlong$HP == 'Richall',]
levels(PAResid_FINlong2$HP)[4] = 'forbs'

MaprichU2 <- ggplot(data = PAResid_FINlong2, aes(x=X, y=Y, color=S)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Richness residuals") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = PArichmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "right",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 

Combined_PA = rbind(PAResid_FINlong, PAResid_SITElong)
Combined_PAA = Combined_PA[Combined_PA$HP == 'Richall',]
Combined_PAA$HP = 'All'

RichAU2.gg = ggplot(Combined_PAA, aes(x=spatial, y=S, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits=c(0.9, 1.2)) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

Combined_PA2 = Combined_PA[!Combined_PA$HP == 'Richall',]
levels(Combined_PA2$HP)[4] = 'forbs'

RichU2.gg = ggplot(Combined_PA2, aes(x=spatial, y=S, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits=c(0.9, 1.2)) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

## plot all ....... Horiz



MeanPart = plot_grid(MaprichAM2, MaprichM2, 
                     RichAM2.gg, RichM2.gg, 
                     ncol = 2, nrow=2, align = "hv",labels=c('a', 'b', 'c', 'd'),
                     rel_widths = c(0.95, 4), rel_heights = c(1.5, 1))

TrendPart = plot_grid(MaprichAT2, MaprichT2, 
                      RichAT2.gg, RichT2.gg,
                      ncol = 2, nrow=2, align = "hv",labels=c('e', 'f','g', 'h'),
                      rel_widths = c(0.95, 4), rel_heights = c(1.5, 1))

VarPart = plot_grid(MaprichAU2, MaprichU2, 
                    RichAU2.gg, RichU2.gg, 
                    ncol = 2, nrow=2, align = "hv", labels=c('i', 'j','k', 'l'),
                    rel_widths = c(0.95, 4), rel_heights = c(1.5, 1))

grid.arrange(arrangeGrob(MeanPart, left = "1- Average"), 
             arrangeGrob(TrendPart, left = "2- Trend"), 
             arrangeGrob(VarPart, left = "3- Residuals variance"), 
             ncol = 1, heights = c(1,1,1), top='I- Predicted diversity: Richness')




##-...................................................................................................
##-...................................................................................................

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



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


ABmean_FINlongA = ABmean_FINlong[ABmean_FINlong$HP == 'eveall_mean',]
ABmean_FINlongA$HP = 'All'
MapABAM2 <- ggplot(data = ABmean_FINlongA, aes(x=X_mean, y=Y_mean, color=EVE)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Evenness mean") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30), limits=c(0,1)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = ABmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 


ABmean_FINlong2 = ABmean_FINlong[!ABmean_FINlong$HP == 'eveall_mean',]
levels(ABmean_FINlong2$HP) = c("herbaceous", "graminoids", "others", "forbs", "deciduous", "coniferous", "All")

MapABM2 <- ggplot(data = ABmean_FINlong2, aes(x=X_mean, y=Y_mean, color=EVE)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Evenness mean") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30), limits=c(0,1)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = ABmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 



Combined_AB = rbind(ABmean_FINlong, ABmean_SITElong)
Combined_ABA = Combined_AB[Combined_AB$HP == 'eveall_mean',]
Combined_ABA$HP = 'All'
ABAM2.gg = ggplot(Combined_ABA, aes(x=spatial, y=EVE, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits = c(0,1)) + #ylim(-100,100) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

Combined_AB2 = Combined_AB[!Combined_AB$HP == 'eveall_mean',]
levels(Combined_AB2$HP) = c("herbaceous", "graminoids", "others", "forbs", "deciduous", "coniferous", "All")

ABM2.gg = ggplot(Combined_AB2, aes(x=spatial, y=EVE, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits = c(0,1)) + #ylim(-100,100) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())




ABtrend_FIN = coef_trend_FINb[, c(1:7, 12:14)]
colnames(ABtrend_FIN)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
ABtrend_SITE = coef_trend_SITEb[, c(1:7, 9:10, 8)]
colnames(ABtrend_SITE)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')

ABtrend_FINlong = gather(ABtrend_FIN, HP, EVE, herbaceous:eveall, factor_key=TRUE)

colnames(ABtrend_SITE)[7] = 'eveall'
ABtrend_SITElong = gather(ABtrend_SITE, HP, EVE, herbaceous:eveall, factor_key=TRUE)

#ABtrend_FINlong$Slog = log(ABtrend_FINlong$S)
#ABtrend_SITElong$Slog = log(ABtrend_SITElong$S)
colnames(ABtrend_FINlong) = colnames(ABtrend_SITElong)

Combined_AB = rbind(ABtrend_FINlong, ABtrend_SITElong)

ABtrend_FINlongA = ABtrend_FINlong[ABtrend_FINlong$HP == 'eveall',]
ABtrend_FINlongA$HP = 'All'

MapABAT2 <- ggplot(data = ABtrend_FINlongA, aes(x=X, y=Y, color=EVE)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Evenness 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.035, 0.035)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank())+
  geom_point(data = ABmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 


ABtrend_FINlong2 = ABtrend_FINlong[!ABtrend_FINlong$HP == 'eveall',]

MapABT2 <- ggplot(data = ABtrend_FINlong2, aes(x=X, y=Y, color=EVE)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Evenness 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.035, 0.035)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = ABmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 



Combined_AB = rbind(ABtrend_FINlong, ABtrend_SITElong)
Combined_ABA = Combined_AB[Combined_AB$HP == 'eveall',]
Combined_ABA$HP = 'All'

ABAT2.gg = ggplot(Combined_ABA, aes(x=spatial, y=EVE, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits=c(-0.035, 0.035)) +# ylim(-0.5, 0.5) + 
  theme_classic() + xlab('') + ylab('') +
  geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 0.6) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

Combined_AB2 = Combined_AB[!Combined_AB$HP == 'eveall',]

ABT2.gg = ggplot(Combined_AB2, aes(x=spatial, y=EVE, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits = c(-0.035, 0.035)) + 
  theme_classic() + xlab('') + ylab('') +
  geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 0.6) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())


##

ABResid_FIN = Resid2_site_FINb[, c(1:7, 12:14)]

Resid2_site_SITEb$X = Moth_knn2$XEUREF
Resid2_site_SITEb$Y = Moth_knn2$YEUREF
ABResid_SITE = Resid2_site_SITEb[, c(1:7, 9:10, 8)]
colnames(ABResid_FIN)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')
colnames(ABResid_SITE)[1:6] = c('herbaceous', 'graminoids', 'others', 'forbs', 'deciduous', 'coniferous')


ABResid_FINlong = gather(ABResid_FIN, HP, EVE, herbaceous:eveall, factor_key=TRUE)
colnames(ABResid_SITE)[7] = 'eveall'
ABResid_SITElong = gather(ABResid_SITE, HP, EVE, herbaceous:eveall, factor_key=TRUE)

#ABResid_FINlong$Slog = log(ABResid_FINlong$S)
#ABResid_SITElong$Slog = log(ABResid_SITElong$S)
colnames(ABResid_FINlong) = colnames(ABResid_SITElong)
Combined_PA = rbind(ABResid_FINlong, ABResid_SITElong)

### Now plot everything

ABResid_FINlongA = ABResid_FINlong[ABResid_FINlong$HP == 'eveall',]
ABResid_FINlongA$HP = 'All'
MapABAU2 <- ggplot(data = ABResid_FINlongA, aes(x=X, y=Y, color=EVE)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Evenness residuals") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = ABmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "bottom",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 


ABResid_FINlong2 = ABResid_FINlong[!ABResid_FINlong$HP == 'eveall',]

MapABU2 <- ggplot(data = ABResid_FINlong2, aes(x=X, y=Y, color=EVE)) + geom_point(size=0.6) + 
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  #ggtitle("Evenness residuals") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_color_gradientn(colours = cividis(30)) + 
  theme_classic()  +
  theme(legend.position = "left",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=8), strip.background = element_blank()) +
  geom_point(data = ABmean_SITE, aes(x=X_mean, y=Y_mean), col='firebrick', size=0.5)+
  guides(color = guide_colorbar(title = "",
                                #label.position = "right",
                                #title.position = "left", title.vjust = 1,
                                # draw border around the legend
                                frame.colour = "black",
                                barwidth = 0.5,
                                barheight = 7)) 

Combined_AB = rbind(ABResid_FINlong, ABResid_SITElong)
Combined_ABA = Combined_AB[Combined_AB$HP == 'eveall',]
Combined_ABA$HP = 'All'

ABAU2.gg = ggplot(Combined_ABA, aes(x=spatial, y=EVE, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits = c(0.9, 1.2)) + #ylim(-100,100) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())

Combined_AB2 = Combined_AB[!Combined_AB$HP == 'eveall',]

ABU2.gg = ggplot(Combined_AB2, aes(x=spatial, y=EVE, fill= spatial)) +
  facet_wrap(~ HP, nrow=1, scales = "free", strip.position = "top") +
  geom_violin(trim=FALSE, draw_quantiles = c(0.25, 0.5, 0.75)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 3), limits = c(0.9, 1.2)) + #ylim(-100,100) + 
  theme_classic() + xlab('') + ylab('') +
  #geom_hline(yintercept = 0, linetype="dashed", colour='gray60', size= 1) +
  coord_flip() + scale_fill_manual(values=c("beige", "firebrick")) + 
  theme(legend.position = "none", axis.text=element_text(size=7), 
        strip.text = element_text(size=8), strip.background = element_blank())


## plot all ....... Horiz


MeanPart = plot_grid(MapABAM2, MapABM2, 
                     ABAM2.gg, ABM2.gg, 
                    ncol = 2, nrow=2, align = "hv", labels=c('a', 'b', 'c', 'd'),
                    rel_widths = c(0.95, 4), rel_heights = c(1.5, 1))

TrendPart = plot_grid(MapABAT2, MapABT2, 
                      ABAT2.gg, ABT2.gg, 
                     ncol = 2, nrow=2, align = "hv", labels=c('e', 'f','g', 'h'),
                     rel_widths = c(0.95, 4), rel_heights = c(1.5, 1))

VarPart = plot_grid(MapABAU2, MapABU2, 
                    ABAU2.gg, ABU2.gg, 
                     ncol = 2, nrow=2, align = "hv", labels=c('i', 'j','k', 'l'),
                     rel_widths = c(0.95, 4), rel_heights = c(1.5, 1))


grid.arrange(arrangeGrob(MeanPart, left = "1- Average"), 
             arrangeGrob(TrendPart, left = "2- Trend"), 
             arrangeGrob(VarPart, left = "3- Residuals variance"), 
             ncol = 1, heights = c(1,1,1), top='II- Predicted diversity: Evenness')



