LPJ-GUESS-OutputPlots / Biomass_Ensemble.R
Biomass_Ensemble.R
Raw
library(tidyverse)
library(lubridate)
library(patchwork)
library(RColorBrewer) #color schemes
library(rio)
library(rgdal)
library(fields)
library(maptools)
library(sf)
library(sp)
library(raster)
library(maps)
library(metR)
library(ggmap)
library(data.table)
library(geofacet)
library(matrixStats)

borders <- sf::st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))

setwd("C:/Users/CMCC-ESM2/")
gall16 <- read.table("cmass.out", sep="", header=T)
gall16$Year <- as.POSIXct(paste(gall16$Year), format = "%Y")
gall16 <- as.data.frame(gall16)
All_yr1 <- gall16 %>% group_by(year=floor_date(Year, "year"))
Per_yr1 <- All_yr1 %>% mutate(year=year(Year)) %>%
  group_by(Year, year)
Per_yr1 <- Per_yr1 %>% filter(year < 2100)


CM <- Per_yr1 %>% filter(year <= "1990")
# CM <- Per_yr1 %>% filter(year >= "2071")
# CM <- Per_yr1 %>% filter(year >= "2031" & year <= "2060")
# CM <- Per_yr1 %>% filter(year >= "1991" & year <= "2020")


Mean_Ab <- CM %>% group_by(Lon, Lat) %>% summarise(Ab_mean = mean(Abi_nor)) 

Mean_Ced <- CM %>% group_by(Lon, Lat) %>% summarise(Ced_mean = mean(Ced_lib)) 

Mean_BES <- CM %>% group_by(Lon, Lat) %>% summarise(BES_mean = mean(BES)) 

Mean_BP <-  CM %>% group_by(Lon, Lat) %>% summarise(BP_mean = mean(Bet_pen))

Mean_BPU <- CM %>% group_by(Lon, Lat) %>% summarise(BPU_mean = mean(Bet_pub))

Mean_Car <- CM %>% group_by(Lon, Lat) %>% summarise(Car_mean = mean(Car_bet))

Mean_Cor <- CM %>% group_by(Lon, Lat) %>% summarise(Cor_mean = mean(Cor_ave))

Mean_Fag <- CM %>% group_by(Lon, Lat) %>% summarise(Fag_mean = mean(Fag_syl))

Mean_Fra <- CM %>% group_by(Lon, Lat) %>% summarise(Fra_mean = mean(Fra_exc))

Mean_Jun <- CM %>% group_by(Lon, Lat) %>% summarise(Jun_mean = mean(Jun_oxy))

Mean_MRS <- CM %>% group_by(Lon, Lat) %>% summarise(MRS_mean = mean(MRS))

Mean_PA <-  CM %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Pic_abi))

Mean_AG <-  CM %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Aln_glu))

Mean_PS <-  CM %>%  group_by(Lon, Lat) %>% summarise(PS_mean = mean(Pin_syl))

Mean_PN <-  CM %>%  group_by(Lon, Lat) %>% summarise(PN_mean = mean(Pin_nig))

Mean_PB <-  CM %>%  group_by(Lon, Lat) %>% summarise(PB_mean = mean(Pin_bru))

Mean_PH <-  CM %>%  group_by(Lon, Lat) %>% summarise(PH_mean = mean(Pin_hal))

Mean_Pop <- CM %>% group_by(Lon, Lat) %>% summarise(Pop_mean = mean(Pop_tre))

Mean_QC <-  CM %>% group_by(Lon, Lat) %>% summarise(QC_mean = mean(Que_coc))

Mean_QI <-  CM %>% group_by(Lon, Lat) %>% summarise(QI_mean = mean(Que_ile))

Mean_QP <-  CM %>% group_by(Lon, Lat) %>% summarise(QP_mean = mean(Que_pub))

Mean_QR <-  CM %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_rob))

Mean_QM <-  CM %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_mac))

Mean_Til <- CM %>% group_by(Lon, Lat) %>% summarise(Til_mean = mean(Til_cor))

Mean_Ulm <- CM %>% group_by(Lon, Lat) %>% summarise(Ulm_mean = mean(Ulm_gla))

Mean_C3 <-  CM %>% group_by(Lon, Lat) %>% summarise(C3_mean = mean(C3_gr))


CM_Species <- cbind.data.frame(Mean_Ab[3], Mean_Ced[3], Mean_BES[3], Mean_BP[3], Mean_BPU[3], Mean_Car[3], Mean_Cor[3], Mean_Fag[3], 
                            Mean_Fra[3], Mean_Jun[3], Mean_MRS[3], Mean_PA[3], Mean_AG[3], Mean_PS[3], Mean_PN[3], Mean_PB[3], 
                            Mean_PH[3], Mean_Pop[3], Mean_QC[3], Mean_QI[3], Mean_QP[3], Mean_QR[3], Mean_QM[3], Mean_Til[3], 
                            Mean_Ulm[3], Mean_C3[3])
colnames(CM_Species) <- c("Abies spp.", "C. libani", "Boreal shrub", "B. pendula", "B. pubescens", "Carpinus spp.", "C. avellana",
                       "Fagus spp.","F. excelsior", "J. oxycedrus", "Med. shrub", "Picea spp.", "A. glutinosa", "P. sylvestris",
                       "P. nigra", "P. brutia", "P. halepensis", "P. tremula", "Q. coccifera", "Q. ilex", "Q. pubescens", 
                       "Q. robur", "Q. macranthera", "T. cordata",  "U. glabra", "Grass")



setwd("C:/Users/NorESM/")
gall16 <- read.table("cmass.out", sep="", header=T)
gall16$Year <- as.POSIXct(paste(gall16$Year), format = "%Y")
gall16 <- as.data.frame(gall16)
All_yr2 <- gall16 %>% group_by(year=floor_date(Year, "year"))
Per_yr2 <- All_yr2 %>% mutate(year=year(Year)) %>%
  group_by(Year, year)
Per_yr2 <- Per_yr2 %>% filter(year < 2100)

Nor <- Per_yr2 %>% filter(year <= "1990")
# Nor <- Per_yr2 %>% filter(year >= "2071")
# Nor <- Per_yr2 %>% filter(year >= "2031" & year <= "2060")
# Nor <- Per_yr2 %>% filter(year >= "1991" & year <= "2020")

Mean_Ab <- Nor %>% group_by(Lon, Lat) %>% summarise(Ab_mean = mean(Abi_nor)) 

Mean_Ced <- Nor %>% group_by(Lon, Lat) %>% summarise(Ced_mean = mean(Ced_lib)) 

Mean_BES <- Nor %>% group_by(Lon, Lat) %>% summarise(BES_mean = mean(BES)) 

Mean_BP <-  Nor %>% group_by(Lon, Lat) %>% summarise(BP_mean = mean(Bet_pen))

Mean_BPU <- Nor %>% group_by(Lon, Lat) %>% summarise(BPU_mean = mean(Bet_pub))

Mean_Car <- Nor %>% group_by(Lon, Lat) %>% summarise(Car_mean = mean(Car_bet))

Mean_Cor <- Nor %>% group_by(Lon, Lat) %>% summarise(Cor_mean = mean(Cor_ave))

Mean_Fag <- Nor %>% group_by(Lon, Lat) %>% summarise(Fag_mean = mean(Fag_syl))

Mean_Fra <- Nor %>% group_by(Lon, Lat) %>% summarise(Fra_mean = mean(Fra_exc))

Mean_Jun <- Nor %>% group_by(Lon, Lat) %>% summarise(Jun_mean = mean(Jun_oxy))

Mean_MRS <- Nor %>% group_by(Lon, Lat) %>% summarise(MRS_mean = mean(MRS))

Mean_PA <-  Nor %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Pic_abi))

Mean_AG <-  Nor %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Aln_glu))

Mean_PS <-  Nor %>%  group_by(Lon, Lat) %>% summarise(PS_mean = mean(Pin_syl))

Mean_PN <-  Nor %>%  group_by(Lon, Lat) %>% summarise(PN_mean = mean(Pin_nig))

Mean_PB <-  Nor %>%  group_by(Lon, Lat) %>% summarise(PB_mean = mean(Pin_bru))

Mean_PH <-  Nor %>%  group_by(Lon, Lat) %>% summarise(PH_mean = mean(Pin_hal))

Mean_Pop <- Nor %>% group_by(Lon, Lat) %>% summarise(Pop_mean = mean(Pop_tre))

Mean_QC <-  Nor %>% group_by(Lon, Lat) %>% summarise(QC_mean = mean(Que_coc))

Mean_QI <-  Nor %>% group_by(Lon, Lat) %>% summarise(QI_mean = mean(Que_ile))

Mean_QP <-  Nor %>% group_by(Lon, Lat) %>% summarise(QP_mean = mean(Que_pub))

Mean_QR <-  Nor %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_rob))

Mean_QM <-  Nor %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_mac))

Mean_Til <- Nor %>% group_by(Lon, Lat) %>% summarise(Til_mean = mean(Til_cor))

Mean_Ulm <- Nor %>% group_by(Lon, Lat) %>% summarise(Ulm_mean = mean(Ulm_gla))

Mean_C3 <-  Nor %>% group_by(Lon, Lat) %>% summarise(C3_mean = mean(C3_gr))


Nor_Species <- cbind.data.frame(Mean_Ab[3], Mean_Ced[3], Mean_BES[3], Mean_BP[3], Mean_BPU[3], Mean_Car[3], Mean_Cor[3], Mean_Fag[3], 
                            Mean_Fra[3], Mean_Jun[3], Mean_MRS[3], Mean_PA[3], Mean_AG[3], Mean_PS[3], Mean_PN[3], Mean_PB[3], 
                            Mean_PH[3], Mean_Pop[3], Mean_QC[3], Mean_QI[3], Mean_QP[3], Mean_QR[3], Mean_QM[3], Mean_Til[3], 
                            Mean_Ulm[3], Mean_C3[3])
colnames(Nor_Species) <- c("Abies spp.", "C. libani", "Boreal shrub", "B. pendula", "B. pubescens", "Carpinus spp.", "C. avellana",
                       "Fagus spp.","F. excelsior", "J. oxycedrus", "Med. shrub", "Picea spp.", "A. glutinosa", "P. sylvestris",
                       "P. nigra", "P. brutia", "P. halepensis", "P. tremula", "Q. coccifera", "Q. ilex", "Q. pubescens", 
                       "Q. robur", "Q. macranthera", "T. cordata",  "U. glabra", "Grass")


setwd("C:/Users/EC-Earth/")
gall16 <- read.table("cmass.out", sep="", header=T)
gall16$Year <- as.POSIXct(paste(gall16$Year), format = "%Y")
gall16 <- as.data.frame(gall16)
All_yr3 <- gall16 %>% group_by(year=floor_date(Year, "year"))
Per_yr3 <- All_yr3 %>% mutate(year=year(Year)) %>%
  group_by(Year, year)

EE<- Per_yr3 %>% filter(year <= "1990")
# EE<- Per_yr3 %>% filter(year >= "2071")
# EE <- Per_yr3 %>% filter(year >= "2031" & year <= "2060")
# EE <- Per_yr3 %>% filter(year >= "1991" & year <= "2020")

Mean_Ab <- EE %>% group_by(Lon, Lat) %>% summarise(Ab_mean = mean(Abi_nor)) 

Mean_Ced <- EE %>% group_by(Lon, Lat) %>% summarise(Ced_mean = mean(Ced_lib)) 

Mean_BES <- EE %>% group_by(Lon, Lat) %>% summarise(BES_mean = mean(BES)) 

Mean_BP <-  EE %>% group_by(Lon, Lat) %>% summarise(BP_mean = mean(Bet_pen))

Mean_BPU <- EE %>% group_by(Lon, Lat) %>% summarise(BPU_mean = mean(Bet_pub))

Mean_Car <- EE %>% group_by(Lon, Lat) %>% summarise(Car_mean = mean(Car_bet))

Mean_Cor <- EE %>% group_by(Lon, Lat) %>% summarise(Cor_mean = mean(Cor_ave))

Mean_Fag <- EE %>% group_by(Lon, Lat) %>% summarise(Fag_mean = mean(Fag_syl))

Mean_Fra <- EE %>% group_by(Lon, Lat) %>% summarise(Fra_mean = mean(Fra_exc))

Mean_Jun <- EE %>% group_by(Lon, Lat) %>% summarise(Jun_mean = mean(Jun_oxy))

Mean_MRS <- EE %>% group_by(Lon, Lat) %>% summarise(MRS_mean = mean(MRS))

Mean_PA <-  EE %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Pic_abi))

Mean_AG <-  EE %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Aln_glu))

Mean_PS <-  EE %>%  group_by(Lon, Lat) %>% summarise(PS_mean = mean(Pin_syl))

Mean_PN <-  EE %>%  group_by(Lon, Lat) %>% summarise(PN_mean = mean(Pin_nig))

Mean_PB <-  EE %>%  group_by(Lon, Lat) %>% summarise(PB_mean = mean(Pin_bru))

Mean_PH <-  EE %>%  group_by(Lon, Lat) %>% summarise(PH_mean = mean(Pin_hal))

Mean_Pop <- EE %>% group_by(Lon, Lat) %>% summarise(Pop_mean = mean(Pop_tre))

Mean_QC <-  EE %>% group_by(Lon, Lat) %>% summarise(QC_mean = mean(Que_coc))

Mean_QI <-  EE %>% group_by(Lon, Lat) %>% summarise(QI_mean = mean(Que_ile))

Mean_QP <-  EE %>% group_by(Lon, Lat) %>% summarise(QP_mean = mean(Que_pub))

Mean_QR <-  EE %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_rob))

Mean_QM <-  EE %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_mac))

Mean_Til <- EE %>% group_by(Lon, Lat) %>% summarise(Til_mean = mean(Til_cor))

Mean_Ulm <- EE %>% group_by(Lon, Lat) %>% summarise(Ulm_mean = mean(Ulm_gla))

Mean_C3 <-  EE %>% group_by(Lon, Lat) %>% summarise(C3_mean = mean(C3_gr))


EE_Species <- cbind.data.frame(Mean_Ab[3], Mean_Ced[3], Mean_BES[3], Mean_BP[3], Mean_BPU[3], Mean_Car[3], Mean_Cor[3], Mean_Fag[3], 
                            Mean_Fra[3], Mean_Jun[3], Mean_MRS[3], Mean_PA[3], Mean_AG[3], Mean_PS[3], Mean_PN[3], Mean_PB[3], 
                            Mean_PH[3], Mean_Pop[3], Mean_QC[3], Mean_QI[3], Mean_QP[3], Mean_QR[3], Mean_QM[3], Mean_Til[3], 
                            Mean_Ulm[3], Mean_C3[3])
colnames(EE_Species) <- c("Abies spp.", "C. libani", "Boreal shrub", "B. pendula", "B. pubescens", "Carpinus spp.", "C. avellana",
                       "Fagus spp.","F. excelsior", "J. oxycedrus", "Med. shrub", "Picea spp.", "A. glutinosa", "P. sylvestris",
                       "P. nigra", "P. brutia", "P. halepensis", "P. tremula", "Q. coccifera", "Q. ilex", "Q. pubescens", 
                       "Q. robur", "Q. macranthera", "T. cordata",  "U. glabra", "Grass")


setwd("C:/Users/MPI/")
gallfin <- read.table("cmass.out", sep="", header=T)
gallfin$Year <- as.POSIXct(paste(gallfin$Year), format = "%Y")
gallfin <- as.data.frame(gallfin)

All_yr4 <- gallfin %>% group_by(year=floor_date(Year, "year"))
Per_yr4 <- All_yr4 %>% mutate(year=year(Year)) %>%
  group_by(Year, year)

Per_yr4 <- Per_yr4 %>% filter(year < 2100)

MP <- Per_yr4 %>% filter(year <= "1990")
# MP <- Per_yr4 %>% filter(year >= "2071")
# MP <- Per_yr4 %>% filter(year >= "2031" & year <= "2060")
# MP <- Per_yr4 %>% filter(year >= "1991" & year <= "2020")

Mean_Ab <- MP %>% group_by(Lon, Lat) %>% summarise(Ab_mean = mean(Abi_nor)) 

Mean_Ced <- MP %>% group_by(Lon, Lat) %>% summarise(Ced_mean = mean(Ced_lib)) 

Mean_BES <- MP %>% group_by(Lon, Lat) %>% summarise(BES_mean = mean(BES)) 

Mean_BP <-  MP %>% group_by(Lon, Lat) %>% summarise(BP_mean = mean(Bet_pen))

Mean_BPU <- MP %>% group_by(Lon, Lat) %>% summarise(BPU_mean = mean(Bet_pub))

Mean_Car <- MP %>% group_by(Lon, Lat) %>% summarise(Car_mean = mean(Car_bet))

Mean_Cor <- MP %>% group_by(Lon, Lat) %>% summarise(Cor_mean = mean(Cor_ave))

Mean_Fag <- MP %>% group_by(Lon, Lat) %>% summarise(Fag_mean = mean(Fag_syl))

Mean_Fra <- MP %>% group_by(Lon, Lat) %>% summarise(Fra_mean = mean(Fra_exc))

Mean_Jun <- MP %>% group_by(Lon, Lat) %>% summarise(Jun_mean = mean(Jun_oxy))

Mean_MRS <- MP %>% group_by(Lon, Lat) %>% summarise(MRS_mean = mean(MRS))

Mean_PA <-  MP %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Pic_abi))

Mean_AG <-  MP %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Aln_glu))

Mean_PS <-  MP %>%  group_by(Lon, Lat) %>% summarise(PS_mean = mean(Pin_syl))

Mean_PN <-  MP %>%  group_by(Lon, Lat) %>% summarise(PN_mean = mean(Pin_nig))

Mean_PB <-  MP %>%  group_by(Lon, Lat) %>% summarise(PB_mean = mean(Pin_bru))

Mean_PH <-  MP %>%  group_by(Lon, Lat) %>% summarise(PH_mean = mean(Pin_hal))

Mean_Pop <- MP %>% group_by(Lon, Lat) %>% summarise(Pop_mean = mean(Pop_tre))

Mean_QC <-  MP %>% group_by(Lon, Lat) %>% summarise(QC_mean = mean(Que_coc))

Mean_QI <-  MP %>% group_by(Lon, Lat) %>% summarise(QI_mean = mean(Que_ile))

Mean_QP <-  MP %>% group_by(Lon, Lat) %>% summarise(QP_mean = mean(Que_pub))

Mean_QR <-  MP %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_rob))

Mean_QM <-  MP %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_mac))

Mean_Til <- MP %>% group_by(Lon, Lat) %>% summarise(Til_mean = mean(Til_cor))

Mean_Ulm <- MP %>% group_by(Lon, Lat) %>% summarise(Ulm_mean = mean(Ulm_gla))

Mean_C3 <-  MP %>% group_by(Lon, Lat) %>% summarise(C3_mean = mean(C3_gr))


MP_Species <- cbind.data.frame(Mean_Ab[3], Mean_Ced[3], Mean_BES[3], Mean_BP[3], Mean_BPU[3], Mean_Car[3], Mean_Cor[3], Mean_Fag[3], 
                            Mean_Fra[3], Mean_Jun[3], Mean_MRS[3], Mean_PA[3], Mean_AG[3], Mean_PS[3], Mean_PN[3], Mean_PB[3], 
                            Mean_PH[3], Mean_Pop[3], Mean_QC[3], Mean_QI[3], Mean_QP[3], Mean_QR[3], Mean_QM[3], Mean_Til[3], 
                            Mean_Ulm[3], Mean_C3[3])
colnames(MP_Species) <- c("Abies spp.", "C. libani", "Boreal shrub", "B. pendula", "B. pubescens", "Carpinus spp.", "C. avellana",
                       "Fagus spp.","F. excelsior", "J. oxycedrus", "Med. shrub", "Picea spp.", "A. glutinosa", "P. sylvestris",
                       "P. nigra", "P. brutia", "P. halepensis", "P. tremula", "Q. coccifera", "Q. ilex", "Q. pubescens", 
                       "Q. robur", "Q. macranthera", "T. cordata",  "U. glabra", "Grass")


setwd("C:/Users/INM-CM5/")
gall16 <- read.table("cmass.out", sep="", header=T)
gall16$Year <- as.POSIXct(paste(gall16$Year), format = "%Y")
gall16 <- as.data.frame(gall16)
All_yr5 <- gall16 %>% group_by(year=floor_date(Year, "year"))
Per_yr5 <- All_yr5 %>% mutate(year=year(Year)) %>%
  group_by(Year, year)
Per_yr5 <- Per_yr5 %>% filter(year < 2100)

IN <- Per_yr5 %>% filter(year <= "1990")
# IN <- Per_yr5 %>% filter(year >= "2071")
# IN <- Per_yr5 %>% filter(year >= "2031" & year <= "2060")
# IN <- Per_yr5 %>% filter(year >= "1991" & year <= "2020")


Mean_Ab <- IN %>% group_by(Lon, Lat) %>% summarise(Ab_mean = mean(Abi_nor)) 

Mean_Ced <- IN %>% group_by(Lon, Lat) %>% summarise(Ced_mean = mean(Ced_lib)) 

Mean_BES <- IN %>% group_by(Lon, Lat) %>% summarise(BES_mean = mean(BES)) 

Mean_BP <-  IN %>% group_by(Lon, Lat) %>% summarise(BP_mean = mean(Bet_pen))

Mean_BPU <- IN %>% group_by(Lon, Lat) %>% summarise(BPU_mean = mean(Bet_pub))

Mean_Car <- IN %>% group_by(Lon, Lat) %>% summarise(Car_mean = mean(Car_bet))

Mean_Cor <- IN %>% group_by(Lon, Lat) %>% summarise(Cor_mean = mean(Cor_ave))

Mean_Fag <- IN %>% group_by(Lon, Lat) %>% summarise(Fag_mean = mean(Fag_syl))

Mean_Fra <- IN %>% group_by(Lon, Lat) %>% summarise(Fra_mean = mean(Fra_exc))

Mean_Jun <- IN %>% group_by(Lon, Lat) %>% summarise(Jun_mean = mean(Jun_oxy))

Mean_MRS <- IN %>% group_by(Lon, Lat) %>% summarise(MRS_mean = mean(MRS))

Mean_PA <-  IN %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Pic_abi))

Mean_AG <-  IN %>%  group_by(Lon, Lat) %>% summarise(PA = mean(Aln_glu))

Mean_PS <-  IN %>%  group_by(Lon, Lat) %>% summarise(PS_mean = mean(Pin_syl))

Mean_PN <-  IN %>%  group_by(Lon, Lat) %>% summarise(PN_mean = mean(Pin_nig))

Mean_PB <-  IN %>%  group_by(Lon, Lat) %>% summarise(PB_mean = mean(Pin_bru))

Mean_PH <-  IN %>%  group_by(Lon, Lat) %>% summarise(PH_mean = mean(Pin_hal))

Mean_Pop <- IN %>% group_by(Lon, Lat) %>% summarise(Pop_mean = mean(Pop_tre))

Mean_QC <-  IN %>% group_by(Lon, Lat) %>% summarise(QC_mean = mean(Que_coc))

Mean_QI <-  IN %>% group_by(Lon, Lat) %>% summarise(QI_mean = mean(Que_ile))

Mean_QP <-  IN %>% group_by(Lon, Lat) %>% summarise(QP_mean = mean(Que_pub))

Mean_QR <-  IN %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_rob))

Mean_QM <-  IN %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_mac))

Mean_Til <- IN %>% group_by(Lon, Lat) %>% summarise(Til_mean = mean(Til_cor))

Mean_Ulm <- IN %>% group_by(Lon, Lat) %>% summarise(Ulm_mean = mean(Ulm_gla))

Mean_C3 <-  IN %>% group_by(Lon, Lat) %>% summarise(C3_mean = mean(C3_gr))


IN_Species <- cbind.data.frame(Mean_Ab[3], Mean_Ced[3], Mean_BES[3], Mean_BP[3], Mean_BPU[3], Mean_Car[3], Mean_Cor[3], Mean_Fag[3], 
                            Mean_Fra[3], Mean_Jun[3], Mean_MRS[3], Mean_PA[3], Mean_AG[3], Mean_PS[3], Mean_PN[3], Mean_PB[3], 
                            Mean_PH[3], Mean_Pop[3], Mean_QC[3], Mean_QI[3], Mean_QP[3], Mean_QR[3], Mean_QM[3], Mean_Til[3], 
                            Mean_Ulm[3], Mean_C3[3])
colnames(IN_Species) <- c("Abies spp.", "C. libani", "Boreal shrub", "B. pendula", "B. pubescens", "Carpinus spp.", "C. avellana",
                       "Fagus spp.","F. excelsior", "J. oxycedrus", "Med. shrub", "Picea spp.", "A. glutinosa", "P. sylvestris",
                       "P. nigra", "P. brutia", "P. halepensis", "P. tremula", "Q. coccifera", "Q. ilex", "Q. pubescens", 
                       "Q. robur", "Q. macranthera", "T. cordata",  "U. glabra", "Grass")


Abies <- cbind.data.frame(IN_Species$`Abies spp.`, Nor_Species$`Abies spp.`, EE_Species$`Abies spp.`, CM_Species$`Abies spp.`, MP_Species$`Abies spp.`)
Abies$row_median = rowMedians(as.matrix(Abies))

Cedrus <- cbind.data.frame(IN_Species$`C. libani`, Nor_Species$`C. libani`, EE_Species$`C. libani`, CM_Species$`C. libani`, MP_Species$`C. libani`)
Cedrus$row_median = rowMedians(as.matrix(Cedrus))

BES <- cbind.data.frame(IN_Species$`Boreal shrub`, Nor_Species$`Boreal shrub`, EE_Species$`Boreal shrub`, CM_Species$`Boreal shrub`, MP_Species$`Boreal shrub`)
BES$row_median = rowMedians(as.matrix(BES))

Bet_pen <- cbind.data.frame(IN_Species$`B. pendula`, Nor_Species$`B. pendula`, EE_Species$`B. pendula`, CM_Species$`B. pendula`, MP_Species$`B. pendula`)
Bet_pen$row_median = rowMedians(as.matrix(Bet_pen))

Bet_pub <- cbind.data.frame(IN_Species$`B. pubescens`, Nor_Species$`B. pubescens`, EE_Species$`B. pubescens`, CM_Species$`B. pubescens`, MP_Species$`B. pubescens`)
Bet_pub$row_median = rowMedians(as.matrix(Bet_pub))

Carpinus <- cbind.data.frame(IN_Species$`Carpinus spp.`, Nor_Species$`Carpinus spp.`, EE_Species$`Carpinus spp.`, CM_Species$`Carpinus spp.`, MP_Species$`Carpinus spp.`) 
Carpinus$row_median = rowMedians(as.matrix(Carpinus))

Corylus <- cbind.data.frame(IN_Species$`C. avellana`, Nor_Species$`C. avellana`, EE_Species$`C. avellana`, CM_Species$`C. avellana`, MP_Species$`C. avellana`) 
Corylus$row_median = rowMedians(as.matrix(Corylus))

Fagus <- cbind.data.frame(IN_Species$`Fagus spp.`, Nor_Species$`Fagus spp.`, EE_Species$`Fagus spp.`, CM_Species$`Fagus spp.`, MP_Species$`Fagus spp.`)
Fagus$row_median = rowMedians(as.matrix(Fagus))

Frax <- cbind.data.frame(IN_Species$`F. excelsior`, Nor_Species$`F. excelsior`, EE_Species$`F. excelsior`, CM_Species$`F. excelsior`, MP_Species$`F. excelsior`)
Frax$row_median = rowMedians(as.matrix(Frax))

Juniper <- cbind.data.frame(IN_Species$`J. oxycedrus`, Nor_Species$`J. oxycedrus`, EE_Species$`J. oxycedrus`, CM_Species$`J. oxycedrus`, MP_Species$`J. oxycedrus`)
Juniper$row_median = rowMedians(as.matrix(Juniper))

MRS <- cbind.data.frame(IN_Species$`Med. shrub`, Nor_Species$`Med. shrub`, EE_Species$`Med. shrub`, CM_Species$`Med. shrub`, MP_Species$`Med. shrub`)
MRS$row_median = rowMedians(as.matrix(MRS))

Picea <- cbind.data.frame(IN_Species$`Picea spp.`, Nor_Species$`Picea spp.`, EE_Species$`Picea spp.`, CM_Species$`Picea spp.`, MP_Species$`Picea spp.`)
Picea$row_median = rowMedians(as.matrix(Picea))

Alnus <- cbind.data.frame(IN_Species$`A. glutinosa`, Nor_Species$`A. glutinosa`, EE_Species$`A. glutinosa`, CM_Species$`A. glutinosa`, MP_Species$`A. glutinosa`)
Alnus$row_median = rowMedians(as.matrix(Alnus))

PS <- cbind.data.frame(IN_Species$`P. sylvestris`, Nor_Species$`P. sylvestris`, EE_Species$`P. sylvestris`, CM_Species$`P. sylvestris`, MP_Species$`P. sylvestris`)
PS$row_median = rowMedians(as.matrix(PS))

PN <- cbind.data.frame(IN_Species$`P. nigra`, Nor_Species$`P. nigra`, EE_Species$`P. nigra`, CM_Species$`P. nigra`, MP_Species$`P. nigra`)
PN$row_median = rowMedians(as.matrix(PN))

PH <-  cbind.data.frame(IN_Species$`P. halepensis`, Nor_Species$`P. halepensis`, EE_Species$`P. halepensis`, CM_Species$`P. halepensis`, MP_Species$`P. halepensis`)
PH$row_median = rowMedians(as.matrix(PH))

PB <- cbind.data.frame(IN_Species$`P. brutia`, Nor_Species$`P. brutia`, EE_Species$`P. brutia`, CM_Species$`P. brutia`, MP_Species$`P. brutia`)
PB$row_median = rowMedians(as.matrix(PB))

Populus <- cbind.data.frame(IN_Species$`P. tremula`, Nor_Species$`P. tremula`, EE_Species$`P. tremula`, CM_Species$`P. tremula`, MP_Species$`P. tremula`)
Populus$row_median = rowMedians(as.matrix(Populus))

QC <- cbind.data.frame(IN_Species$`Q. coccifera`, Nor_Species$`Q. coccifera`, EE_Species$`Q. coccifera`, CM_Species$`Q. coccifera`, MP_Species$`Q. coccifera`)
QC$row_median = rowMedians(as.matrix(QC))

QI <- cbind.data.frame(IN_Species$`Q. ilex`, Nor_Species$`Q. ilex`, EE_Species$`Q. ilex`, CM_Species$`Q. ilex`, MP_Species$`Q. ilex`)
QI$row_median = rowMedians(as.matrix(QI))

QP <- cbind.data.frame(IN_Species$`Q. pubescens`, Nor_Species$`Q. pubescens`, EE_Species$`Q. pubescens`, CM_Species$`Q. pubescens`, MP_Species$`Q. pubescens`)
QP$row_median = rowMedians(as.matrix(QP))

QR <- cbind.data.frame(IN_Species$`Q. robur`, Nor_Species$`Q. robur`, EE_Species$`Q. robur`, CM_Species$`Q. robur`, MP_Species$`Q. robur`)
QR$row_median = rowMedians(as.matrix(QR))

QM <- cbind.data.frame(IN_Species$`Q. macranthera`, Nor_Species$`Q. macranthera`, EE_Species$`Q. macranthera`, CM_Species$`Q. macranthera`, MP_Species$`Q. macranthera`) 
QM$row_median = rowMedians(as.matrix(QM))

Tilia <- cbind.data.frame(IN_Species$`T. cordata`, Nor_Species$`T. cordata`, EE_Species$`T. cordata`, CM_Species$`T. cordata`, MP_Species$`T. cordata`)
Tilia$row_median = rowMedians(as.matrix(Tilia))

Ulmus <- cbind.data.frame(IN_Species$`U. glabra`, Nor_Species$`U. glabra`, EE_Species$`U. glabra`, CM_Species$`U. glabra`, MP_Species$`U. glabra`)
Ulmus$row_median = rowMedians(as.matrix(Ulmus))

C3 <- cbind.data.frame(IN_Species$`Grass`, Nor_Species$`Grass`, EE_Species$`Grass`, CM_Species$`Grass`, MP_Species$`Grass`)
C3$row_median = rowMedians(as.matrix(C3))



Ensemble <- cbind.data.frame(Mean_Ulm[1:2], Abies[6], Cedrus[6], BES[6], Bet_pen[6], Bet_pub[6], Carpinus[6], Corylus[6], Fagus[6], 
                             Frax[6], Juniper[6], MRS[6], Picea[6], Alnus[6], 
                             PS[6], PN[6], PB[6], PH[6], Populus[6], QC[6], QI[6], QP[6], QR[6], QM[6], 
                             Tilia[6], Ulmus[6], C3[6])


################################################################################



BM_sp <- Ensemble[c(3:28)] #taking only the CMASS fields
colnames(BM_sp) <- c("Abies spp.", "C. libani", "Boreal shrub", "B. pendula", "B. pubescens", "Carpinus spp.", "C. avellana",
                      "Fagus spp.","F. excelsior", "J. oxycedrus", "Med. shrub", "Picea spp.", "A. glutinosa", "P. sylvestris",
                      "P. nigra", "P. brutia", "P. halepensis", "P. tremula", "Q. coccifera", "Q. ilex", "Q. pubescens", 
                      "Q. robur", "Q. macranthera", "T. cordata",  "U. glabra", "Grass")


BM_spm <- reshape2::melt(BM_sp) #reshaping the dataframe


theme_set(theme_minimal())

BM_box <- ggplot(BM_spm, aes(x = variable, y = value, color = variable)) +           
  geom_boxplot(fatten = 7.5)+ # fatten = 4   lwd=1
  # labs(title =paste("1961-1990 mean biomass per species from ensemble median"))+
  # labs(title =paste("1991-2020 mean biomass per species from ensemble median"))+
  labs(title =paste("2031-2060 mean biomass per species from ensemble median"))+
  # labs(title =paste("2071-2100 mean biomass per species from ensemble median"))+
  
  scale_color_manual(values = c("Abies spp."= "#799c11","C. libani" = "#b7c50e","Boreal shrub" = "blue", "B. pendula"="red","B. pubescens"="steel blue",
                                "Carpinus spp."="#90eb10", "C. avellana" =  "#0cb8e8", "Fagus spp."="#134d19", "F. excelsior" = "#b595e6", 
                                "J. oxycedrus" = "#207852", "Med. shrub" = "#122da6", "Picea spp."="#e02b73",  "A. glutinosa" = "#F87F65", 
                                "P. sylvestris"="#8d49cc", "P. nigra" = "#1c83ba", "P. brutia" = "grey", "P. halepensis"="#f772b0", "P. tremula"= "#b7eb34",
                                "Q. coccifera"= "#e8550c", "Q. ilex"="#57c27e", "Q. pubescens" = "#e8940c", "Q. robur"="#DB6C05", "Q. macranthera" = "#fcba03",
                                "T. cordata"= "#ad1320", "U. glabra" = "#611975", "Grass" = "#89dc69"))+ 
  ylab("")+
  xlab("")+
  scale_y_continuous(limits=c(0,13.5))+
  theme(legend.position="none", 
        plot.title = element_text(size=25, face="bold"),
        axis.text.x=element_text(size=30, angle = -60, hjust = 0, vjust = 0), axis.text.y = element_text(size=30))

print(BM_box)


setwd("C:/Users/ITU/Documents/DATA/GCMs/LPJ_GUESS runs/ENSEMBLE_local_AB/New_Dominance/")


png(file=paste("Biomass_1961-1990.png"), width=1980, height=932)
# png(file=paste("Biomass_1991-2020.png"), width=1980, height=932)
# png(file=paste("Biomass_2031-2060.png"), width=1980, height=932)
# png(file=paste("Biomass_2071_2100.png"), width=1980, height=932)

print(BM_box)
dev.off()