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/") Density <- read.table("dens.out", sep="", header=T) Density$Year <- as.POSIXct(paste(Density$Year), format = "%Y") Density <- as.data.frame(Density) All_yr <- gall16 %>% group_by(year=floor_date(Year, "year")) Per_yr <- All_yr %>% mutate(year=year(Year)) %>% group_by(Year, year) ########################## TIMESLICES ########################################## # # This section is optional if you wish to calculate stand density dominance # for a specific period mean; adjust taxa according to your list of species. # ################################################################################ # i.e. 1991 - 2020 DENS1 <- Per_yr %>% filter(year >= "1991" & year <= "2020") Mean_Ab <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Ab_mean = mean(Abi_nor)) Mean_Ced <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Ced_mean = mean(Ced_lib)) Mean_BES <- DENS1 %>% group_by(Lon, Lat) %>% summarise(BES_mean = mean(BES)) Mean_BP <- DENS1 %>% group_by(Lon, Lat) %>% summarise(BP_mean = mean(Bet_pen)) Mean_BPU <- DENS1 %>% group_by(Lon, Lat) %>% summarise(BPU_mean = mean(Bet_pub)) Mean_Car <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Car_mean = mean(Car_bet)) Mean_Cor <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Cor_mean = mean(Cor_ave)) Mean_Fag <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Fag_mean = mean(Fag_syl)) Mean_Fra <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Fra_mean = mean(Fra_exc)) Mean_Jun <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Jun_mean = mean(Jun_oxy)) Mean_MRS <- DENS1 %>% group_by(Lon, Lat) %>% summarise(MRS_mean = mean(MRS)) Mean_PA <- DENS1 %>% group_by(Lon, Lat) %>% summarise(PA = mean(Pic_abi)) Mean_AG <- DENS1 %>% group_by(Lon, Lat) %>% summarise(PA = mean(Aln_glu)) Mean_PS <- DENS1 %>% group_by(Lon, Lat) %>% summarise(PS_mean = mean(Pin_syl)) Mean_PN <- DENS1 %>% group_by(Lon, Lat) %>% summarise(PN_mean = mean(Pin_nig)) Mean_PB <- DENS1 %>% group_by(Lon, Lat) %>% summarise(PB_mean = mean(Pin_bru)) Mean_PH <- DENS1 %>% group_by(Lon, Lat) %>% summarise(PH_mean = mean(Pin_hal)) Mean_Pop <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Pop_mean = mean(Pop_tre)) Mean_QC <- DENS1 %>% group_by(Lon, Lat) %>% summarise(QC_mean = mean(Que_coc)) Mean_QI <- DENS1 %>% group_by(Lon, Lat) %>% summarise(QI_mean = mean(Que_ile)) Mean_QP <- DENS1 %>% group_by(Lon, Lat) %>% summarise(QP_mean = mean(Que_pub)) Mean_QR <- DENS1 %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_rob)) Mean_QM <- DENS1 %>% group_by(Lon, Lat) %>% summarise(QR_mean = mean(Que_mac)) Mean_Til <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Til_mean = mean(Til_cor)) Mean_Ulm <- DENS1 %>% group_by(Lon, Lat) %>% summarise(Ulm_mean = mean(Ulm_gla)) ###### Calculate stand density dominance ######################################### 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]) colnames(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") Species_max <- Species %>% rowwise() %>% mutate(row_max = names(.)[which.max(c_across(everything()))]) Dominant_sp <- cbind.data.frame(Mean_Ab$Lon, Mean_Ab$Lat, Species_max$row_max) colnames(Dominant_sp) <- c("Lon", "Lat", "Dominant_Sp") ###### Plot stand density dominance as a map #################################### dir.create("Dominance") setwd("C:/Users/Dominance/") borders <- sf::st_as_sf(maps::map("world", plot = FALSE, fill = TRUE)) All_sp <- Dominant_sp xmin1 = range(All_sp$Lon)[1] xmax1 = range(All_sp$Lon)[2] ymin1 = range(All_sp$Lat)[1] ymax1 = range(All_sp$Lat)[2] Dominance<- ggplot(All_sp) + geom_sf(data = borders, fill= "white")+ geom_tile(aes(x=Lon, y=Lat, fill=Dominant_Sp))+ geom_sf(data = borders,fill = NA, colour = "dark green", size=0.2)+ coord_sf(xlim = c(xmin1, xmax1), ylim = c(ymin1, ymax1))+ labs(title =paste("Dominant species 1991-2020\n"), fill="Species density")+ scale_fill_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"), limits=force)+ scale_x_longitude(ticks = 2) + scale_y_latitude(ticks = 2) + theme(legend.position="right", legend.title= element_text(size = 25, color = "dark green"), legend.text = element_text(size = 22), panel.background = element_rect(fill = "#edf0f5"), plot.title = element_text(color="dark green", size=25, face="bold"), axis.line = element_line(colour = "black"), axis.text.x=element_text(size=24), axis.text.y = element_text(size=24), panel.border = element_rect(colour = "black",fill = NA, size = 0.5)) # print(Dominance) png(file=paste("Dominance_1991-2020.png"), width=1980, height=932) print(Dominance) dev.off() ####### Calculate and plot stand density dominance per year (for the entire simulation period) ######## Species <- Per_yr[c(4:28)] colnames(Species) <- c("Abies nord.", "Cedrus lib.", "Boreal Shrub", "Betula pen.", "Betula pub.", "Carpinus b.", "Corylus a.", "Fagus s.","Fraxinus e.", "Juniperus o.", "Med. Shrub", "Picea a.", "Alnus g.", "Pinus s.", "Pinus n.", "Pinus b.", "Pinus h.", "Populus t.", "Quercus c.", "Quercus i.", "Quercus p.", "Quercus r.", "Quercus m.", "Tilia c.", "Ulmus g.") Species_max <- Species %>% rowwise() %>% mutate(row_max = names(.)[which.max(c_across(everything()))]) Dominant_sp <- cbind.data.frame(Per_yr$Lon, Per_yr$Lat, Per_yr$Year, Species_max$row_max, Per_yr$year) colnames(Dominant_sp) <- c("Lon", "Lat", "Year", "Dominant_Sp", "year") ############# BONUS: add Bioclimatic zones for further analysis ############ Bio_zones <- Dominant_sp %>% mutate(Zone = case_when(Dominant_Sp == "Quercus p." ~ "Temperate", Dominant_Sp == "Abies nord." ~ "Temperate", Dominant_Sp == "Cedrus lib." ~ "Temperate", Dominant_Sp == "Betula pen." ~ "Temperate", Dominant_Sp == "Carpinus b."~ "Temperate", Dominant_Sp == "Corylus a."~ "Temperate", Dominant_Sp == "Fagus s."~ "Temperate", Dominant_Sp == "Fraxinus e."~ "Temperate", Dominant_Sp == "Juniperus o."~ "Temperate", Dominant_Sp == "Med. Shrub"~ "Temperate", Dominant_Sp == "Pinus h." ~ "Temperate", Dominant_Sp == "Pinus n." ~ "Temperate", Dominant_Sp == "Pinus b." ~ "Temperate", Dominant_Sp == "Populus t."~ "Temperate", Dominant_Sp == "Alnus g."~ "Temperate", Dominant_Sp == "Quercus m."~ "Temperate", Dominant_Sp == "Quercus c."~ "Temperate", Dominant_Sp == "Quercus i."~ "Temperate", Dominant_Sp == "Quercus r."~ "Temperate", Dominant_Sp == "Tilia c."~ "Temperate", Dominant_Sp == "Ulmus g." ~ "Temperate", Dominant_Sp == "Boreal Shrub" ~ "Boreal", Dominant_Sp == "Betula pub." ~"Boreal", Dominant_Sp == "Pinus s." ~"Boreal", Dominant_Sp == "Picea a."~"Boreal")) ################################################################################ setwd("C:/Users/Dominance/") dir.create("Dominance_per_year") setwd("C:/Users/Dominance/Dominance_per_year") Years <- c(1961:2100) for(i in Years){ All_sp <- Dominant_sp %>% filter(year == i) xmin1 = range(All_sp$Lon)[1] xmax1 = range(All_sp$Lon)[2] ymin1 = range(All_sp$Lat)[1] ymax1 = range(All_sp$Lat)[2] #ploting the trends Dominance<- ggplot(All_sp) + geom_sf(data = borders, fill= "white")+ geom_tile(aes(x=Lon, y=Lat, fill=Dominant_Sp))+ geom_sf(data = borders,fill = NA, colour = "dark green", size=0.2)+ geom_polygon(data = shp, aes(x = long, y = lat, group = group), colour = "dark green", fill = "#edf0f5")+ geom_polygon(data = shp2, aes(x = long, y = lat, group = group), colour = "dark green", fill = "#edf0f5")+ coord_sf(xlim = c(xmin1, xmax1), ylim = c(ymin1, ymax1))+ # coord_sf(xlim = c(25.7, 42), ylim = c(40.5, 41))+ labs(title =paste("Dominant species", i), fill="Species density")+ scale_fill_manual(values = c("Quercus p." = "#e8940c", "Quercus r."="#DB6C05","Quercus i."="#57c27e", "Fagus s."="#134d19", "Fraxinus e." = "#b595e6", "Betula pub."="steel blue", "Carpinus b."="#90eb10","Betula pen."="red", "Populus t."= "#b7eb34", "Pinus h."="#f772b0", "Pinus s."="#8d49cc", "Quercus c."= "#e8550c", "Boreal Shrub" = "blue", "Med. Shrub" = "#122da6", "Abies nord."= "#799c11", "Abies cil."= "#89dc69", "Ulmus g." = "#611975", "Tilia c."= "#ad1320", "Picea a."="#e02b73", "Alnus g." = "#F87F65", "Juniperus o." = "#207852", "Corylus a." = "#0cb8e8", "Pinus n." = "#1c83ba", "Pinus b." = "grey", "Cedrus lib." = "#b7c50e", "Quercus m." = "#fcba03"), limits=force)+ scale_x_longitude(ticks = 2) + scale_y_latitude(ticks = 2) + theme(legend.position="right", legend.title= element_text(size = 25, color = "dark green"), legend.text = element_text(size = 22), panel.background = element_rect(fill = "#edf0f5"), plot.title = element_text(color="dark green", size=25, face="bold"), axis.line = element_line(colour = "black"), axis.text.x=element_text(size=24), axis.text.y = element_text(size=24), panel.border = element_rect(colour = "black",fill = NA, size = 0.5)) # print(Dominance) png(file=paste(i,"Dominance.png"), width=1980, height=932) print(Dominance) dev.off() } # Plot Zones setwd("C:/Users/ITU/Documents/DATA/GCMs/LPJ_GUESS runs/CMCC-ESM/CM_local_Ced/") dir.create("Zones") setwd("C:/Users/ITU/Documents/DATA/GCMs/LPJ_GUESS runs/CMCC-ESM/CM_local_Ced/Zones/") for(i in Years){ All_sp2 <- Bio_zones %>% filter(year == i) xmin1 = range(All_sp2$Lon)[1] xmax1 = range(All_sp2$Lon)[2] ymin1 = range(All_sp2$Lat)[1] ymax1 = range(All_sp2$Lat)[2] #ploting the trends Zones <- ggplot(All_sp2) + geom_sf(data = borders, fill= "white")+ geom_tile(aes(x=Lon, y=Lat, fill=Zone), alpha=0.6)+ geom_sf(data = borders,fill = NA, colour = "dark green", size=0.2)+ geom_polygon(data = shp, aes(x = long, y = lat, group = group), colour = "dark green", fill = "#edf0f5")+ geom_polygon(data = shp2, aes(x = long, y = lat, group = group), colour = "dark green", fill = "#edf0f5")+ coord_sf(xlim = c(xmin1, xmax1), ylim = c(ymin1, ymax1))+ labs(title =paste("Climate zones", i), fill="zones")+ scale_fill_manual(values = c("Boreal" = "steel blue", "Temperate"="pink"))+ scale_x_longitude(ticks = 2) + scale_y_latitude(ticks = 2) + theme(legend.position="right", legend.title= element_text(size = 25, color = "dark green"), legend.text = element_text(size = 22), panel.background = element_rect(fill = "#edf0f5"), plot.title = element_text(color="dark green", size=25, face="bold"), axis.line = element_line(colour = "black"), axis.text.x=element_text(size=24), axis.text.y = element_text(size=24), panel.border = element_rect(colour = "black",fill = NA, size = 0.5)) # print(Dominance) png(file=paste(i,"Zones_Cedrus1.png"), width=1980, height=932) print(Zones) dev.off() }