Raster4GIS_and_SAD.R · LPJ-GUESS-OutputPlots

##########################################################################################################################
#
# This code complements the bivariate spatio-temporal analysis done for the article: 
# Bikem Ekberzade et al., Latitude or altitude will save us? A centennial case for forests in Asia Minor and its immediate 
# surroundings
# 
# It builds a large dataframe from dens.out and cmass.out files from multiple LPJ-GUESS simulations, calculates the mean
# of certain time periods, calculates the median of period means, calculates Sum of Absolute Deviations (SAD) for projections
# (time periods involving the scenario phase, and converts the final output to a .csv file to run spatial analysis in 
# ArcGIS Pro (or other GIS software)
# 
###########################################################################################################################

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)


setwd("C:/Users/CMCC-ESM2")   # Directory should contain the results of the specified GCM simulation

gall16DENS <- read.table("dens.out", sep="", header=T)
gall16DENS_K <- cbind.data.frame(gall16DENS[1:3], 10000*(gall16DENS[4:30]))   # build dataframe and calculate density/ha
gall16BM <- read.table("cmass.out", sep="", header=T)

##### run this code for each df changing gall16 for biomass or density depending on which df you are building 
##### for all GCMs first (later you will call and take their median):

gall16 <- gall16BM     #or gall16DENS_K

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)    # EC-Earth final year is 2099
Per_yr1$T_C <- Per_yr1$Total - Per_yr1$C3_gr
Per_yr1$shrubs <- (Per_yr1$MRS + Per_yr1$BES + Per_yr1$Cor_ave + Per_yr1$Que_coc + Per_yr1$Jun_oxy)
Per_yr1$temperate <- (Per_yr1$Abi_nor + Per_yr1$Bet_pen + Per_yr1$Car_bet + Per_yr1$Ced_lib + Per_yr1$Cor_ave + Per_yr1$Fag_syl +
                        Per_yr1$Fra_exc + Per_yr1$Jun_oxy + Per_yr1$Pin_bru + Per_yr1$Pin_nig + Per_yr1$Pin_hal + Per_yr1$Pop_tre +
                        Per_yr1$Que_coc + Per_yr1$Que_ile + Per_yr1$Que_mac + Per_yr1$Que_pub + Per_yr1$Que_rob + Per_yr1$Aln_glu +
                        Per_yr1$Til_cor + Per_yr1$Ulm_gla + Per_yr1$MRS)

Per_yr1$boreal <- (Per_yr1$BES + Per_yr1$Pin_syl + Per_yr1$Pic_abi + Per_yr1$Bet_pub)

# un-commentout for which time period you are running the script for:

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))

Mean_Tot <- CM %>% group_by(Lon, Lat) %>% summarise(Tot_mean = mean(Total))

Mean_TC <- CM %>% group_by(Lon, Lat) %>% summarise(TC_mean = mean(T_C))

Mean_shrubs <- CM %>% group_by(Lon, Lat) %>% summarise(Shrub_mean = mean(shrubs))

Mean_temp <- CM %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(temperate))

Mean_bor <- CM %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(boreal))


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_shrubs[3], Mean_temp[3], Mean_bor[3], Mean_C3[3], Mean_Tot[3], Mean_TC[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", "Shrubs", "Temperate", "Boreal","Grass", "Total", "T-C")

################################## NORESM ##############################

setwd("C:/Users/NorESM2-MM/")
gall17DENS <- read.table("dens.out", sep="", header=T)
gall17DENS_K <- cbind.data.frame(gall17DENS[1:3], 10000*(gall17DENS[4:30]))
gall17BM <- read.table("cmass.out", sep="", header=T)

gall17 <- gall17BM # this changes depending (see first GCM)

gall17$Year <- as.POSIXct(paste(gall17$Year), format = "%Y")
gall17 <- as.data.frame(gall17)
All_yr2 <- gall17 %>% 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)
Per_yr2$T_C <- Per_yr2$Total - Per_yr2$C3_gr
Per_yr2$shrubs <- (Per_yr2$MRS + Per_yr2$BES + Per_yr2$Cor_ave + Per_yr2$Que_coc + Per_yr2$Jun_oxy)
Per_yr2$temperate <- (Per_yr2$Abi_nor + Per_yr2$Bet_pen + Per_yr2$Car_bet + Per_yr2$Ced_lib + Per_yr2$Cor_ave + Per_yr2$Fag_syl +
                        Per_yr2$Fra_exc + Per_yr2$Jun_oxy + Per_yr2$Pin_bru + Per_yr2$Pin_nig + Per_yr2$Pin_hal + Per_yr2$Pop_tre +
                        Per_yr2$Que_coc + Per_yr2$Que_ile + Per_yr2$Que_mac + Per_yr2$Que_pub + Per_yr2$Que_rob + Per_yr2$Aln_glu +
                        Per_yr2$Til_cor + Per_yr2$Ulm_gla + Per_yr2$MRS)

Per_yr2$boreal <- (Per_yr2$BES + Per_yr2$Pin_syl + Per_yr2$Pic_abi + Per_yr2$Bet_pub)


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))

Mean_Tot <- Nor %>% group_by(Lon, Lat) %>% summarise(Tot_mean = mean(Total))

Mean_TC <- Nor %>% group_by(Lon, Lat) %>% summarise(TC_mean = mean(T_C))

Mean_shrubs <- Nor %>% group_by(Lon, Lat) %>% summarise(Shrub_mean = mean(shrubs))

Mean_temp <- Nor %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(temperate))

Mean_bor <- Nor %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(boreal))

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_shrubs[3], Mean_temp[3], Mean_bor[3], Mean_C3[3], Mean_Tot[3], Mean_TC[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", "Shrubs", "Temperate", "Boreal","Grass", "Total", "T-C")

################################## EC-Earth ####################################


setwd("C:/Users/EC-Earth-3P-HR/")
gall18DENS <- read.table("dens.out", sep="", header=T)
gall18DENS_K <- cbind.data.frame(gall18DENS[1:3], 10000*(gall18DENS[4:30]))
gall18BM <- read.table("cmass.out", sep="", header=T)

gall18 <- gall18BM     #this changes to DENS_K or BM, see first GCM

gall18$Year <- as.POSIXct(paste(gall18$Year), format = "%Y")
gall18 <- as.data.frame(gall18)
All_yr3 <- gall18 %>% group_by(year=floor_date(Year, "year"))

Per_yr3 <- All_yr3 %>% mutate(year=year(Year)) %>%
  group_by(Year, year)

Per_yr3$T_C <- Per_yr3$Total - Per_yr3$C3_gr
Per_yr3$shrubs <- (Per_yr3$MRS + Per_yr3$BES + Per_yr3$Cor_ave + Per_yr3$Que_coc + Per_yr3$Jun_oxy)
Per_yr3$temperate <- (Per_yr3$Abi_nor + Per_yr3$Bet_pen + Per_yr3$Car_bet + Per_yr3$Ced_lib + Per_yr3$Cor_ave + Per_yr3$Fag_syl +
                        Per_yr3$Fra_exc + Per_yr3$Jun_oxy + Per_yr3$Pin_bru + Per_yr3$Pin_nig + Per_yr3$Pin_hal + Per_yr3$Pop_tre +
                        Per_yr3$Que_coc + Per_yr3$Que_ile + Per_yr3$Que_mac + Per_yr3$Que_pub + Per_yr3$Que_rob + Per_yr3$Aln_glu +
                        Per_yr3$Til_cor + Per_yr3$Ulm_gla + Per_yr3$MRS)

Per_yr3$boreal <- (Per_yr3$BES + Per_yr3$Pin_syl + Per_yr3$Pic_abi + Per_yr3$Bet_pub)

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))

Mean_Tot <- EE %>% group_by(Lon, Lat) %>% summarise(Tot_mean = mean(Total))

Mean_TC <- EE %>% group_by(Lon, Lat) %>% summarise(TC_mean = mean(T_C))

Mean_shrubs <- EE %>% group_by(Lon, Lat) %>% summarise(Shrub_mean = mean(shrubs))

Mean_temp <- EE %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(temperate))

Mean_bor <- EE %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(boreal))

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_shrubs[3], Mean_temp[3], Mean_bor[3], Mean_C3[3], Mean_Tot[3], Mean_TC[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", "Shrubs", "Temperate", "Boreal","Grass", "Total", "T-C")


################################## INM #############################################

setwd("C:/Users/INM-CM5/")
gall19DENS <- read.table("dens.out", sep="", header=T)
gall19DENS_K <- cbind.data.frame(gall19DENS[1:3], 10000*(gall19DENS[4:30]))
gall19BM <- read.table("cmass.out", sep="", header=T)


gall19 <- gall19BM #this changes to DENS_K or BM, see first GCM

gall19$Year <- as.POSIXct(paste(gall19$Year), format = "%Y")
gall19 <- as.data.frame(gall19)
All_yr5 <- gall19 %>% 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)

Per_yr5$T_C <- Per_yr5$Total - Per_yr5$C3_gr
Per_yr5$shrubs <- (Per_yr5$MRS + Per_yr5$BES + Per_yr5$Cor_ave + Per_yr5$Que_coc + Per_yr5$Jun_oxy)
Per_yr5$temperate <- (Per_yr5$Abi_nor + Per_yr5$Bet_pen + Per_yr5$Car_bet + Per_yr5$Ced_lib + Per_yr5$Cor_ave + Per_yr5$Fag_syl +
                        Per_yr5$Fra_exc + Per_yr5$Jun_oxy + Per_yr5$Pin_bru + Per_yr5$Pin_nig + Per_yr5$Pin_hal + Per_yr5$Pop_tre +
                        Per_yr5$Que_coc + Per_yr5$Que_ile + Per_yr5$Que_mac + Per_yr5$Que_pub + Per_yr5$Que_rob + Per_yr5$Aln_glu +
                        Per_yr5$Til_cor + Per_yr5$Ulm_gla + Per_yr5$MRS)

Per_yr5$boreal <- (Per_yr5$BES + Per_yr5$Pin_syl + Per_yr5$Pic_abi + Per_yr5$Bet_pub)


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))

Mean_Tot <- IN %>% group_by(Lon, Lat) %>% summarise(Tot_mean = mean(Total))

Mean_TC <- IN %>% group_by(Lon, Lat) %>% summarise(TC_mean = mean(T_C))

Mean_shrubs <- IN %>% group_by(Lon, Lat) %>% summarise(Shrub_mean = mean(shrubs))

Mean_temp <- IN %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(temperate))

Mean_bor <- IN %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(boreal))


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_shrubs[3], Mean_temp[3], Mean_bor[3], Mean_C3[3], Mean_Tot[3], Mean_TC[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", "Shrubs", "Temperate", "Boreal","Grass", "Total", "T-C")

################################## MPI ####################################################

setwd("C:/Users/MPI/")
gall12 <- read.table("cmass.out", sep="", header=T)
setwd("C:/Users/ITU/Documents/DATA/GCMs/LPJ_GUESS runs/MPI/MPI_local_AB/MPI_local_additional_AB/")
gall14 <- read.table("cmass.out", sep="", header=T)
gall14 <- gall14 %>% filter(Lon >= 49.7)

gallBM <- rbind.data.frame(gall12, gall14)

setwd("C:/Users/ITU/Documents/DATA/GCMs/LPJ_GUESS runs/MPI/MPI_local_AB/")
gall10 <- read.table("dens.out", sep="", header=T)
setwd("C:/Users/ITU/Documents/DATA/GCMs/LPJ_GUESS runs/MPI/MPI_local_AB/MPI_local_additional_AB/")
gall11 <- read.table("dens.out", sep="", header=T)
gall11 <- gall14 %>% filter(Lon >= 49.7)

gallDENS <- rbind.data.frame(gall10, gall11)
gallDENS_K <- cbind.data.frame(gallDENS[1:3], 10000*(gallDENS[4:30]))


gallfin <- gall19BM # this changes to DENS_K or BM, see first GCM

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$T_C <- Per_yr4$Total - Per_yr4$C3_gr
Per_yr4$shrubs <- (Per_yr4$MRS + Per_yr4$BES + Per_yr4$Cor_ave + Per_yr4$Que_coc + Per_yr4$Jun_oxy)
Per_yr4$temperate <- (Per_yr4$Abi_nor + Per_yr4$Bet_pen + Per_yr4$Car_bet + Per_yr4$Ced_lib + Per_yr4$Cor_ave + Per_yr4$Fag_syl +
                        Per_yr4$Fra_exc + Per_yr4$Jun_oxy + Per_yr4$Pin_bru + Per_yr4$Pin_nig + Per_yr4$Pin_hal + Per_yr4$Pop_tre +
                        Per_yr4$Que_coc + Per_yr4$Que_ile + Per_yr4$Que_mac + Per_yr4$Que_pub + Per_yr4$Que_rob + Per_yr4$Aln_glu +
                        Per_yr4$Til_cor + Per_yr4$Ulm_gla + Per_yr4$MRS)

Per_yr4$boreal <- (Per_yr4$BES + Per_yr4$Pin_syl + Per_yr4$Pic_abi + Per_yr4$Bet_pub)

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))

Mean_Tot <- MP %>% group_by(Lon, Lat) %>% summarise(Tot_mean = mean(Total))

Mean_TC <- MP %>% group_by(Lon, Lat) %>% summarise(TC_mean = mean(T_C))

Mean_shrubs <- MP %>% group_by(Lon, Lat) %>% summarise(Shrub_mean = mean(shrubs))

Mean_temp <- MP %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(temperate))

Mean_bor <- MP %>% group_by(Lon, Lat) %>% summarise(Temp_mean = mean(boreal))


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_shrubs[3], Mean_temp[3], Mean_bor[3], Mean_C3[3], Mean_Tot[3], Mean_TC[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", "Shrubs", "Temperate", "Boreal","Grass", "Total", "T-C")

################################## Calculate medians ######################################################################

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$AB_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$CED_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$BES_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$BPE_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$BPU_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$CAR_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$COR_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$FAG_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$FRA_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$JUN_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$MRS_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$PA_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$AG_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$PS_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$PN_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$PH_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$PB_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$POP_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$QC_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$QI_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$QP_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$QR_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$QM_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$TIL_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$ULM_median = rowMedians(as.matrix(Ulmus))

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

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

TC <- cbind.data.frame(IN_Species$`T-C`, Nor_Species$`T-C`, EE_Species$`T-C`, CM_Species$`T-C`, MP_Species$`T-C`)
TC$TC_median = rowMedians(as.matrix(TC))

Shrub <- cbind.data.frame(IN_Species$`Shrubs`, Nor_Species$`Shrubs`, EE_Species$`Shrubs`, CM_Species$`Shrubs`, MP_Species$`Shrubs`)
Shrub$SHR_median = rowMedians(as.matrix(Shrub))

Temp <- cbind.data.frame(IN_Species$`Temperate`, Nor_Species$`Temperate`, EE_Species$`Temperate`, CM_Species$`Temperate`, MP_Species$`Temperate`)
Temp$TMP_median = rowMedians(as.matrix(Temp))

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


##############################################################################################
######################### Build dataframes for specified time periods : 1961-1990, ###########
######################### 1991-2020, 2031-2060, 2071-2100.                         ###########
##############################################################################################


 Ensemble_BM1961 <- 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],
                                     Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
 
 
 colnames(Ensemble_BM1961) <- c("Lon","Lat","Abies spp. BM", "C. libani BM", "Boreal shrub BM", "B. pendula BM", "B. pubescens BM", "Carpinus spp. BM", "C. avellana BM", "Fagus spp. BM",
                                "F. excelsior BM", "J. oxycedrus BM", "Med. shrub BM", "Picea spp. BM", "A. glutinosa BM", "P. sylvestris BM", "P. nigra BM", "P. brutia BM", "P. halepensis BM",
                                "P. tremula BM", "Q. coccifera BM", "Q. ilex BM", "Q. pubescens BM", "Q. robur BM", "Q. macranthera BM", "T. cordata BM", "U. glabra BM",
                                "Grass_BM", "Total_BM", "T-C_BM", "Shrubs_BM", "Temperate_BM", "Boreal_BM")
 
# Ensemble_BM2070 <- 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],
#                              Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
# 
# 
# colnames(Ensemble_BM2070) <- c("Lon","Lat","Abies spp. BM", "C. libani BM", "Boreal shrub BM", "B. pendula BM", "B. pubescens BM", "Carpinus spp. BM", "C. avellana BM", "Fagus spp. BM",
#                         "F. excelsior BM", "J. oxycedrus BM", "Med. shrub BM", "Picea spp. BM", "A. glutinosa BM", "P. sylvestris BM", "P. nigra BM", "P. brutia BM", "P. halepensis BM",
#                         "P. tremula BM", "Q. coccifera BM", "Q. ilex BM", "Q. pubescens BM", "Q. robur BM", "Q. macranthera BM", "T. cordata BM", "U. glabra BM",
#                         "Grass_BM", "Total_BM", "T-C_BM", "Shrubs_BM", "Temperate_BM", "Boreal_BM")

# Ensemble_BM2031 <- 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],
#                                   Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
#
#
#colnames(Ensemble_BM2031) <- c("Lon","Lat","Abies spp. BM", "C. libani BM", "Boreal shrub BM", "B. pendula BM", "B. pubescens BM", "Carpinus spp. BM", "C. avellana BM", "Fagus spp. BM",
#                               "F. excelsior BM", "J. oxycedrus BM", "Med. shrub BM", "Picea spp. BM", "A. glutinosa BM", "P. sylvestris BM", "P. nigra BM", "P. brutia BM", "P. halepensis BM",
#                               "P. tremula BM", "Q. coccifera BM", "Q. ilex BM", "Q. pubescens BM", "Q. robur BM", "Q. macranthera BM", "T. cordata BM", "U. glabra BM",
#                               "Grass_BM", "Total_BM", "T-C_BM", "Shrubs_BM", "Temperate_BM", "Boreal_BM")

# Ensemble_BM1991 <- 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],
#                                     Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
# 
# 
# colnames(Ensemble_BM1991) <- c("Lon","Lat","Abies spp. BM", "C. libani BM", "Boreal shrub BM", "B. pendula BM", "B. pubescens BM", "Carpinus spp. BM", "C. avellana BM", "Fagus spp. BM",
#                                "F. excelsior BM", "J. oxycedrus BM", "Med. shrub BM", "Picea spp. BM", "A. glutinosa BM", "P. sylvestris BM", "P. nigra BM", "P. brutia BM", "P. halepensis BM",
#                                "P. tremula BM", "Q. coccifera BM", "Q. ilex BM", "Q. pubescens BM", "Q. robur BM", "Q. macranthera BM", "T. cordata BM", "U. glabra BM",
#                                "Grass_BM", "Total_BM", "T-C_BM", "Shrubs_BM", "Temperate_BM", "Boreal_BM")

# Ensemble_DENS1961 <- 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],
#                              Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
# 
# 
# colnames(Ensemble_DENS1961) <- c("Lon","Lat","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", "Total", "T-C", "Shrubs", "Temperate", "Boreal")

# Ensemble_DENS1991 <- 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],
#                              Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
# 
# 
# colnames(Ensemble_DENS1991) <- c("Lon","Lat","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", "Total", "T-C", "Shrubs", "Temperate", "Boreal")

# Ensemble_DENS2031 <- 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],
#                                       Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
# 
# 
# colnames(Ensemble_DENS2031) <- c("Lon","Lat","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", "Total", "T-C", "Shrubs", "Temperate", "Boreal")


# Ensemble_DENS2070 <- 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],
#                                       Grass[6], Total[6], TC[6], Shrub[6], Temp[6], Bor[6])
# 
# 
# colnames(Ensemble_DENS2070) <- c("Lon","Lat","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", "Total", "T-C", "Shrubs", "Temperate", "Boreal")


##############################################################################################
######################### Merge and write .csv for GIS analysis for : 1961-1990,   ###########
######################### 1991-2020, 2031-2060, 2071-2100.                         ###########
##############################################################################################

setwd("C:/Users/")
dir.create("ENSEMBLE")
setwd("C:/Users/ENSEMBLE")

########## Density and Biomass for beginning and end, respectively, written as separate .csv files for back up
########## alternately, you may run all the codes above, create all the necessary Ensemble datasets, and then just 
########## call them from memory to merge as below (Line 782):

write.csv(Ensemble_DENS1961, "DENS_ENS_1961.csv", row.names = FALSE)
write.csv(Ensemble_DENS2070, "DENS_ENS_2070.csv", row.names = FALSE)

write.csv(Ensemble_BM1961, "BM_ENS_1961.csv", row.names = FALSE)
write.csv(Ensemble_BM2070, "CMASS_ENS_2070.csv", row.names = FALSE)


Merged <- cbind.data.frame(Ensemble_DENS1961, Ensemble_BM1961[3:33])
write.csv(Merged, "DENS_CMASS_ENS1961.csv", row.names = FALSE)

# Merged <- cbind.data.frame(Ensemble_DENS2070, Ensemble_BM2070[3:33])
# write.csv(Merged, "DENS_CMASS_ENS2070.csv", row.names = FALSE)


###### Sum of Absolute Deviations for forecast period means:

# SAD <- (abs(Ensemble_DENS1991[3:30] - CM_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS1991[3:30] - EE_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS1991[3:30] - IN_Species[c(1:25, 29:31)])
#         + abs(Ensemble_DENS1991[3:30] - MP_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS1991[3:30] - Nor_Species[c(1:25, 29:31)]))
# SAD <- (abs(Ensemble_DENS2031[3:30] - CM_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS2031[3:30] - EE_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS2031[3:30] - IN_Species[c(1:25, 29:31)])
#         + abs(Ensemble_DENS2031[3:30] - MP_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS2031[3:30] - Nor_Species[c(1:25, 29:31)]))
# SAD <- (abs(Ensemble_DENS2070[3:30] - CM_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS2070[3:30] - EE_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS2070[3:30] - IN_Species[c(1:25, 29:31)])
#         + abs(Ensemble_DENS2070[3:30] - MP_Species[c(1:25, 29:31)]) + abs(Ensemble_DENS2070[3:30] - Nor_Species[c(1:25, 29:31)]))
 

# SAD <- (abs(Ensemble_BM1991[3:30] - CM_Species[c(1:25, 29:31)]) + abs(Ensemble_BM1991[3:30] - EE_Species[c(1:25, 29:31)]) + abs(Ensemble_BM1991[3:30] - IN_Species[c(1:25, 29:31)])
#         + abs(Ensemble_BM1991[3:30] - MP_Species[c(1:25, 29:31)]) + abs(Ensemble_BM1991[3:30] - Nor_Species[c(1:25, 29:31)]))
# SAD <- (abs(Ensemble_BM2031[3:30] - CM_Species[c(1:25, 29:31)]) + abs(Ensemble_BM2031[3:30] - EE_Species[c(1:25, 29:31)]) + abs(Ensemble_BM2031[3:30] - IN_Species[c(1:25, 29:31)])
#        + abs(Ensemble_BM2031[3:30] - MP_Species[c(1:25, 29:31)]) + abs(Ensemble_BM2031[3:30] - Nor_Species[c(1:25, 29:31)]))
# SAD <- (abs(Ensemble_BM2070[3:30] - CM_Species[c(1:25, 29:31)]) + abs(Ensemble_BM2070[3:30] - EE_Species[c(1:25, 29:31)]) + abs(Ensemble_BM2070[3:30] - IN_Species[c(1:25, 29:31)])
#         + abs(Ensemble_BM2070[3:30] - MP_Species[c(1:25, 29:31)]) + abs(Ensemble_BM2070[3:30] - Nor_Species[c(1:25, 29:31)]))

SAD$Lon <- Ensemble_DENS2070$Lon  #This is a fixed list across datasets
SAD$Lat <- Ensemble_DENS2070$Lat

setwd("C:/Users/ENSEMBLE/")
write.csv(SAD, "SAD_BM_1991.csv", row.names = FALSE)