LPJ-GUESS-OutputPlots / Biomass_singleGCM.R
Biomass_singleGCM.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)

##############################################################################
#
# Biomass calculated for individual GCM runs
#
##############################################################################

setwd("C:/Users/")         #change direcory per GCM

gall15 <- read.table("cmass.out", sep="", header=T)
gall15$Year <- as.POSIXct(paste(gall15$Year), format = "%Y")
gall15 <- as.data.frame(gall15)
All_yr <- gall15 %>% group_by(year=floor_date(Year, "year"))
Per_yr <- All_yr %>% mutate(year=year(Year)) %>%
  group_by(Year, year)

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

########################## TIMESLICES ##########################################
# just un-commentout the respective period:

BM1961_1990 <- Per_yrBM %>% filter(year <= "1990")
BM2 <- BM1961_1990
# BM1991_2020 <- Per_yr %>% filter(year >= "1991" & year <= "2020")
# BM2 <- BM1991_2020
# BM2031_2060 <- Per_yr %>% filter(year >= "2031" & year <= "2060")
# BM2 <- BM2031_2060
# BM2071_2100 <- Per_yr %>% filter(year >= "2070")
# BM2 <- BM2071_2100


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


BM2mean <- cbind.data.frame(Mean_Ab, 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], Mean_Tot[3])
BM2_sp <- BM2mean[c(3:28)] #taking only the CMASS fields
colnames(BM2_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")


BM2_spm <- reshape2::melt(BM2_sp) #reshaping the dataframe
theme_set(theme_minimal())

BM_box <- ggplot(BM2_spm, aes(x = variable, y = value, color = variable)) +           
  geom_boxplot()+
  labs(title =paste("1961-1990 mean biomass per species"))+
  # labs(title =paste("1991-2020 mean biomass per species"))+
  # labs(title =paste("2031-2060 mean biomass per species"))+
  # labs(title =paste("2071-2100 mean biomass per species"))+
  
  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))+
  theme(legend.position="none", 
        plot.title = element_text(size=25, face="bold"),
        axis.text.x=element_text(size=24, angle = -60, hjust = 0, vjust = 0), axis.text.y = element_text(size=24))

# print(BM_box)

# png(file=paste("Biomass_1991-2020.png"), width=1980, height=932)
png(file=paste("Biomass_1961-1990.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()