
library(ggplot2)
library(ggridges)
library(grid)
library(data.table)
library(reshape2)

library(Hmsc)
library(reshape2)
library(abind)
library(tidyr)
library(dplyr)
#library(berryFunctions)
#library(lemon)
library(cowplot)
library(patchwork)
library(ggh4x)
library(ggpubr)

library(vioplot)




# ------------------------------------------------------------

## Simplified plot ------------------------------------------------------------
# ------------------------------------------------------------

setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/3.ModelRun")

load('VPall_traits2FIN.RDATA') ###  FINLAND
##                                    OR
load('VPall_traits_Sitesprep.RDATA') ### SITES


# -----------------------------------------------------------------
# Table of species and traits
data.table.traits = unique(VP_PA.all[, c(1,5,6)])

levels(data.table.traits$HostPlantForm)[1] = 'Herbaceous'
levels(data.table.traits$HostPlantForm)[4] = 'Dwarf shrubs'
levels(data.table.traits$HostPlantForm)[5] = 'Deciduous'


data.table.traits$Wingspan = NA
data.table.traits$Wingspan[data.table.traits$WingspanF<10] ='w = 0'
data.table.traits$Wingspan[data.table.traits$WingspanF>=17.9 & data.table.traits$WingspanF<24.1] ='17.9<w<24.1'
data.table.traits$Wingspan[data.table.traits$WingspanF>=24.1 & data.table.traits$WingspanF<30.3] ='24.1<w<30.3'
data.table.traits$Wingspan[data.table.traits$WingspanF>=30.3 & data.table.traits$WingspanF<36.5] ='30.3<w<36.5'
data.table.traits$Wingspan[data.table.traits$WingspanF>=36.5 & data.table.traits$WingspanF<43] ='36.5<w<43'
data.table.traits$Wingspan[data.table.traits$WingspanF>=48] ='w=55.5'
table(data.table.traits$Wingspan)

data.table.traits$Wingspan = factor(data.table.traits$Wingspan, 
                                    levels = c('w = 0', '17.9<w<24.1', '24.1<w<30.3', '30.3<w<36.5',
                                               '36.5<w<43', 'w=55.5'))

data.table.traits = data.table.traits[order(data.table.traits$HostPlantForm, data.table.traits$WingspanF),]
library(dplyr)

remotes::install_github("ewenme/ghibli")
library(ghibli)
ghibli_palettes -> ghibli_palettes_list
ghibli_palettes$MononokeMedium

c("#838A90FF", '#9FA7BEFF', '#B3B8B1FF',
  '#E7A79BFF', '#F2C695FF', '#F5EDC9FF')

data.table.traits$color = NA
data.table.traits= data.table.traits %>%
  group_by(Wingspan) %>%
  mutate(color = case_when( Wingspan %in% c('w = 0') ~ "#838A90FF", 
                            Wingspan %in% c('17.9<w<24.1') ~ '#9FA7BEFF',
                            Wingspan %in% c('24.1<w<30.3') ~ '#B3B8B1FF',
                            Wingspan %in% c('30.3<w<36.5') ~ '#E7A79BFF', 
                            Wingspan %in% c('36.5<w<43') ~ '#F2C695FF',
                            Wingspan %in% c('w=55.5') ~ '#F5EDC9FF',
  )) 


data.table.traits = arrange(data.table.traits, desc(HostPlantForm), desc(WingspanF))

# library(DT)
# datatable(data.table.traits) %>% formatStyle(
#   'HostPlantForm',
#   backgroundColor = styleEqual(c('Graminoids', 'Herbaceous',  
#                                  'Others', 'Forbs',
#                                  'Coniferous', 'Deciduous'), 
#                                c("cyan", "salmon4", 'purple', 'pink',
#                                  "forestgreen", "orange"))
# )


setwd("D:/Helsinki/Manuscript_moth/Figures")
library(gt)

data.table.traits_s1 = as.data.frame(data.table.traits[1:39,])
data.table.traits_s2 = as.data.frame(data.table.traits[40:78,])

gt1 = data.table.traits_s1 |>
  #dplyr::select(-contains("WingspanF")) |>
  gt() |>
  tab_style(
    style = cell_fill(color = from_column(column = "color")),
    locations = cells_body(columns = WingspanF)
  ) |>
  # tab_style(
  #   style = cell_fill(color = from_column(column = "colorHP")),
  #   locations = cells_body(columns = HostPlantForm)
  # ) |>
  # data_color(
  #   columns = c(WingspanF), method = "factor",
  #   palette = c("#838A90FF", '#BA968AFF', 
  #               '#9FA7BEFF', '#B3B8B1FF',
  #               '#E7A79BFF', '#F2C695FF',
  #               '#F5EDC9FF')
  # )|>
  data_color(
    columns = c(HostPlantForm), method = "factor",
    palette = c("cyan", "salmon4", 'purple', 'pink',
                "forestgreen", "orange")
  ) |>
  tab_options(table.width = pct(80),
              table.font.size = "small",
              data_row.padding = px(1)
  ) 
gt1 |>  cols_hide(columns = c(color, Wingspan))


gt2 = data.table.traits_s2 |>
  #dplyr::select(-contains("WingspanF")) |>
  gt() |>
  tab_style(
    style = cell_fill(color = from_column(column = "color")),
    locations = cells_body(columns = WingspanF)
  ) |>
  # tab_style(
  #   style = cell_fill(color = from_column(column = "colorHP")),
  #   locations = cells_body(columns = HostPlantForm)
  # ) |>
  # data_color(
  #   columns = c(WingspanF), method = "factor",
  #   palette = c("#838A90FF", '#BA968AFF', 
  #               '#9FA7BEFF', '#B3B8B1FF',
  #               '#E7A79BFF', '#F2C695FF',
  #               '#F5EDC9FF')
  # )|>
  data_color(
    columns = c(HostPlantForm), method = "factor",
    palette = c("cyan", "salmon4", 'purple', 'pink',
                "forestgreen", "orange")
  ) |>
  tab_options(table.width = pct(80),
              table.font.size = "small",
              data_row.padding = px(1)
  ) 
gt2 |>  cols_hide(columns = c(color, Wingspan))



## ------------------------------------------------------------------------

VP.PA = VP_PA.all

class(VP.PA$HostPlantForm)
levels(VP.PA$HostPlantForm)[4] = 'Dwarf shrubs'
levels(VP.PA$HostPlantForm)[5] = 'Deciduous'

colnames(VP.PA)[6] = 'HP'
gpa1 = ggplot(VP.PA, aes(x = VPhab , y = VP, fill = HP)) +
  #geom_point()+
  stat_density_2d(aes(alpha = ..level.., fill = HP), geom = "polygon", 
                  position = "identity", h=0.2, n=50, bins=10,
                  show.legend = TRUE) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange"))+
  theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+
  ylim(0,0.86) + xlim(0,0.86) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Hostplant form"
  ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none') +
  #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
  theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86),
        legend.position = "top") # <- specify position of legend


# legend_b <- get_legend(
#   gpa1 +
#     guides(color = guide_legend(nrow = 1)) +
#     theme(legend.position = "bottom",
#           legend.text = element_text(size=13),
#           legend.title = element_text(size=16))
# )

legend_b = cowplot::get_plot_component(gpa1, 'guide-box-top', return_all = TRUE)



gpa1 = ggplot(VP.PA, aes(x = VPhab , y = VP, fill = HP)) +
  #geom_point()+
  stat_density_2d(aes(alpha = ..level.., fill = HP), geom = "polygon", 
                  position = "identity", h=0.2, n=50, bins=10,
                  show.legend = FALSE) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange"))+
  theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+
  ylim(0,0.86) + xlim(0,0.86) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Hostplant form"
  ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none') +
  #ggtitle('PA model') +
  #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
  theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))


xmargpa1 <- ggplot(VP.PA, aes(x = VPhab, y=..density.., fill = HP)) + 
  geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange")) + 
  theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
  theme(legend.position = "none")+ ggtitle('PA model') +
  theme(plot.title = element_text(hjust = 1.2, vjust = -10))


ymargpa1 <- ggplot(VP.PA) + 
  geom_density(alpha = 0.4, aes(x = VP, fill = HP)) + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange"))+ 
  theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
  theme(legend.position = "none") + 
  coord_flip()


### AB
VP.AB = VP_AB.all
colnames(VP.AB)[6] = 'HP'
levels(VP.AB$HP)[4] = 'Forbs'
levels(VP.AB$HP)[5] = 'Deciduous'

gab1 = ggplot(VP.AB, aes(x = VPhab , y = VP, fill = HP)) +
  #geom_point()+
  stat_density_2d(aes(alpha = ..level.., fill = HP), geom = "polygon", 
                  position = "identity", h=0.2, n=100, bins=10,
                  show.legend = FALSE) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange"))+
  theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+
  ylim(0,0.86) + xlim(0,0.86) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Hostplant form"
  ) + geom_abline(intercept = 0, slope = 1) + #guides(alpha='none') +
  #ggtitle('AB model') +
  #labs(subtitle=paste('  SST = ', round(SST.yAB.Heter,1),'\n SSE = ', round(SSE.yAB.Heter,1))) + 
  theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))



xmargab1 <- ggplot(VP.AB, aes(x = VPhab, y=..density.., fill = HP)) + 
  geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange")) + 
  theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
  theme(legend.position = "none")+ ggtitle('AB model') +
  theme(plot.title = element_text(hjust = 1.2, vjust = -10))


ymargab1 <- ggplot(VP.AB) + 
  geom_density(alpha = 0.4, aes(x = VP, fill = HP)) + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('Graminoids'="cyan", 'Herbaceous'="salmon4", 
                               'Others'='purple', 'Dwarf shrubs'='pink',
                               'Coniferous'="forestgreen", 'Deciduous'="orange"))+ 
  theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
  theme(legend.position = "none") + 
  coord_flip()


plot_grid(xmargpa1, NULL, xmargab1, NULL, 
          gpa1, ymargpa1, gab1, ymargab1, 
          ncol = 4, nrow=2, align = "hv", #labels = c('PA model', NULL, 'AB model',NULL),
          rel_widths = c(2, 1), rel_heights = c(1, 2))


# 
# 
# ### ..........     wingF  3 groups
# VP.PA = VP_PA.all[VP_PA.all$WingspanF>0,]
# summary(VP.PA$WingspanF)
# VP.PA$wing2 = NA
# VP.PA$wing2[VP.PA$WingspanF<28] ='w<28'
# #VP.PA$wing2[VP.PA$WingspanF>=25 & VP.PA$WingspanF<35] ='25>w>35'
# VP.PA$wing2[VP.PA$WingspanF>=28 & VP.PA$WingspanF<36] ='28<w<34'
# VP.PA$wing2[VP.PA$WingspanF>=34] ='w>34'
# table(VP.PA$wing2)
# 
# gpa2 = ggplot(VP.PA, aes(x = VPhab , y = VP)) +
#   #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
#   stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
#                   position = "identity", h=0.2, n=50, bins=20,
#                   show.legend = TRUE) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
#   theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
#   ylim(0,0.86) + xlim(0,0.86) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Wing size"
#   ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none') +
#   #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
#   theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))
# 
# 
# 
# legend_c <- get_legend(
#   gpa2 + 
#     guides(color = guide_legend(nrow = 1)) +
#     theme(legend.position = "bottom", 
#           legend.text = element_text(size=13), 
#           legend.title = element_text(size=16))
# )
# 
# 
# 
# gpa2 = ggplot(VP.PA, aes(x = VPhab , y = VP)) +
#   #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
#   stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
#                   position = "identity", h=0.2, n=50, bins=20,
#                   show.legend = TRUE) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
#   theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
#   ylim(0,0.86) + xlim(0,0.86) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Wing size"
#   ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none', fill = 'none') +
#   #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
#   theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))
# 
# 
# 
# xmargpa2 <- ggplot(VP.PA, aes(x = VPhab, y=..density.., fill = wing2)) + 
#   geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
#   theme(legend.position = "none")+ 
#   theme(plot.title = element_text(hjust = 1.2, vjust = -10))
# 
# 
# ymargpa2 <- ggplot(VP.PA) + 
#   geom_density(alpha = 0.4, aes(x = VP, fill = wing2)) + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
#   theme(legend.position = "none") + 
#   coord_flip()
# 
# 
# 
# ### AB
# VP.AB = VP_AB.all[VP_AB.all$WingspanF>0,]
# VP.AB$wing2 = NA
# 
# VP.AB$wing2[VP.AB$WingspanF<28] ='w<28'
# #VP.AB$wing2[VP.AB$WingspanF>=25 & VP.AB$WingspanF<35] ='25>w>35'
# VP.AB$wing2[VP.AB$WingspanF>=28 & VP.AB$WingspanF<34] ='28<w<34'
# VP.AB$wing2[VP.AB$WingspanF>=34] ='w>34'
# 
# gab2 = ggplot(VP.AB, aes(x = VPhab,  y = VP)) +
#   #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
#   stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
#                   position = "identity", h=0.2, n=50, bins=20,
#                   show.legend = TRUE) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
#   theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
#   ylim(0,0.86) + xlim(0,0.86) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Wing size"
#   ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none', fill = 'none') +
#   #labs(subtitle=paste('  SST = ', round(SST.yAB.Heter,1),'\n SSE = ', round(SSE.yAB.Heter,1))) + 
#   theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))
# 
# xmargab2 <- ggplot(VP.AB, aes(x = VPhab, y=..density.., fill = wing2)) + 
#   geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
#   theme(legend.position = "none")+
#   theme(plot.title = element_text(hjust = 1.2, vjust = -10))
# 
# 
# ymargab2 <- ggplot(VP.AB) + 
#   geom_density(alpha = 0.4, aes(x = VP, fill = wing2)) + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<28'="yellow", #'25>w>35'="#fcd368", 
#                                '28<w<34'='#fea043', 'w>34'='red'))+
#   theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
#   theme(legend.position = "none") + 
#   coord_flip()
# 
# 
# 
# plot_grid(#xmargpa2, NULL, xmargab2, NULL, 
#                     gpa2, #ymargpa2, 
#                     gab2, #ymargab2, 
#                     ncol = 2, nrow=1, align = "hv", labels = c('C - ', NULL, 'D - ',NULL),
#                     rel_widths = c(2, 2), rel_heights = c(2, 2))
# 

## 7 groups to match host plant


### ..........     wingF  7 groups - regular
VP.PA = VP_PA.all#[VP_PA.all$WingspanF>0,]
summary(VP.PA$WingspanF)
hist(VP.PA$WingspanF)


VP.PA$wing2 = NA
VP.PA$wing2[VP.PA$WingspanF<10] ='w = 0'
VP.PA$wing2[VP.PA$WingspanF>=17.9 & VP.PA$WingspanF<24.1] ='17.9<w<24.1'
VP.PA$wing2[VP.PA$WingspanF>=24.1 & VP.PA$WingspanF<30.3] ='24.1<w<30.3'
VP.PA$wing2[VP.PA$WingspanF>=30.3 & VP.PA$WingspanF<36.5] ='30.3<w<36.5'
VP.PA$wing2[VP.PA$WingspanF>=36.5 & VP.PA$WingspanF<43] ='36.5<w<43'
VP.PA$wing2[VP.PA$WingspanF>=48] ='w=55.5'
table(VP.PA$wing2)

VP.PA$wing2 = factor(VP.PA$wing2, levels = c('w = 0', '17.9<w<24.1', '24.1<w<30.3', '30.3<w<36.5',
                                               '36.5<w<43', 'w=55.5'))

remotes::install_github("ewenme/ghibli")
library(ghibli)
ghibli_palettes -> ghibli_palettes_list
ghibli_palettes$MononokeMedium
c("#838A90FF", '#9FA7BEFF', '#B3B8B1FF',
  '#E7A79BFF', '#F2C695FF', '#F5EDC9FF')

gpa2leg = ggplot(VP.PA, aes(x = VPhab , y = VP)) +
  #geom_point(aes(size = wing2F, colour = wing2F))+ scale_size(range = c(1, 2))+
  stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon",
                  position = "identity", h=0.2, n=50, bins=20,
                  show.legend = TRUE) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
  theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
  ylim(0,0.86) + xlim(0,0.86) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Wing size"
  ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none') +
  #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) +
  theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86),
        legend.position = "top") # <- specify position of legend

legend = cowplot::get_plot_component(gpa2leg, 'guide-box-top', return_all = TRUE)
cowplot::ggdraw(legend)


gpa2leg2 = gpa2leg +
  guides(color = guide_legend(nrow = 1)) +
  theme(legend.position = "bottom",
        legend.text = element_text(size=13),
        legend.title = element_text(size=16))

# legend_c <- cowplot::get_plot_component(gpa2leg2, 'guide-box-top', return_all = TRUE)
# cowplot::ggdraw(legend_c)
# 
# 
# legend_c <- get_legend(
#   gpa2leg + 
#     guides(color = guide_legend(nrow = 1)) +
#     theme(legend.position = "bottom", 
#           legend.text = element_text(size=13), 
#           legend.title = element_text(size=16))
# )


gpa2 = ggplot(VP.PA, aes(x = VPhab , y = VP)) +
  #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
  stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
                  position = "identity", h=0.2, n=50, bins=20,
                  show.legend = TRUE) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
  theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
  ylim(0,0.86) + xlim(0,0.86) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Wing size"
  ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none', fill = 'none') +
  #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
  theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))



xmargpa2 <- ggplot(VP.PA, aes(x = VPhab, y=..density.., fill = wing2)) + 
  geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
  theme(legend.position = "none")+ 
  theme(plot.title = element_text(hjust = 1.2, vjust = -10))


ymargpa2 <- ggplot(VP.PA) + 
  geom_density(alpha = 0.4, aes(x = VP, fill = wing2)) + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
  theme(legend.position = "none") + 
  coord_flip()



### AB
VP.AB = VP_AB.all#[VP_AB.all$WingspanF>0,]
summary(VP.AB$WingspanF)
hist(VP.AB$WingspanF)


VP.AB$wing2 = NA
VP.AB$wing2[VP.AB$WingspanF<10] ='w = 0'
VP.AB$wing2[VP.AB$WingspanF>=17.9 & VP.AB$WingspanF<24.1] ='17.9<w<24.1'
VP.AB$wing2[VP.AB$WingspanF>=24.1 & VP.AB$WingspanF<30.3] ='24.1<w<30.3'
VP.AB$wing2[VP.AB$WingspanF>=30.3 & VP.AB$WingspanF<36.5] ='30.3<w<36.5'
VP.AB$wing2[VP.AB$WingspanF>=36.5 & VP.AB$WingspanF<43] ='36.5<w<43'
VP.AB$wing2[VP.AB$WingspanF>=48] ='w=55.5'
table(VP.AB$wing2)

VP.AB$wing2 = factor(VP.AB$wing2, levels = c('w = 0', '17.9<w<24.1', '24.1<w<30.3', '30.3<w<36.5',
                                             '36.5<w<43', 'w=55.5'))


gab2 = ggplot(VP.AB, aes(x = VPhab,  y = VP)) +
  #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
  stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
                  position = "identity", h=0.2, n=50, bins=20,
                  show.legend = TRUE) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
  theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
  ylim(0,0.86) + xlim(0,0.86) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Wing size"
  ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none', fill = 'none') +
  #labs(subtitle=paste('  SST = ', round(SST.yAB.Heter,1),'\n SSE = ', round(SSE.yAB.Heter,1))) + 
  theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))

xmargab2 <- ggplot(VP.AB, aes(x = VPhab, y=..density.., fill = wing2)) + 
  geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
  theme(legend.position = "none")+
  theme(plot.title = element_text(hjust = 1.2, vjust = -10))


ymargab2 <- ggplot(VP.AB) + 
  geom_density(alpha = 0.4, aes(x = VP, fill = wing2)) + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('w = 0'="#838A90FF", 
                               '17.9<w<24.1'='#9FA7BEFF', '24.1<w<30.3'= '#B3B8B1FF',
                               '30.3<w<36.5' = '#E7A79BFF', '36.5<w<43'='#F2C695FF',
                               'w=55.5'='#F5EDC9FF'))+
  theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
  theme(legend.position = "none") + 
  coord_flip()



plot_grid(xmargpa2, NULL, xmargab2, NULL, 
          gpa2, ymargpa2, 
          gab2, ymargab2, 
          ncol = 4, nrow=2, align = "hv", labels = c('C - ', NULL, 'D - ',NULL),
          rel_widths = c(2, 1), rel_heights = c(1, 2))


# 
# ### ..........     wingF  7 groups - based on equal number of species
# VP.PA = VP_PA.all[VP_PA.all$WingspanF>0,]
# summary(VP.PA$WingspanF)
# hist(VP.PA$WingspanF)
# 
# # create 7 groups
# WG = unique(VP.PA$WingspanF)[order(unique(VP.PA$WingspanF))]
# dplyr::ntile(WG, 7)
# 
# cbind(WG, dplyr::ntile(WG, 7))
# 
# VP.PA$wing2 = NA
# VP.PA$wing2[VP.PA$WingspanF<22] ='w<22'
# VP.PA$wing2[VP.PA$WingspanF>=22 & VP.PA$WingspanF<26.3] ='22<w<26.3'
# VP.PA$wing2[VP.PA$WingspanF>=26.3 & VP.PA$WingspanF<30.2] ='26.3<w<30.2'
# VP.PA$wing2[VP.PA$WingspanF>=30.2 & VP.PA$WingspanF<32.6] ='30.2<w<32.6'
# VP.PA$wing2[VP.PA$WingspanF>=32.6 & VP.PA$WingspanF<34.5] ='32.6<w<34.5'
# VP.PA$wing2[VP.PA$WingspanF>=34.5 & VP.PA$WingspanF<37.1] ='34.5<w<37.1'
# VP.PA$wing2[VP.PA$WingspanF>=37.1] ='w>37.1'
# table(VP.PA$wing2)
# 
# #remotes::install_github("ewenme/ghibli")
# library(ghibli)
# ghibli_palettes -> ghibli_palettes_list
# ghibli_palettes$MononokeMedium
# 
# gpa2leg = ggplot(VP.PA, aes(x = VPhab , y = VP)) +
#   #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
#   stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
#                   position = "identity", h=0.2, n=50, bins=20,
#                   show.legend = TRUE) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
#   theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
#   ylim(0,0.86) + xlim(0,0.86) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Wing size"
#   ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none') +
#   #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
#   theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86), 
#         legend.position = "top") # <- specify position of legend
# 
# legend = cowplot::get_plot_component(plot, 'guide-box-top', return_all = TRUE)
# cowplot::ggdraw(legend)
# 
# 
# gpa2leg2 = gpa2leg + 
#   guides(color = guide_legend(nrow = 1)) +
#   theme(legend.position = "bottom", 
#         legend.text = element_text(size=13), 
#         legend.title = element_text(size=16))
# 
# legend_c <- cowplot::get_plot_component(gpa2leg2, 'guide-box-top', return_all = TRUE)
# cowplot::ggdraw(legend_c)
# 
# get_legend2 <- function(plot, legend = NULL) {
#   if (is.ggplot(plot)) {
#     gt <- ggplotGrob(plot)
#   } else {
#     if (is.grob(plot)) {
#       gt <- plot
#     } else {
#       stop("Plot object is neither a ggplot nor a grob.")
#     }
#   }
#   pattern <- "guide-box"
#   if (!is.null(legend)) {
#     pattern <- paste0(pattern, "-", legend)
#   }
#   indices <- grep(pattern, gt$layout$name)
#   not_empty <- !vapply(
#     gt$grobs[indices], 
#     inherits, what = "zeroGrob", 
#     FUN.VALUE = logical(1)
#   )
#   indices <- indices[not_empty]
#   if (length(indices) > 0) {
#     return(gt$grobs[[indices[1]]])
#   }
#   return(NULL)
# }
# 
# 
# legend_c <- get_legend2(
#   gpa2leg + 
#     guides(color = guide_legend(nrow = 1)) +
#     theme(legend.position = "bottom", 
#           legend.text = element_text(size=13), 
#           legend.title = element_text(size=16))
# )
# 
# 
# 
# gpa2 = ggplot(VP.PA, aes(x = VPhab , y = VP)) +
#   #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
#   stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
#                   position = "identity", h=0.2, n=50, bins=20,
#                   show.legend = TRUE) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
#   theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
#   ylim(0,0.86) + xlim(0,0.86) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Wing size"
#   ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none', fill = 'none') +
#   #labs(subtitle=paste('  SST = ', round(SST.yPA.Wat,1),'\n SSE = ', round(SSE.yPA.Wat,1))) + 
#   theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))
# 
# 
# 
# xmargpa2 <- ggplot(VP.PA, aes(x = VPhab, y=..density.., fill = wing2)) + 
#   geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
#   theme(legend.position = "none")+ 
#   theme(plot.title = element_text(hjust = 1.2, vjust = -10))
# 
# 
# ymargpa2 <- ggplot(VP.PA) + 
#   geom_density(alpha = 0.4, aes(x = VP, fill = wing2)) + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
#   theme(legend.position = "none") + 
#   coord_flip()
# 
# 
# 
# ### AB
# VP.AB = VP_AB.all[VP_AB.all$WingspanF>0,]
# VP.AB$wing2 = NA
# VP.AB$wing2[VP.AB$WingspanF<22] ='w<22'
# VP.AB$wing2[VP.AB$WingspanF>=22 & VP.AB$WingspanF<26.3] ='22<w<26.3'
# VP.AB$wing2[VP.AB$WingspanF>=26.3 & VP.AB$WingspanF<30.2] ='26.3<w<30.2'
# VP.AB$wing2[VP.AB$WingspanF>=30.2 & VP.AB$WingspanF<32.6] ='30.2<w<32.6'
# VP.AB$wing2[VP.AB$WingspanF>=32.6 & VP.AB$WingspanF<34.5] ='32.6<w<34.5'
# VP.AB$wing2[VP.AB$WingspanF>=34.5 & VP.AB$WingspanF<37.1] ='34.5<w<37.1'
# VP.AB$wing2[VP.AB$WingspanF>=37.1] ='w>37.1'
# table(VP.AB$wing2)
# 
# gab2 = ggplot(VP.AB, aes(x = VPhab,  y = VP)) +
#   #geom_point(aes(size = WingspanF, colour = WingspanF))+ scale_size(range = c(1, 2))+
#   stat_density_2d(aes(alpha = ..level.., fill = wing2), geom = "polygon", 
#                   position = "identity", h=0.2, n=50, bins=20,
#                   show.legend = TRUE) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   #scale_colour_gradient(low = "yellow", high = "red", na.value = 0) +
#   theme_pubr(legend = c(1,0.9)) + theme(plot.margin = unit(c(0,0,0,0), "cm"))+
#   ylim(0,0.86) + xlim(0,0.86) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Wing size"
#   ) + geom_abline(intercept = 0, slope = 1) + guides(alpha='none', fill = 'none') +
#   #labs(subtitle=paste('  SST = ', round(SST.yAB.Heter,1),'\n SSE = ', round(SSE.yAB.Heter,1))) + 
#   theme(plot.subtitle = element_text(size = 10, face = "bold", hjust = 0.86))
# 
# xmargab2 <- ggplot(VP.AB, aes(x = VPhab, y=..density.., fill = wing2)) + 
#   geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
#   theme(legend.position = "none")+
#   theme(plot.title = element_text(hjust = 1.2, vjust = -10))
# 
# 
# ymargab2 <- ggplot(VP.AB) + 
#   geom_density(alpha = 0.4, aes(x = VP, fill = wing2)) + xlim(0.05, 0.78) +
#   scale_fill_manual(values = c('w<22'="#838A90FF", '22<w<26.3' = '#BA968AFF', 
#                                '26.3<w<30.2'='#9FA7BEFF', '30.2<w<32.6'= '#B3B8B1FF',
#                                '32.6<w<34.5' = '#E7A79BFF', '34.5<w<37.1'='#F2C695FF',
#                                'w>37.1'='#F5EDC9FF'))+
#   theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
#   theme(legend.position = "none") + 
#   coord_flip()
# 
# 
# 
# plot_grid(xmargpa2, NULL, xmargab2, NULL, 
#           gpa2, ymargpa2, 
#           gab2, ymargab2, 
#           ncol = 4, nrow=2, align = "hv", labels = c('C - ', NULL, 'D - ',NULL),
#           rel_widths = c(2, 1), rel_heights = c(1, 2))
# 
# 


# Arranging the plot using cowplot
TopPart = plot_grid(xmargpa1, NULL, xmargab1, NULL, 
                    gpa1, ymargpa1, gab1, ymargab1, 
                    ncol = 4, nrow=2, align = "hv", labels = c('a) ', ' ', ' ', 'b) '),
                    rel_widths = c(2, 1), rel_heights = c(1, 2), label_size = 16)

BotPart = plot_grid(xmargpa2, NULL, xmargab2, NULL, 
                    gpa2, ymargpa2, gab2, ymargab2, 
                    ncol = 4, nrow=2, align = "hv", labels = c('c) ', ' ', ' ', 'd) '),
                    rel_widths = c(2, 1), rel_heights = c(1, 2), label_size = 16)


plot_grid(TopPart, legend_b, BotPart, legend, 
          ncol = 1, rel_heights = c(0.7,.1,0.7, 0.1))

plot_grid(BotPart, legend, 
          ncol = 1, rel_heights = c(0.7, 0.1))
