################################################################################################ ### Title: Functional susceptibility of tropical forests to climate change ### Journal: Nature Ecology and Evolution March, 2022 ### Author: Jesus Aguirre Gutierrez et al. ### Institution: Environmental Change Institute, University of Oxford ### Summary of the code and functions used in the development of the manuscript specified above ### the original data is not provided and may be accessed following the Data Availability statement #################################################################################################### library(readr);library(raster);library(tidyverse); library(reshape2);library(ggridges);library(conflicted);library(data.table);library(Weighted.Desc.Stat);library(corrplot) library(ggcorrplot);library(ggrepel);library(ggiraphExtra);library(performance);library(ape);library(sf) library(mapview);library(gstat);library(ncf);library(spdep);library(stats4); library("factoextra") library(lme4);library(visreg);library(DescTools);library(visreg);library(corrplot);library(glmulti) library(tidyverse); library(reshape2);library(factoextra); library(ggpubr); library(directlabels) library(rstanarm); library(loo);library('mgcv'); library('brms'); library('schoenberg');library(bayestestR) library(brmstools);library(bayestestR);library(broom); library(modelr);library(viridis);library(tidybayes) library(RColorBrewer); library(gtools);library("bayesplot");library(rasterVis);library(gridExtra) library(utils);library(gdalUtils);library(rworldmap);library("rnaturalearth");library("rnaturalearthdata") library(ggplot2);library(elevatr) options(mc.cores = 7) conflict_prefer("select", "dplyr") conflict_prefer("mutate", "dplyr") conflict_prefer("rename", "dplyr") conflict_prefer("melt", "reshape2") conflict_prefer("filter", "dplyr") conflict_prefer("levelplot", "rasterVis") hsTheme <- modifyList(GrTheme(), list(regions=list(alpha=1))) colr <- colorRampPalette((brewer.pal(11, 'RdBu'))) colr1<- colorRampPalette((brewer.pal(11, 'YlGnBu'))) colr2<- colorRampPalette((brewer.pal(11, 'BrBG'))) col3 <- viridis col4<- colorRampPalette((brewer.pal(11, 'Spectral'))) is_outlier <- function(x) { return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x)) } setwd("/Users/your environment") ############################# #### Get Environmental data # ############################# si<-read_csv('PlotData.csv')%>%unique() #Plot code data env<-read_csv('ClimateData.csv')%>% #Environmental data group_by(Loc,Plot)%>%summarise_all(.funs=mean)%>%merge(si) ################################################# # Get diversity data and merge with environment # ################################################# div<-read_csv('DiversityTable.csv')%>%unique() #Functional Diversity and Redundancy data div_env<-merge(div, env, by="Plot") ################################################# #Moran I test spatial autocorrelation ########## ################################################# cen<-read_csv('PlotsCentroids.csv')%>%unique() #Centroids of plots Fd<-div_env%>%merge(cen[,c(1,3,4)], by='Plot') Fd.dists <- as.matrix(dist(cbind(Fd$x, Fd$y))) Fd.dists.inv <- 1/Fd.dists diag(Fd.dists.inv) <- 0 Fd.dists.inv[1:5, 1:5] Moran.I(Fd$Value, Fd.dists.inv) ######################################################################### #Plot locations and creating column of Plot Group based on the Spatial #Autocorrelation observed from above. The spatial autocorrelation #decreases at 2 km ######################################################################### cen_s <- st_as_sf(x = cen, coords = c("x", "y"),crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") mdist <- dist(cen[c("x","y")]) hc <- hclust(mdist, method="complete") d=2 ##distance of ~2km cen$Clust_2km = cutree(hc, h=d) div_env<- merge(div_env,cen[,c(1,5)])%>%mutate(Clust_2km=as.factor(Clust_2km)) ################################################## # Soil data from SOILGRIDS #Tuparro falls into water in the SoilGrids # include it in the main data frame #So I took the values of the surrounding cells ##################################################for all soil variables---- sg<- read_csv('SoilGrids.csv') div_env<-div_env%>%merge(sg, by=c('Plot','Loc'), all.x=T, all.y=F) ############################################################ # DATA FOR PANTROPIAL PREDICTIONS AFTER STATITSICAL MODELS # extracted from Google Earth Engine ############################################################ for_crs<-'+proj=longlat +datum=WGS84 +no_defs' pantro_all<-read_csv('mcwdvpd_results_TABLE.csv')%>%select(1,2,5,6,12,13) names(pantro_all)<-c('x','y', 'FT_VPD','Ch_VPD','FT_MCWD','Ch_MCWD') pantro_all1<-pantro_all%>%na.omit()%>%mutate(Size=1,Clust_2km='No') onlymask<-rasterFromXYZ(pantro_all1[1:3])#to mask soil and climate below crs(onlymask)<-'+proj=longlat +datum=WGS84 +no_defs' treecover<-raster('TreeCover_Hansen.tif')>30 treecover<-projectRaster(treecover,onlymask, method='ngb')%>%mask(onlymask) treecover[treecover<1]<-NA pantro_extract<-pantro_all1 coordinates(pantro_extract)<-~x+y clay_pantropical<-raster('Clay0_30cm_soilgrids_pantropical.tif') clay_pantropical<-projectRaster(clay_pantropical, treecover, method="bilinear") names(clay_pantropical)<-'clay_0to30cm' clay_pantropical<-clamp(clay_pantropical, lower=8, upper=51 , useValues=TRUE) cec_pantropical1<-raster('CEC0_30cm_soilgrids_pantropical.tif') cec_pantropical1<-projectRaster(cec_pantropical1, treecover, method="bilinear") names(cec_pantropical1)<-'cec_0to30cm' cec_pantropical1<-clamp(cec_pantropical1, upper=34 , useValues=TRUE) cec_clay<-stack(clay_pantropical,cec_pantropical1) real_ftmcwd<-stack('MCWDfullterm5817.tif') real_ftmcwd<-projectRaster(real_ftmcwd[[4]]*-1, treecover)#%>%mask(treecover) real_ftmcwd<-clamp(real_ftmcwd, upper=941.6207, useValues=TRUE)#1000 names(real_ftmcwd)<-'FT_MCWD' #VPD vpd_real<-stack('VPD_realTerraClimate2021.tif') names(vpd_real)<-c("vpd_5887","vpd_8817","FT_VPD","Ch_VPD") vpd_real<-projectRaster(vpd_real, treecover, method="bilinear")#%>%mask(treecover) #stack soil and new MCWD cec_clay_mcwd1<-stack(cec_clay,real_ftmcwd,real_chmcwd,vpd_real[[3:4]])%>%mask(treecover) plot(cec_clay_mcwd1) #Extract soil and new MCWD to table for the tropics extracted<-raster::extract(cec_clay_mcwd1, pantro_extract) pantro_all1<-cbind(pantro_all1, extracted) pantro_all1<-pantro_all1%>% select(1,2,13,14,11,12,7:10)%>%na.omit() ############################################# # INCLUDE TRUE MCWD AND DELTA MCWD from TerraClimate Monthly data ################## cen2<-cen[,c(3,4,1)] coordinates(cen2)<-~x+y cen3<-raster::extract(stack(clay_pantropical,cec_pantropical1),cen2)%>%cbind(as.data.frame(cen2)) cen4<-raster::extract(stack(real_ftmcwd,real_chmcwd,vpd_real[[3:4]]),cen2)%>%cbind(as.data.frame(cen2)) xxx<-div_env%>%select(1,2)%>%unique()%>%merge(cen3[,c(1,2,5)], by='Plot')%>% merge(cen4[,c(1:4,7)], by='Plot')%>%rename(clay_0to30cm=Clay0_30cm_soilgrids_pantropical, cec_0to30cm=CEC0_30cm_soilgrids_pantropical) div_env<-div_env%>%merge(xxx[,c(1,3:8)],by='Plot',all.x=T, all.y=F) div_env<-div_env%>%select(1:23,47,25:27,45,29:32,44,46,35:37,43,42,40,41) div_env<-div_env%>%rename(FT_MCWD=FT_MCWD.y,Ch_MCWD=Ch_MCWD.y,Ch_VPD=Ch_VPD.y,FT_VPD=FT_VPD.y,cec_0to30cm=cec_0to30cm.y, clay_0to30cm=clay_0to30cm.y) div_env%>%select(1,2,24,25,28,32,33:35,38,39)%>%unique()%>% write_csv('TableClimatePlots.csv') ####################################################### # MAKE PLOTS OF ENVIRONMENTAL VARIABLES FOR SI ######## ####################################################### cec_clay_mcwd2<-mask(cec_clay_mcwd1,treecover) names(cec_clay_mcwd2)<-c('Clay %)', 'CEC (mmol c/Kg)', 'MCWD (mm)', expression(paste(Delta,'MCWD (mm)')), 'VPD (kPa)', expression(paste(Delta,'VPD (kPa)'))) si_maps_climate<- list() for(i in 1:6){ figname=names(cec_clay_mcwd2[[i]]) si_maps_climate[[i]]<-levelplot(cec_clay_mcwd2[[i]],maxpixels = 2e5, margin=FALSE, colorkey=list( #space='bottom', axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = figname) + latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))} do.call(grid.arrange,c(si_maps_climate, ncol=2)) ###################################### # CORRELATION analysis of covariates # # for SI material # ###################################### c<-div_env%>%select(1,24,25,28,33:35, 38:41)%>%unique()#%>%select(1,8:10) names(c)<-c('Plot','ΔVPD','ΔCV','ΔMCWD','MCWD','VPD','CV', 'CEC', 'Clay', 'pH', 'Sand') cc = cor(c[,2:11], method = "spearman") ggcorrplot(cc, hc.order = TRUE, type = "lower",lab = TRUE) ggsave("Correlation_covariates.png", width = 20, height = 12, units = "cm") ####################################### # GRAPGH where the plots fall in # # environmental space # ####################################### env_pca<- div_env%>%select(1,2,24,25,28,33:35)%>%unique()%>% arrange(Loc)%>%select(1,3:8)%>%column_to_rownames('Plot') names(env_pca)<-c('ΔVPD','ΔCV','ΔMCWD','MCWD','VPD','CV') pca1<-prcomp(env_pca[,c(1,3,4:5)], scale=T) summary(pca1) plot(pca1) group <- c(rep("Australia", times=7),rep("Colombia", times=11), rep("Gabon", times=3), rep("Ghana", times=9),rep("Malaysia", times=4),rep("Mexico", times=14), rep("NovaX", times=4),rep("Panama", times=12), rep("Peru", times=10), rep("Santarem", times=5)) fviz_pca_biplot(pca1, repel=T, pointsize=8, pointshape=20, col.var="black", arrowsize=0.9, labelsize=4,label = "var", col.ind = group, palette=getPalette, addEllipses=T, ellipse.type="confidence",#ellipse.level=.95, gradient.cols = group, max.overlaps=30, alpha.ind=0.7, legend.title = "Study areas")+theme_pubclean()#+theme_ridges() ggsave("PCA_climate_and_plots.png", width = 20, height = 14, units = "cm") pca2<-prcomp(env_pca[,c(4:5)], scale=T) fviz_pca_biplot(pca2, repel=F, pointsize=2, pointshape=20, col.var="black", arrowsize=0.6, labelsize=3, col.ind=group, palette=getPalette, addEllipses=TRUE, ellipse.type="confidence")+theme_ridges() ggsave("PCA_climateMCWDVPD_and_plots.png", width = 20, height = 14, units = "cm") #PCA FOR SOIL env_pca3<- div_env%>%select(1,2,38:41)%>%unique()%>% arrange(Loc)%>%select(1,3:6)%>%column_to_rownames('Plot') names(env_pca2)<-c('CEC','Clay','pH','Sand') pca5<-prcomp(env_pca3, scale=T) summary(pca5) plot(pca5) group2 <- c(rep("Australia", times=7),rep("Colombia", times=8), rep("Gabon", times=3), rep("Ghana", times=9),rep("Malaysia", times=4),rep("Mexico", times=13), rep("NovaX", times=4),rep("Panama", times=11), rep("Peru", times=10), rep("Santarem", times=5)) fviz_pca_biplot(pca5, repel=TRUE, pointsize=4, pointshape=20, col.var="black", arrowsize=0.9, labelsize=4,label = "var", col.ind = group2, palette=getPalette, addEllipses=T, ellipse.type="confidence",#ellipse.level=.95, gradient.cols = group, max.overlaps=30, alpha.ind=0.3, legend.title = "Study areas")#+theme_pubclean()#+theme_ridges() ggsave("Soil_PCA_plots.png", width = 20, height = 14, units = "cm") ################################### #Select colour palettes for lines ################################### myColorScale_1 <- brewer.pal(11,"RdBu");myColorScale_1<-myColorScale_1[c(7:11)] myColorScale_2 <- brewer.pal(11,"Spectral");myColorScale_2<-myColorScale_2[c(1:5)] myColorScale_2 <- c(myColorScale_1[c(1:5)],'gold4',myColorScale_2[c(1:5)],'orange2','orange3') levels_traits<-div_env%>%select(5,36)%>%unique()%>%arrange(traittype,Trait)%>%slice(c(1:12,25:42,44:47))%>%select(1,2) levels_traits<-levels_traits%>%mutate(Trait=as.factor(Trait)) names(myColorScale_2)<- fct_reorder(levels_traits$Trait, levels_traits$traittype) colScale <- scale_colour_manual(name = "Trait",values = myColorScale_2) ################################################################### # COMPARE outcome of Models with and without RANDOM FACTOR ######## ################################################################### # Example: Fd<-div_env%>%subset(Metric=='FD'&is.na(Value)==F)%>%filter(str_detect(Trait, "\\_BA"))%>% subset(Trait=="FDis_mor_BA"|Trait=="FDis_nut_BA")%>%dplyr::select(1:5,24,25,28,32:39)%>% mutate_at(.vars =c(6:12,15,16), .funs = list("scaled" = scale)) Fd <-Fd%>%mutate_at(vars(names(Fd[17:25])), ~as.numeric(as.character(.))) trait_mod<-'FDis_mor_BA' Fd_mod<-Fd%>%subset(Trait==trait_mod) F1<- stan_glm(Value ~ Ch_MCWD_scaled*Ch_VPD_scaled + Size_scaled+FT_MCWD_scaled*FT_VPD_scaled+ cec_0to30cm_scaled+ clay_0to30cm_scaled, data=Fd_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=2000, control = list(adapt_delta = 0.99)) F1R<-stan_glmer(Value ~ Ch_MCWD_scaled*Ch_VPD_scaled + Size_scaled+FT_MCWD_scaled*FT_VPD_scaled + cec_0to30cm_scaled+ clay_0to30cm_scaled + (1|Clust_2km), data=Fd_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99)) draws=500 F1_loo <-loo(F1, cores = 12, k_threshold = 0.7) F1R_loo <-loo(F1R, cores = 12, k_threshold = 0.7) F1_best <-loo_compare(F1_loo,F1R_loo) F1_best ########################################################### # FUNCTIONAL DIVERSITY as a function of MCWD, VPD and SOIL # I do not do interactions for the Photosynthesis models # because the lower number of records ########################################################## Fd<-div_env%>%subset(Metric=='FD'&is.na(Value)==F)%>%filter(str_detect(Trait, "\\_BA"))%>% subset(Trait=="FDis_mor_BA"|Trait=="FDis_nut_BA"|Trait=='FDis_pho_BA')%>% dplyr::select(1:5,24,25,28,32:39)%>%subset(Value!=0) ##### Check spatial autocorrelation FD models ##### -Repeat the same for FDis of Morphological, ##### Nutrients and Photosynthesis traits and for the same of FRed draws=700 trial_fd<-Fd trait_mod<-unique(trial_fd$Trait)[1] FD_mod<-trial_fd%>%subset(Trait==trait_mod) hist(FD_mod$Value) outliers <- boxplot(FD_mod$Value, plot=FALSE)$out FD_mod <- FD_mod[-which(FD_mod$Value %in% outliers),] if (nrow(FD_mod)>0){ print('All good') }else { FD_mod <- trial_fd%>%subset(Trait==trait_mod) } set.seed(103) scaled_FD <-scale(FD_mod[,c(6,8,9,10,11, 15, 16)]) FD_mod <- FD_mod%>%select(1:5,13,14)%>%cbind(scaled_FD) phoFD_Model_Climate<- stan_glmer(log(Value) ~ FT_MCWD+FT_VPD+ Ch_MCWD+Ch_VPD + Size + (1|Clust_2km), data=FD_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99)) phoNoRF_FD_Model_Climate<- stan_glm(log(Value) ~ FT_MCWD+FT_VPD+ Ch_MCWD+Ch_VPD + Size , data=FD_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99)) plot(phoFD_Model_Climate$residuals,phoFD_Model_Climate$fitted.values,main='phoFD_Model_Climate') plot(phoNoRF_FD_Model_Climate$residuals,phoNoRF_FD_Model_Climate$fitted.values,main='phoNoRF_FD_Model_Climate') co<-read_csv('PlotsLocations.csv') phoNRF_tri<-merge(FD_mod,co, by=c('Plot','Loc'),all.x=T, all.y=F) dists <- as.matrix(dist(phoNRF_tri[,15:16])) dists.inv <- 1/dists diag(dists.inv) <- 0 phoNRF_MI<-Moran.I(phoNoRF_FD_Model_Climate$residuals, dists.inv)%>%as.data.frame()%>%mutate(Index='phoNRF_MI') phoRF_MI <-Moran.I(phoFD_Model_Climate$residuals, dists.inv)%>%as.data.frame()%>%mutate(Index='phoRF_MI') #After carrying out the above for FDis of Morphological, Nutrients and Photosynthesis traits and for the same of FRed #Put all results above in dataframe rbind(phoNRF_MI,phoRF_MI,nutNRF_MI,nutRF_MI,morNRF_MI,morRF_MI,phoNRF_fr_MI,phoRF_fr_MI,nutNRF_fr_MI,nutRF_fr_MI,morNRF_fr_MI,morRF_fr_MI)%>% write_csv('MoranI_results.csv') #Check semivariogram to see distance for autocorrelation trialsac<-cbind(nutNRF_fr_models_trial$residuals,nutNRF_tri)%>%rename(res=`nutNRF_fr_models_trial$residuals`)%>%st_as_sf(coords = c("Long", "Lat"),crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") sacnut_var<-variogram(res~1, trialsac, width=.1) plot(sacnut_var,pch=20,cex=1.5,col="black", ylab=expression("Semivariance ("*gamma*")"), xlab="Distance (m)",main = "FD nutrients Moran's I") ############################# # Running the models for FD ############################# trial_fd<-Fd Models_table<-list() ForPlotting<-list() scaleListFD<-list() fd_mapspred_done<-list() pdtfigs_fd<-list() ndraws = 700 for (i in 1:length(unique(trial_fd$Trait))){ trait_mod<-unique(trial_fd$Trait)[i] FD_mod<-trial_fd%>%subset(Trait==trait_mod) hist(FD_mod$Value) if (nrow(FD_mod)>0){ print('All good') }else { FD_mod <- trial_fd%>%subset(Trait==trait_mod) } set.seed(103) scaled_FD <-scale(FD_mod[,c(6,8,9,10,11, 15, 16)]) scaleListFD[[i]]<- data.frame(name=trait_mod, scale = attr(scaled_FD, "scaled:scale"), center = attr(scaled_FD, "scaled:center"))%>%rownames_to_column()%>%rename(varia=rowname) FD_mod <- FD_mod%>%select(1:5,13,14)%>%cbind(scaled_FD) if (trait_mod=='FDis_pho_BA'){ FD_models_trial<- stan_glmer(log(Value) ~ FT_MCWD+FT_VPD+ Ch_MCWD+Ch_VPD+ Size + (1|Clust_2km), data=FD_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=10000, control = list(adapt_delta = 0.99)) }else { if (trait_mod=='FDis_nut_BA'){ FD_models_trial<- stan_glmer(log(Value) ~ FT_MCWD*FT_VPD + Ch_MCWD*Ch_VPD + Size + cec_0to30cm + clay_0to30cm + (1|Clust_2km), data=FD_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=10000, control = list(adapt_delta = 0.99)) }else{ FD_models_trial<- stan_glm(log(Value) ~ FT_MCWD*FT_VPD+ Ch_MCWD*Ch_VPD + Size + cec_0to30cm + clay_0to30cm, data=FD_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=10000, control = list(adapt_delta = 0.99)) } } pdp<-as.matrix(FD_models_trial)%>%mcmc_areas(prob = 0.9)+ggtitle(trait_mod,"Posterior distributions with medians and 90% intervals") pdtfigs_fd[[i]]<-list(pdp) FD_trial_posterior<-describe_posterior(FD_models_trial,rope_ci=.9, ci=.9)%>%as.data.frame()%>% mutate(Trait=trait_mod, R2=bayes_R2(FD_models_trial) %>% median(), R2cond=as.numeric(r2_bayes(FD_models_trial,ci = 0.90)[1]), R2AdjLoo=as.numeric(r2_loo(FD_models_trial)))%>% mutate(NPlots=nrow(FD_mod)) Models_table[[i]]<-FD_trial_posterior FD_models_trial_draws1<-FD_mod %>% data_grid(FT_MCWD =seq_range(FT_MCWD, n = 101),cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_VPD = mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='FT_MCWD') FD_models_trial_draws2<-FD_mod %>% data_grid(FT_VPD =seq_range(FT_VPD, n = 101), cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_MCWD = mean(FT_MCWD),Ch_VPD=mean(Ch_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='FT_VPD') FD_models_trial_draws3<-FD_mod %>% data_grid(Ch_VPD =seq_range(Ch_VPD, n = 101), cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='Ch_VPD') FD_models_trial_draws4<-FD_mod %>% data_grid(Ch_MCWD =seq_range(Ch_MCWD, n = 101),cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='Ch_MCWD') FD_models_trial_draws5<-FD_mod %>% data_grid(cec_0to30cm=seq_range(cec_0to30cm, n = 101), Ch_MCWD = mean(Ch_MCWD), clay_0to30cm = mean(clay_0to30cm), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='cec_0to30cm') FD_models_trial_draws6<-FD_mod %>% data_grid(clay_0to30cm=seq_range(clay_0to30cm, n = 101),cec_0to30cm= mean(cec_0to30cm), Ch_MCWD = mean(Ch_MCWD), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='clay_0to30cm') FD_models_trial_draws7<-FD_mod %>% group_by(FT_VPD=ifelse(FT_VPD>=mean(FT_VPD), max(FT_VPD),min(FT_VPD)))%>% data_grid(FT_MCWD = seq_range(FT_MCWD, n = 101), Ch_VPD=mean(Ch_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size), cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='FT_MCWD:FT_VPD') FD_models_trial_draws8<-FD_mod %>% group_by(Ch_VPD=ifelse(Ch_VPD>=mean(Ch_VPD), max(Ch_VPD),min(Ch_VPD)))%>% data_grid(Ch_MCWD = seq_range(Ch_MCWD, n = 101), FT_VPD=mean(FT_VPD), FT_MCWD=mean(FT_MCWD), Size=mean(Size), cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), Clust_2km='Nowhere') %>% add_fitted_draws(FD_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='Ch_MCWD:Ch_VPD') FD_models_trial_draws<-bind_rows(FD_models_trial_draws1,FD_models_trial_draws2,FD_models_trial_draws3,FD_models_trial_draws4,FD_models_trial_draws5,FD_models_trial_draws6,FD_models_trial_draws7,FD_models_trial_draws8) ForPlotting[[i]]<-FD_models_trial_draws print(paste('yes! number',i,'out of', length(unique(trial_fd$Trait)))) ################################## ########## PREDICT TO FULL TROPICS ################################## pantro_all_FD<-pantro_all1%>%mutate(Ch_VPD= scale(Ch_VPD, center = attr(scaled_FD, "scaled:center")[1],scale =attr(scaled_FD, "scaled:scale")[1]), Ch_MCWD=scale(Ch_MCWD,center = attr(scaled_FD, "scaled:center")[2],scale =attr(scaled_FD, "scaled:scale")[2]), Size= scale(Size, center = attr(scaled_FD, "scaled:center")[3],scale =attr(scaled_FD, "scaled:scale")[3]), FT_MCWD=scale(FT_MCWD,center = attr(scaled_FD, "scaled:center")[4],scale =attr(scaled_FD, "scaled:scale")[4]), FT_VPD= scale(FT_VPD, center = attr(scaled_FD, "scaled:center")[5],scale =attr(scaled_FD, "scaled:scale")[5]), cec_0to30cm=scale(cec_0to30cm,center = attr(scaled_FD, "scaled:center")[6],scale =attr(scaled_FD, "scaled:scale")[6]), clay_0to30cm=scale(clay_0to30cm,center = attr(scaled_FD, "scaled:center")[7],scale =attr(scaled_FD, "scaled:scale")[7]))%>%na.omit() names(pantro_all_FD)<-c("x","y","FT_VPD", "Ch_VPD", "FT_MCWD", "Ch_MCWD", "Size", "Clust_2km","clay_0to30cm", "cec_0to30cm") fd_mapspred<-posterior_predict(FD_models_trial, pantro_all_FD, draws = 700) fd_mapspred1<-as.data.frame(fd_mapspred) fd_mapspred2<-fd_mapspred1%>%colMeans()%>%melt()%>%cbind(pantro_all_FD)%>%rename(prediction=value) fd_mapspred3<-fd_mapspred2%>%select(2,3,1)%>%mutate(prediction=exp(prediction))%>%rasterFromXYZ() names(fd_mapspred3)<-trait_mod plot(fd_mapspred3) fd_mapspred_done[[i]]<-fd_mapspred3 } ############################################################################################################## # FUNCTIONAL REDUNDANCY as a function of continuous MCWD and VPD ############################################# # I do not do interactions for the Photosynthesis models because the low number of records ################### ######continuous models ###################################################################################### Fr<-div_env%>%subset(Metric=='FD'&is.na(Value)==F)%>% dplyr::select(1:5,24,25,28,32:39)%>%subset(Value!=0) trial_fr<-Fr Models_table<-list() ForPlotting<-list() scaleListFRed<-list() fred_mapspred_done<-list() pdtfigs<-list() ndraws = 700 for (i in 1:length(unique(trial_fr$Trait))){ trait_mod<-unique(trial_fr$Trait)[i] fr_mod<-trial_fr%>%subset(Trait==trait_mod) hist(fr_mod$Value) if (nrow(fr_mod)>0){ print('All good') }else { fr_mod <- trial_fr%>%subset(Trait==trait_mod) } set.seed(104) scaled_Fred <-scale(fr_mod[,c(6,8,9,10,11,15,16)]) scaleListFRed[[i]]<- data.frame(name=trait_mod, scale = attr(scaled_Fred, "scaled:scale"), center = attr(scaled_Fred, "scaled:center"))%>%rownames_to_column()%>%rename(varia=rowname) fr_mod <- fr_mod%>%select(1:5,13,14)%>%cbind(scaled_Fred) if (trait_mod=='RStar_pho_BA'){ fr_models_trial<- stan_glm(log(Value) ~ FT_MCWD + FT_VPD + Ch_MCWD + Ch_VPD + Size, data=fr_mod, algorithm = "sampling", cores = 13, seed = 18,chains=3, iter=10000, control = list(adapt_delta = 0.99)) }else { if (trait_mod=='RStar_nut_BA'){ fr_models_trial<- stan_glmer(log(Value) ~ FT_MCWD*FT_VPD + Ch_MCWD*Ch_VPD + Size + cec_0to30cm + clay_0to30cm + (1|Clust_2km), data=fr_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=10000, control = list(adapt_delta = 0.99)) }else{ fr_models_trial<- stan_glm(log(Value) ~ FT_MCWD*FT_VPD + Ch_MCWD*Ch_VPD + Size + cec_0to30cm + clay_0to30cm, data=fr_mod, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=10000, control = list(adapt_delta = 0.99)) } } pdp<-as.matrix(fr_models_trial)%>%mcmc_areas(prob = 0.9)+ggtitle(trait_mod,"Posterior distributions with medians and 90% intervals") pdtfigs[[i]]<-list(pdp) fr_trial_posterior<-describe_posterior(fr_models_trial,rope_ci=.9, ci=.9)%>%as.data.frame()%>%mutate(Trait=trait_mod, R2=bayes_R2(fr_models_trial) %>% median(), R2cond=as.numeric(r2_bayes(fr_models_trial,ci = 0.90)[1]), R2AdjLoo=as.numeric(r2_loo(fr_models_trial)))%>%mutate(NPlots=nrow(fr_mod)) Models_table[[i]]<-fr_trial_posterior fr_models_trial_draws1<-fr_mod %>% data_grid(FT_MCWD = seq_range(FT_MCWD, n = 101),cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_VPD = mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='FT_MCWD') fr_models_trial_draws2<-fr_mod %>% data_grid(FT_VPD = seq_range(FT_VPD, n = 101) ,cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_MCWD = mean(FT_MCWD),Ch_VPD=mean(Ch_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='FT_VPD') fr_models_trial_draws3<-fr_mod %>% data_grid(Ch_VPD = seq_range(Ch_VPD, n = 101) ,cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_MCWD=mean(Ch_MCWD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='Ch_VPD') fr_models_trial_draws4<-fr_mod %>% data_grid(Ch_MCWD = seq_range(Ch_MCWD, n = 101),cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm),FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='Ch_MCWD') fr_models_trial_draws5<-fr_mod %>% data_grid(cec_0to30cm=seq_range(cec_0to30cm, n = 101), clay_0to30cm = mean(clay_0to30cm), Ch_MCWD = mean(Ch_MCWD), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='cec_0to30cm') fr_models_trial_draws6<-fr_mod %>% data_grid(clay_0to30cm=seq_range(clay_0to30cm, n = 101),cec_0to30cm= mean(cec_0to30cm), Ch_MCWD = mean(Ch_MCWD), FT_MCWD = mean(FT_MCWD),FT_VPD=mean(FT_VPD), Ch_VPD=mean(Ch_VPD), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='clay_0to30cm') fr_models_trial_draws7<-fr_mod %>% group_by(FT_VPD=ifelse(FT_VPD>=mean(FT_VPD), max(FT_VPD),min(FT_VPD)))%>% data_grid(FT_MCWD = seq_range(FT_MCWD, n = 101), Ch_VPD=mean(Ch_VPD), Ch_MCWD=mean(Ch_MCWD), cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='FT_MCWD:FT_VPD') fr_models_trial_draws8<-fr_mod %>% group_by(Ch_VPD=ifelse(Ch_VPD>=mean(Ch_VPD), max(Ch_VPD),min(Ch_VPD)))%>% data_grid(Ch_MCWD = seq_range(Ch_MCWD, n = 101), FT_VPD=mean(FT_VPD), FT_MCWD=mean(FT_MCWD), cec_0to30cm = mean(cec_0to30cm), clay_0to30cm = mean(clay_0to30cm), Size=mean(Size),Clust_2km='Nowhere') %>% add_fitted_draws(fr_models_trial, n = ndraws, allow_new_levels = T) %>% mutate(median=median(.value), Trait=trait_mod, MappingVar='Ch_MCWD:Ch_VPD') fr_models_trial_draws<-bind_rows(fr_models_trial_draws1,fr_models_trial_draws2,fr_models_trial_draws3,fr_models_trial_draws4,fr_models_trial_draws5,fr_models_trial_draws6,fr_models_trial_draws7,fr_models_trial_draws8) ForPlotting[[i]]<-fr_models_trial_draws print(paste('yes! number',i,'out of', length(unique(trial_fr$Trait)))) ################################## ########## PREDICT TO FULL TROPICS ################################## pantro_all_FRed<-pantro_all1%>%mutate(Ch_VPD= scale(Ch_VPD, center = attr(scaled_Fred, "scaled:center")[1],scale =attr(scaled_Fred, "scaled:scale")[1]), Ch_MCWD=scale(Ch_MCWD,center = attr(scaled_Fred, "scaled:center")[2],scale =attr(scaled_Fred, "scaled:scale")[2]), Size= scale(Size, center = attr(scaled_Fred, "scaled:center")[3],scale =attr(scaled_Fred, "scaled:scale")[3]), FT_MCWD=scale(FT_MCWD,center = attr(scaled_Fred, "scaled:center")[4],scale =attr(scaled_Fred, "scaled:scale")[4]), FT_VPD= scale(FT_VPD, center = attr(scaled_Fred, "scaled:center")[5],scale =attr(scaled_Fred, "scaled:scale")[5]), cec_0to30cm=scale(cec_0to30cm,center = attr(scaled_Fred, "scaled:center")[6],scale =attr(scaled_Fred, "scaled:scale")[6]), clay_0to30cm=scale(clay_0to30cm,center = attr(scaled_Fred, "scaled:center")[7],scale =attr(scaled_Fred, "scaled:scale")[7]))%>% na.omit() names(pantro_all_FD)<-c("x","y","FT_VPD", "Ch_VPD", "FT_MCWD", "Ch_MCWD", "Size", "Clust_2km","clay_0to30cm", "cec_0to30cm") fred_mapspred<-posterior_predict(fr_models_trial, pantro_all_FRed, draws = 700) fred_mapspred1<-as.data.frame(fred_mapspred) fred_mapspred2<-fred_mapspred1%>%colMeans()%>%melt()%>%cbind(pantro_all_FRed)%>%rename(prediction=value) fred_mapspred3<-fred_mapspred2%>%select(2,3,1)%>%mutate(prediction=exp(prediction))%>%rasterFromXYZ()#use 'exp()' because in the model I log transformed the response variable names(fred_mapspred3)<-trait_mod plot(fred_mapspred3) fred_mapspred_done[[i]]<-fred_mapspred3 } ################################################################################################################ # SAVE TABLES WITH RESULTS FROM ALL FD-FRed MODELS ########################################################## ################################################################################################################ bind_rows(Models_table_all_M_FULL_FD, Models_table_all_M_FULL_fr)%>% write_csv('FDFRed_GLMM_models_tablesresults.csv') ################################################################################### # FREQUENCY DISTRIBUTION OF HIGH-MEDIUM-LOW values in FD-FRed shown in Manuscript # ################################################################################### FDFRED_mapspred_done<-stack(stack(fd_mapspred_done),stack(fred_mapspred_done)) FreqDistFDFRed<-list() for (i in 1:6){ IndexName<-names(FDFRED_mapspred_done[[i]]) ras<-FDFRED_mapspred_done[[i]] ra1<-summary(ras)[1] ra2<-summary(ras)[5] ra3<-(ra2-ra1)/3 fd_reclassified<-reclassify(ras, c(0.000001, (ra1+ra3),1,(ra1+ra3+0.00001),(ra1+ra3+ra3),2,(ra1+ra3+ra3+0.00001),Inf,3)) TotalCells <-sum(freq(fd_reclassified, value=1),freq(fd_reclassified, value=2),freq(fd_reclassified, value=3)) CellsLow <-(freq(fd_reclassified, value=1)*100)/TotalCells CellsMedium<-(freq(fd_reclassified, value=2)*100)/TotalCells CellsHigh <-(freq(fd_reclassified, value=3)*100)/TotalCells allinfo<-data.frame(IndexName, TotalCells, CellsLow, CellsMedium, CellsHigh) FreqDistFDFRed[[i]]<- allinfo } FreqDistFDFRed<-bind_rows(FreqDistFDFRed)%>% write_csv('FreqDistFDFRed_wSoils_results.csv') ####################################################################################################################### ####################################################################################################################### ############################################# FD and FRed prediction maps ############################################# ####################################################################################################################### world <- ne_countries(scale = "medium") world1 <- ne_countries(scale = "medium", returnclass = "sf") ggplot(data = world) + theme_classic()+ geom_sf() + coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ") worldmap <- getMap(resolution = "coarse") ############################ PUT ALL PLOTS OF PREDICTIONS TOGETHER---- FDFRED_mapspred_done<-stack(stack(fd_mapspred_done[[3]],fd_mapspred_done[[2]],fd_mapspred_done[[1]]), stack(fred_mapspred_done[[1]],fred_mapspred_done[[3]], fred_mapspred_done[[2]])) names(FDFRED_mapspred_done) <- c('FD -Morphology/Structure','FD -Nutrients','FD -Photosynthesis','FRed -Morphology/Structure', 'FRed -Nutrients','FRed -Photosynthesis') crs(FDFRED_mapspred_done)<- CRS("+proj=longlat +datum=WGS84") new_crs_maps<-'+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'#Molligden projection FDFRED_mapspred_done1<-FDFRED_mapspred_done%>%projectRaster(crs=new_crs_maps) worldmap <- getMap(resolution = "coarse") worldmap1<-spTransform(worldmap,crs(new_crs_maps)) elevation<- get_elev_raster(FDFRED_mapspred_done[[1]], z = 1) elevation[elevation<0]<-NA slope <- terrain(elevation, opt='slope') aspect <- terrain(elevation, opt='aspect') hill <- hillShade(slope, aspect) hill <- crop(hill, FDFRED_mapspred_done) plot(hill, col=grey(0:100/100), legend=FALSE) hsTheme <- modifyList(GrTheme(), list(regions=list(alpha=1))) colr <- colorRampPalette(#rev (brewer.pal(11, 'RdBu'))) colr1<- colorRampPalette((brewer.pal(11, 'YlGnBu'))) colr2<- colorRampPalette((brewer.pal(11, 'BrBG'))) col3 <- viridis col4<- colorRampPalette((brewer.pal(11, 'Spectral'))) ################## maps here library(scico) p <- list() for(i in 1:6){ figname=ifelse(names(FDFRED_mapspred_done)[[i]]=="FD..Photosynthesis",'FD -Photosynthesis', ifelse(names(FDFRED_mapspred_done)[[i]]=="FD..Nutrients",'FD -Nutrients', ifelse(names(FDFRED_mapspred_done)[[i]]=="FD..Morphology.Structure",'FD -Morphology/Structure', ifelse(names(FDFRED_mapspred_done)[[i]]=="FRed..Morphology.Structure",'FRed -Morphology/Structure', ifelse(names(FDFRED_mapspred_done)[[i]]=="FRed..Photosynthesis",'FRed -Photosynthesis','FRed -Nutrients'))))) rastmap<-FDFRED_mapspred_done[[i]] rastmap<-rastmap%>%as.data.frame(xy=T)%>%na.omit() names(rastmap)<-c("x","y","Val") plotx<-ggplot(rastmap, aes(x = x, y = y)) + ggtitle(figname)+ geom_sf(data=world1,fill='gray80',col='gray90',inherit.aes = FALSE)+ geom_raster(aes(fill = Val)) + scico::scale_fill_scico(palette = "roma" ,na.value = "transparent") + #bilbao,roma theme(text = element_text(size = 10, colour = "black"),panel.background = element_rect(fill = "white")) + borders(colour = "gray60", size = 0.05) + coord_sf(expand = FALSE, xlim = c(-120, 181), ylim = c(-23.5, 23.5))+ theme(legend.position = "right", legend.title = element_blank(), plot.background = element_blank(), strip.text = element_text(size = 12, colour = "black"), panel.border = element_blank(), axis.text.y = element_text(angle = 90, hjust = 0.5), axis.text = element_text(size = 9, colour = "black"), axis.title = element_text(size = 9, colour = "black"), axis.title.y = element_text(angle = 90, size = 9, colour = "black")) + labs(x = "Longitude", y = "Latitude")+ scale_y_continuous(breaks = c(-20, 0, 20)) p[[i]]<-plotx } do.call(grid.arrange,c(p[1:3], ncol=1)) do.call(grid.arrange,c(p[4:6], ncol=1)) ################################################################################ ############################## BIVARIATE MAP CONSTRUCTION ###################### ################################################################################ clipExt <- extent(-130, 160, -23.5, 23.5) nquantilesforall<-5 colmat <- function(nquantiles = nquantilesforall, upperleft = "#0096EB", upperright = "#820050", bottomleft= "#BEBEBE", bottomright = "#FFE60F", xlab=expression(FD[Nut]), ylab=expression(FRed[Nut]), plotLeg = TRUE, saveLeg = TRUE) { my.data <- seq(0, 1, .01) my.class <- classInt::classIntervals(my.data, n = nquantiles, style = "quantile" ) my.pal.1 <- findColours(my.class, c(upperleft, bottomleft)) my.pal.2 <- findColours(my.class, c(upperright, bottomright)) col.matrix <- matrix(nrow = 101, ncol = 101, NA) for (i in 1:101) { my.col <- c(paste(my.pal.1[i]), paste(my.pal.2[i])) col.matrix[102 - i, ] <- findColours(my.class, my.col) } col.matrix.plot <- col.matrix %>% as.data.frame(.) %>% mutate("Y" = row_number()) %>% mutate_at(.tbl = ., .vars = vars(starts_with("V")), .funs = list(as.character)) %>% pivot_longer(data = ., cols = -Y, names_to = "X", values_to = "HEXCode") %>% mutate("X" = as.integer(sub("V", "", .$X))) %>% distinct(as.factor(HEXCode), .keep_all = TRUE) %>% mutate(Y = rev(.$Y)) %>% dplyr::select(-c(4)) %>% mutate("Y" = rep(seq(from = 1, to = nquantiles, by = 1), each = nquantiles), "X" = rep(seq(from = 1, to = nquantiles, by = 1), times = nquantiles)) %>% mutate("UID" = row_number()) # Use plotLeg if you want a preview of the legend if (plotLeg) { p <- ggplot(col.matrix.plot, aes(X, Y, fill = HEXCode)) + geom_raster() + scale_fill_identity() + coord_equal(expand = FALSE) + theme_void() + theme(aspect.ratio = 1, axis.title = element_text(size = 12, colour = "black",hjust = 0.5, vjust = 1), axis.title.y = element_text(angle = 90, hjust = 0.5)) + xlab(bquote(.(xlab) ~ symbol("\256"))) + ylab(bquote(.(ylab) ~ symbol("\256"))) print(p) assign( x = "BivLegend", value = p, pos = .GlobalEnv ) } # Use saveLeg if you want to save a copy of the legend if (saveLeg) { ggsave(filename = "bivLegend.pdf", plot = p, device = "pdf", path = "./", width = 4, height = 4, units = "in", dpi = 300) } seqs <- seq(0, 100, (100 / nquantiles)) seqs[1] <- 1 col.matrix <- col.matrix[c(seqs), c(seqs)] } bivariate.map <- function(rasterx, rastery, colormatrix = col.matrix, nquantiles = nquantilesforall, export.colour.matrix = TRUE, outname = paste0("colMatrix_rasValues", names(rasterx))) { quanmean <- getValues(rasterx) temp <- data.frame(quanmean, quantile = rep(NA, length(quanmean))) brks <- with(temp, quantile(temp, na.rm = TRUE, probs = c(seq(0, 1, 1 / nquantiles)) )) brks[-1] <- brks[-1] + seq_along(brks[-1]) * .Machine$double.eps r1 <- within(temp, quantile <- cut(quanmean, breaks = brks, labels = 2:length(brks), include.lowest = TRUE )) quantr <- data.frame(r1[, 2]) quanvar <- getValues(rastery) temp <- data.frame(quanvar, quantile = rep(NA, length(quanvar))) brks <- with(temp, quantile(temp, na.rm = TRUE, probs = c(seq(0, 1, 1 / nquantiles)) )) brks[-1] <- brks[-1] + seq_along(brks[-1]) * .Machine$double.eps r2 <- within(temp, quantile <- cut(quanvar, breaks = brks, labels = 2:length(brks), include.lowest = TRUE )) quantr2 <- data.frame(r2[, 2]) as.numeric.factor <- function(x) { as.numeric(levels(x))[x] } col.matrix2 <- colormatrix cn <- unique(colormatrix) for (i in 1:length(col.matrix2)) { ifelse(is.na(col.matrix2[i]), col.matrix2[i] <- 1, col.matrix2[i] <- which( col.matrix2[i] == cn )[1] ) } if (export.colour.matrix) { exportCols <- as.data.frame(cbind( as.vector(col.matrix2), as.vector(colormatrix), t(col2rgb(as.vector(colormatrix))) )) colnames(exportCols)[1:2] <- c("rasValue", "HEX") assign( x = outname, value = exportCols, pos = .GlobalEnv ) } cols <- numeric(length(quantr[, 1])) for (i in 1:length(quantr[, 1])) { a <- as.numeric.factor(quantr[i, 1]) b <- as.numeric.factor(quantr2[i, 1]) cols[i] <- as.numeric(col.matrix2[b, a]) } r <- rasterx r[1:length(r)] <- cols return(r) } # Define the number of breaks nBreaks <- nquantilesforall # Create the colour matrix col.matrix <- colmat(nquantiles = nBreaks, upperleft = "#0096EB", upperright = "#820050", bottomleft= "#BEBEBE", bottomright = "#FFE60F", xlab=expression(FD[Nut]), ylab=expression(FRed[Nut]), saveLeg = FALSE, plotLeg = TRUE) bivmap_Nut<-bivariate.map(FDFRED_mapspred_done$FD..Nutrients, FDFRED_mapspred_done$FRed..Nutrients, colormatrix=col.matrix,export.colour.matrix = TRUE, outname = "bivMapCols", nquantiles=nBreaks) bivmap_Mor<-bivariate.map(FDFRED_mapspred_done$FD..Morphology.Structure, FDFRED_mapspred_done$FRed..Morphology.Structure, colormatrix=col.matrix,export.colour.matrix = TRUE, outname = "bivMapCols", nquantiles=nBreaks) bivmap_pho<-bivariate.map(FDFRED_mapspred_done$FD..Photosynthesis, FDFRED_mapspred_done$FRed..Photosynthesis, colormatrix=col.matrix,export.colour.matrix = TRUE, outname = "bivMapCols", nquantiles=nBreaks) biv_stacked<-stack(bivmap_Mor,bivmap_Nut,bivmap_pho) p_biv1 <- list() for(i in 1:3){ figname=ifelse(names(biv_stacked)[[i]]=='FD..Photosynthesis','Photosynthesis', ifelse(names(biv_stacked)[[i]]=="FD..Nutrients",'Nutrients','Morphology/Structure')) rastmap<-biv_stacked[[i]] rastmap<-as.data.frame(biv_stacked[[i]], xy = TRUE)%>%tbl_df() %>%dplyr::rename("BivValue" = 3) %>%pivot_longer(., names_to = "Variable", values_to = "bivVal", cols = BivValue) plotx<-ggplot(rastmap, aes(x = x, y = y)) + ggtitle(figname) + geom_sf(data=world1,fill='gray100',col='gray80',inherit.aes = FALSE)+ geom_raster(aes(fill = bivVal)) + scale_fill_gradientn(colours = col.matrix, na.value = "transparent") + theme(text = element_text(size = 10, colour = "black"),panel.background = element_rect(fill = "white")) + borders(colour = "gray60", size = 0.05) + coord_sf(expand = FALSE, xlim = c(-120, 181), ylim = c(-23.5, 23.5))+ theme(legend.position = "none", legend.title = element_blank(), plot.background = element_blank(), strip.text = element_text(size = 12, colour = "black"), panel.border = element_blank(), axis.text.y = element_text(angle = 90, hjust = 0.5), axis.text = element_text(size = 9, colour = "black"), axis.title = element_text(size = 9, colour = "black"), axis.title.y = element_text(angle = 90, size = 9, colour = "black")) + labs(x = "Longitude", y = "Latitude")+ scale_y_continuous(breaks = c(-20, 0, 20)) p_biv1[[i]]<-plotx } #then add this one below to the list.... eraselater <- FDFRED_mapspred_done erase_fdmor<- rescale0to1(FDFRED_mapspred_done[[1]]) erase_fdnut<- rescale0to1(FDFRED_mapspred_done[[2]]) erase_fdpho<- rescale0to1(FDFRED_mapspred_done[[3]]) erase_fd_all<- sum(erase_fdmor,erase_fdnut,erase_fdpho) erase_fd_all_mean<- mean(erase_fdmor,erase_fdnut,erase_fdpho) erase_frmor<- rescale0to1(FDFRED_mapspred_done[[4]]) erase_frnut<- rescale0to1(FDFRED_mapspred_done[[5]]) erase_frpho<- rescale0to1(FDFRED_mapspred_done[[6]]) erase_fr_all<- sum(erase_frmor,erase_frnut,erase_frpho) erase_fr_all_mean<- mean(erase_frmor,erase_frnut,erase_frpho) bivmap_scaled_all_together<-bivariate.map(erase_fd_all, erase_fr_all, colormatrix=col.matrix, export.colour.matrix = TRUE, outname = "bivMapCols", nquantiles=nBreaks) bivmap_scaled_all_together_mean<-bivariate.map(erase_fd_all_mean, erase_fr_all_mean, colormatrix=col.matrix, export.colour.matrix = TRUE, outname = "bivMapCols", nquantiles=nBreaks) figname='Combined FD and FRed' rastmap_scaled<-as.data.frame(bivmap_scaled_all_together, xy = TRUE)%>%tbl_df() %>%dplyr::rename("BivValue" = 3) %>%pivot_longer(., names_to = "Variable", values_to = "bivVal", cols = BivValue) plotx_scaled<-ggplot(rastmap_scaled, aes(x = x, y = y)) + ggtitle(figname) + geom_sf(data=world1,fill='gray100',col='gray80',inherit.aes = FALSE)+ geom_raster(aes(fill = bivVal)) + scale_fill_gradientn(colours = col.matrix, na.value = "transparent") + theme(text = element_text(size = 10, colour = "black"),panel.background = element_rect(fill = "white")) + borders(colour = "gray60", size = 0.05) + coord_sf(expand = FALSE, xlim = c(-120, 181), ylim = c(-23.5, 23.5))+ theme(legend.position = "none", legend.title = element_blank(), plot.background = element_blank(), strip.text = element_text(size = 12, colour = "black"), panel.border = element_blank(), axis.text.y = element_text(angle = 90, hjust = 0.5), axis.text = element_text(size = 9, colour = "black"), axis.title = element_text(size = 9, colour = "black"), axis.title.y = element_text(angle = 90, size = 9, colour = "black")) + labs(x = "Longitude", y = "Latitude")+ scale_y_continuous(breaks = c(-20, 0, 20)) p_biv1[[4]]<-plotx_scaled do.call(grid.arrange,c(p_biv1, ncol=1)) ############################################################################################################################# ##### COMPARISONS BETWEEN BIVARIATE MAPS WHEN MODELLS ARE FITTED BY LEAVING ONE REGION OUT (Americas, Africa or Asia) ############################################################################################################################# library(spatialEco) conflict_prefer("theme_map", "ggthemes") #### CORRELATION BETWEEN FULL MODEL PREDICTIONS AND 1 continent out MODELS ### FD MORPHOLOGY #### And do it in the same way for Nutrients and Photosynthesis and for the FRed group dataFDMorph<-Fd%>%subset(Trait=='FDis_mor_BA'&Value!=0) set.seed(103) scaled_dataFDMorph <-scale(dataFDMorph[,c(6,8,9,10,11,15,16)]) scaleListFD<- data.frame(name='FDmorphology',scale = attr(scaled_dataFDMorph, "scaled:scale"),center = attr(scaled_dataFDMorph, "scaled:center"))%>%rownames_to_column()%>%rename(varia=rowname) dataFDMorph <- dataFDMorph%>%select(1:5,13,14)%>%cbind(scaled_dataFDMorph) fd_mor_prediction_maps_AMout<- stan_glm(log(Value) ~ FT_MCWD*FT_VPD+ Ch_MCWD*Ch_VPD + cec_0to30cm + clay_0to30cm , data=dataFDMorph%>%subset(Loc=="Australia"|Loc=="Ghana"|Loc=="Malaysia"|Loc=="Gabon"), algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99))#, fd_mor_prediction_maps_AFout<- stan_glm(log(Value) ~ FT_MCWD*FT_VPD+ Ch_MCWD*Ch_VPD + cec_0to30cm + clay_0to30cm + Size , data=dataFDMorph%>%subset(Loc=="Australia"|Loc=="Santarem"|Loc=="Malaysia"|Loc=="Peru"|Loc=="Colombia"|Loc=="NovaX"|Loc=="Mexico"|Loc=="Panama"), algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99))#, fd_mor_prediction_maps_ASout<- stan_glm(log(Value) ~ FT_MCWD*FT_VPD+ Ch_MCWD*Ch_VPD + cec_0to30cm + clay_0to30cm + Size , data=dataFDMorph%>%subset(Loc=="Peru"|Loc=="Ghana"|Loc=="Santarem"|Loc=="Gabon"|Loc=="Colombia"|Loc=="NovaX"|Loc=="Mexico"|Loc=="Panama"), algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99))#, pantro_all_FDmor<-pantro_all1%>%mutate(Ch_VPD= scale(Ch_VPD, center = attr(scaled_dataFDMorph, "scaled:center")[1],scale =attr(scaled_dataFDMorph, "scaled:scale")[1]), Ch_MCWD=scale(Ch_MCWD,center = attr(scaled_dataFDMorph, "scaled:center")[2],scale =attr(scaled_dataFDMorph, "scaled:scale")[2]), Size= scale(Size, center = attr(scaled_dataFDMorph, "scaled:center")[3],scale =attr(scaled_dataFDMorph, "scaled:scale")[3]), FT_MCWD=scale(FT_MCWD,center = attr(scaled_dataFDMorph, "scaled:center")[4],scale =attr(scaled_dataFDMorph, "scaled:scale")[4]), FT_VPD= scale(FT_VPD, center = attr(scaled_dataFDMorph, "scaled:center")[5],scale =attr(scaled_dataFDMorph, "scaled:scale")[5]), cec_0to30cm= scale(cec_0to30cm, center = attr(scaled_dataFDMorph, "scaled:center")[6],scale =attr(scaled_dataFDMorph, "scaled:scale")[6]), clay_0to30cm= scale(clay_0to30cm, center = attr(scaled_dataFDMorph, "scaled:center")[7],scale =attr(scaled_dataFDMorph, "scaled:scale")[7]))%>% na.omit() fd_mor_mapspred_AMout<-posterior_predict(fd_mor_prediction_maps_AMout, pantro_all_FDmor, draws = 700)%>%as.data.frame()%>%colMeans()%>%melt()%>%cbind(pantro_all_FDmor)%>%rename(AMout_prediction=value)%>%select(2,3,1)%>%rasterFromXYZ()%>%as.data.frame(xy=T) fd_mor_mapspred_AFout<-posterior_predict(fd_mor_prediction_maps_AFout, pantro_all_FDmor, draws = 700)%>%as.data.frame()%>%colMeans()%>%melt()%>%cbind(pantro_all_FDmor)%>%rename(AFout_prediction=value)%>%select(2,3,1)%>%rasterFromXYZ()%>%as.data.frame(xy=T) fd_mor_mapspred_ASout<-posterior_predict(fd_mor_prediction_maps_ASout, pantro_all_FDmor, draws = 700)%>%as.data.frame()%>%colMeans()%>%melt()%>%cbind(pantro_all_FDmor)%>%rename(ASout_prediction=value)%>%select(2,3,1)%>%rasterFromXYZ()%>%as.data.frame(xy=T) fd_mor_mapspred_FULL<-fd_mapspred_done[[3]]%>%as.data.frame(xy=T) mor_maps_comparison<-cbind(fd_mor_mapspred_FULL,fd_mor_mapspred_AMout%>%select(3),fd_mor_mapspred_AFout%>%select(3),fd_mor_mapspred_ASout%>%select(3))%>%drop_na()%>% rename(FDmor_FULL_prediction=FDis_mor_BA, FDmor_AMout_prediction=AMout_prediction,FDmor_AFout_prediction=AFout_prediction,FDmor_ASout_prediction=ASout_prediction) cor_fd_mor_maps<-mor_maps_comparison%>%select(3:6)%>%cor(method = "spearman") ggcorrplot(cor_fd_mor_maps, hc.order = TRUE, type = "lower",lab = TRUE) ### CHECK RELATIONSHIP BETWEEN NUMBER OF RECORDS AND CORRELATIONS ######## comparisons<-c('FULL','AMout','AFout','ASout') records<-c(nrow(dataFDMorph),nrow(dataFDMorph%>%subset(Loc=="Australia"|Loc=="Ghana"|Loc=="Malaysia"|Loc=="Gabon")),nrow(dataFDMorph%>%subset(Loc=="Australia"|Loc=="Santarem"|Loc=="Malaysia"|Loc=="Peru"|Loc=="Colombia"|Loc=="NovaX"|Loc=="Mexico"|Loc=="Panama")),nrow(dataFDMorph%>%subset(Loc=="Peru"|Loc=="Ghana"|Loc=="Santarem"|Loc=="Gabon"|Loc=="Colombia"|Loc=="NovaX"|Loc=="Mexico"|Loc=="Panama"))) correl <-c(cor_fd_mor_maps[1],cor_fd_mor_maps[2],cor_fd_mor_maps[3],cor_fd_mor_maps[4]) dfcmor<-data.frame(comparisons,records,correl) cor_fd_mormaps<-cor(dfcmor$records, dfcmor$correl) ############################################################################################################################# ############################################################################################################################# # FD, FRed and AGB analysis ############################################################################################################################# library(readxl) conflict_prefer("extract", "raster") agb_bennet<-read_excel('1_AGB_Plotslocations_ready.xlsx',1)%>%select(1,2,4) names(agb_bennet)<-c('Plot', 'AGB_pre', 'AGB_post') plots_bennet<-read_excel('1_AGB_Plotslocations_ready.xlsx',2)%>% select(1:5,7)%>%rename(Plot=`Plot Code`,Lat=Lat.,Long=Long., PlotArea_ha= `Area (ha)`) all_bennet <- merge(plots_bennet,agb_bennet,by='Plot')%>%mutate(AGB_change=AGB_post-AGB_pre) boxplot(all_bennet$AGB_change) hist(all_bennet$AGB_change) outliers <- boxplot(all_bennet$AGB_change, plot=FALSE)$out all_bennet <- all_bennet[-which(all_bennet$AGB_change %in% outliers),] boxplot(all_bennet$AGB_change) coordinates(all_bennet)<-~Long+Lat ###### GET THE FD, FRed, Bivariate maps names(FDFRED_mapspred_done) stack_for_bennet<-stack(bivmap_Mor, bivmap_Nut, bivmap_pho) names(stack_for_bennet)<-c("BivarMor", "BivarNut","BivarPho") all_stack_for_bennet<-stack(FDFRED_mapspred_done,stack_for_bennet) ##### Extract FD, FRed, Bivariate values to points of plots plots_stack_bennet<-raster::extract(all_stack_for_bennet, all_bennet)%>%as.data.frame()%>%cbind(all_bennet)%>% select(10:18,1:9)%>%drop_na()%>% melt(id.vars=c('Plot','Country','Cluster','Lat','Long',"PlotArea_ha",'AGB_pre','AGB_post','AGB_change'))%>% mutate(AGB_change_rel= ((AGB_post-AGB_pre)/AGB_pre)*100) ggplot(plots_stack_bennet,aes(value,AGB_change))+geom_point()+ geom_smooth(method = lm)+facet_wrap(~variable, scales = 'free') ggplot(plots_stack_bennet,aes(value,AGB_change_rel))+geom_point()+ geom_smooth(method = lm)+facet_wrap(~variable, scales = 'free') ggplot(plots_stack_bennet,aes(value,AGB_change))+geom_point()+ geom_smooth(method = loess)+facet_wrap(~variable, scales = 'free') ggplot(plots_stack_bennet,aes(value,AGB_change_rel))+geom_point()+ geom_smooth(method = loess)+facet_wrap(~variable, scales = 'free') trops<-pantro_all1[,c(1,2,6,4,5,3)]%>%rasterFromXYZ() tri2<-plots_stack_bennet%>%reshape2::dcast(Plot+Country+Cluster+Lat+Long+PlotArea_ha+ AGB_pre+AGB_post+AGB_change+AGB_change_rel ~ variable) tri3<-tri2[,c(1,4,5)] coordinates(tri3)<-~Long+Lat tri3<-extract(trops, tri3)%>%cbind(tri2) #correlations cor_agb<-tri3%>%select(1:4,15:23) cor_agb_m = cor(cor_agb, method = "spearman") ggcorrplot(cor_agb_m, hc.order = TRUE, type = "lower",lab = TRUE) #Because of the high correlations between FD, FRed and the Bivariate maps I excluded Bivariate info in the models #but carry our the models with interactions with Delta MCWD and Delta VPD #I did the models using the % change AGB as suggested by Reviewer 1 scaled_agb <-scale(tri3[,c(1:4,10,15:23)]) scaleList_agb<- data.frame(name='agb', scale = attr(scaled_agb, "scaled:scale"), center = attr(scaled_agb, "scaled:center"))%>%rownames_to_column()%>%rename(varia=rowname) tri3 <- tri3%>%select(c(5:9,11:14))%>%cbind(scaled_agb) scalebackdataframe <-tri3%>%#### THIS DATAFRAME WILL BE USED LATER FOR PLOTTING THE POINTS OF THE RAW DATA ONLY mutate(Ch_VPD= Ch_VPD*scaleList_agb[2,3] + scaleList_agb[2,4], Ch_MCWD= Ch_MCWD*scaleList_agb[1,3] + scaleList_agb[1,4], FT_VPD= FT_VPD*scaleList_agb[4,3] + scaleList_agb[4,4], FT_MCWD= FT_MCWD*scaleList_agb[3,3] + scaleList_agb[3,4], PlotArea_ha= PlotArea_ha* scaleList_agb[5,3] + scaleList_agb[5,4], FD..Morphology.Structure= FD..Morphology.Structure*scaleList_agb[6,3] + scaleList_agb[6,4], FRed..Morphology.Structure= FRed..Morphology.Structure*scaleList_agb[9,3] + scaleList_agb[9,4], FD..Nutrients= FD..Nutrients*scaleList_agb[7,3] + scaleList_agb[7,4], FRed..Nutrients= FRed..Nutrients*scaleList_agb[10,3] + scaleList_agb[10,4], FD..Photosynthesis= FD..Photosynthesis*scaleList_agb[8,3] + scaleList_agb[8,4], FRed..Photosynthesis= FRed..Photosynthesis*scaleList_agb[11,3] + scaleList_agb[11,4], BivarMor=BivarMor*scaleList_agb[12,3] + scaleList_agb[12,4], BivarNut=BivarNut*scaleList_agb[13,3] + scaleList_agb[13,4], BivarPho=BivarPho*scaleList_agb[14,3] + scaleList_agb[14,4]) ############################################ ###### AGB models for MORPHOLOGY/STRUCTURE # ############################################ agb_fd_mor1<-stan_glm(AGB_change ~ ((FD..Morphology.Structure) + (FRed..Morphology.Structure))*Ch_MCWD+ ((FD..Morphology.Structure) + (FRed..Morphology.Structure))*Ch_VPD+ (PlotArea_ha), data=tri3, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99)) agb_fd_mor1_post <-describe_posterior(agb_fd_mor1, rope_ci=.9, ci=.90)%>%as.data.frame()%>%mutate(Trait='Morpho', R2=bayes_R2(agb_fd_mor1) %>% median(), R2cond=as.numeric(r2_bayes(agb_fd_mor1,ci = 0.90)[1]), R2AdjLoo=as.numeric(r2_loo(agb_fd_mor1)))%>% mutate(NPlots=nrow(tri3)) ########################################### ###### AGB models for Nutrients ########################################### agb_fd_nut1<- stan_glm(AGB_change ~ ((FD..Nutrients) + (FRed..Nutrients))*Ch_MCWD + ((FD..Nutrients) + (FRed..Nutrients))*Ch_VPD + (PlotArea_ha), data=tri3, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99)) describe_posterior(agb_fd_nut1, ci=.90) agb_fd_nut1_post <-describe_posterior(agb_fd_nut1, rope_ci=.9, ci=.90)%>%as.data.frame()%>%mutate(Trait='Nutrients', R2=bayes_R2(agb_fd_nut1) %>% median(), R2cond=as.numeric(r2_bayes(agb_fd_nut1,ci = 0.90)[1]), R2AdjLoo=as.numeric(r2_loo(agb_fd_nut1)))%>% mutate(NPlots=nrow(tri3)) ########################################### ###### AGB models for Photosynthesis ########################################### agb_fd_pho1<- stan_glm(AGB_change ~((FD..Photosynthesis) + (FRed..Photosynthesis))*Ch_MCWD+ ((FD..Photosynthesis) + (FRed..Photosynthesis))*Ch_VPD+ (PlotArea_ha), data=tri3, algorithm = "sampling", cores = 13, seed = 17,chains=3, iter=5000, control = list(adapt_delta = 0.99)) describe_posterior(agb_fd_pho1, ci=.90) agb_fd_pho1_post <-describe_posterior(agb_fd_pho1, rope_ci=.9, ci=.90)%>%as.data.frame()%>%mutate(Trait='Photosynthesis', R2=bayes_R2(agb_fd_pho1) %>% median(), R2cond=as.numeric(r2_bayes(agb_fd_pho1,ci = 0.90)[1]), R2AdjLoo=as.numeric(r2_loo(agb_fd_pho1)))%>% mutate(NPlots=nrow(tri3)) ########################################################################## # FD and FRed maps with points used to fit the models # and that also show climatic space not available in model fitting ########################################################################## RecFDmorpho<-trial_fd%>%subset(Trait=="FDis_mor_BA") RecFDmorpho<-left_join(RecFDmorpho, cen[,c(1,3,4)]) coordinates(RecFDmorpho)<-~x+y RecFDnutrients<-trial_fd%>%subset(Trait=="FDis_nut_BA") RecFDnutrients<-left_join(RecFDnutrients, cen[,c(1,3,4)]) coordinates(RecFDnutrients)<-~x+y RecFDPhoto<-trial_fd%>%subset(Trait=="FDis_pho_BA") RecFDPhoto<-left_join(RecFDPhoto, cen[,c(1,3,4)]) coordinates(RecFDPhoto)<-~x+y r2_plots1<-levelplot(FDFRED_mapspred_done[[1]],maxpixels = 2e5, margin=FALSE, colorkey=list( axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = 'FD Morphology/Structure') + latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))+ latticeExtra::layer(sp.polygons(RecFDmorpho, col='blue')) r2_plots2<-levelplot(FDFRED_mapspred_done[[2]],maxpixels = 2e5, margin=FALSE, colorkey=list( axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = 'FD Nutrients') +#(69 records) latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))+ latticeExtra::layer(sp.polygons(RecFDnutrients, col='blue')) r2_plots3<-levelplot(FDFRED_mapspred_done[[3]],maxpixels = 2e5, margin=FALSE, colorkey=list( axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = 'FD Photosynthesis') + #(21 records) latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))+ latticeExtra::layer(sp.polygons(RecFDPhoto, col='blue',offset=5)) grid.arrange(r2_plots1,r2_plots2,r2_plots3, ncol=1) #FRed RecFRedmorpho<-trial_fr%>%subset(Trait=="RStar_mor_BA") RecFRedmorpho<-left_join(RecFRedmorpho, cen[,c(1,3,4)]) coordinates(RecFRedmorpho)<-~x+y RecFRednutrients<-trial_fr%>%subset(Trait=="RStar_nut_BA") RecFRednutrients<-left_join(RecFRednutrients, cen[,c(1,3,4)]) coordinates(RecFRednutrients)<-~x+y RecFRedPhoto<-trial_fr%>%subset(Trait=="RStar_pho_BA") RecFRedPhoto<-left_join(RecFRedPhoto, cen[,c(1,3,4)]) coordinates(RecFRedPhoto)<-~x+y r2_plots4<-levelplot(FDFRED_mapspred_done[[4]],maxpixels = 2e5, margin=FALSE, colorkey=list( axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = 'FRed Morphology/Structure') + latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))+ latticeExtra::layer(sp.polygons(RecFRedmorpho, col='blue')) r2_plots5<-levelplot(FDFRED_mapspred_done[[5]],maxpixels = 2e5, margin=FALSE, colorkey=list( axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = 'FRed Nutrients') + latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))+ latticeExtra::layer(sp.polygons(RecFRednutrients, col='blue')) r2_plots6<-levelplot(FDFRED_mapspred_done[[6]],maxpixels = 2e5, margin=FALSE, colorkey=list( axis.line=list(col='black')), par.settings=list( axis.line=list(col='transparent')), scales=list(draw=T), col.regions=colr2, main = 'FRed Photosynthesis') + latticeExtra::layer(sp.polygons(worldmap, lwd=.2,fill='gray50', col='gray60'), under = TRUE)+ latticeExtra::layer(sp.polygons(worldmap, lwd=.2, col='gray60'))+ latticeExtra::layer(sp.polygons(RecFRedPhoto, col='blue',offset=5)) grid.arrange(r2_plots4,r2_plots5,r2_plots6, ncol=1) #####FIND OUT REPRESENTATIVENES OF CLIMATE IN OBSERVATIONS AND PREDICTIONS- range(RecFDnutrients$FT_MCWD) range(RecFDnutrients$Ch_MCWD) range(RecFDnutrients$FT_VPD) range(RecFDnutrients$Ch_VPD) range(RecFDnutrients$cec_0to30cm) range(RecFDnutrients$clay_0to30cm) range(pantro_all1$FT_MCWD) range(pantro_all1$Ch_MCWD) range(pantro_all1$FT_VPD) range(pantro_all1$Ch_VPD) range(pantro_all1$cec_0to30cm, na.rm = T) range(pantro_all1$clay_0to30cm, na.rm = T) climateoverlap<-rasterFromXYZ(pantro_all1[,c(1:4,6)]) crs(climateoverlap)<-'+proj=longlat +datum=WGS84 +no_defs' real_ftmcwd_X<-stack('MCWDfullterm5817.tif') real_ftmcwd_X<-projectRaster(real_ftmcwd_X[[4]]*-1, climateoverlap) real_ftmcwd_X<-clamp(real_ftmcwd_X, upper=941.6207, useValues=TRUE)#1000 names(real_ftmcwd_X)<-'FT_MCWD' clay_pantropical1<-raster('Clay0_30cm_soilgrids_pantropical.tif') clay_pantropical1<-projectRaster(clay_pantropical1, treecover, method="bilinear") names(clay_pantropical1)<-'clay_0to30cm' cec_pantropical2<-raster('CEC0_30cm_soilgrids_pantropical.tif') cec_pantropical2<-projectRaster(cec_pantropical2, treecover, method="bilinear") names(cec_pantropical2)<-'cec_0to30cm' climateoverlap <-stack(climateoverlap,real_ftmcwd_X)%>%projectRaster(treecover)%>%mask(treecover)%>%stack(cec_pantropical2,clay_pantropical1)%>%mask(treecover) NotOverlap_FTMCWD<-raster::reclassify(climateoverlap$FT_MCWD,c(0,941.6207,1, 941.6208,Inf,0)) NotOverlap_FTVPD<-raster::reclassify(climateoverlap$FT_VPD,c(0.2459565,1.2237805,1 , 0,0.2459564,0, 1.2237806,Inf,0)) NotOverlap_ChMCWD<-raster::reclassify(climateoverlap$Ch_MCWD,c(-190.6885, 156.9886,1 , -800,-190.6885,0, 156.9886,Inf,0)) NotOverlap_ChVPD<-raster::reclassify(climateoverlap$Ch_VPD,c(-0.06395138, 0.08440180,1 , -0.16,-0.06395139,0, 0.08440181,.2,0)) NotOverlap_cec<-raster::reclassify(climateoverlap$cec_0to30cm,c(7.767297, 34.000000,1 , 0,7.767297,0, 34,Inf,0)) NotOverlap_clay<-raster::reclassify(climateoverlap$clay_0to30cm,c(20.57466, 50.71682,1 , 0,20.57466,0, 50.71682,Inf,0)) NoOverlap_all<-stack(NotOverlap_FTMCWD, NotOverlap_FTVPD,NotOverlap_ChMCWD,NotOverlap_ChVPD,NotOverlap_cec,NotOverlap_clay) plot(NotOverlap_clay,col=c('darkred','white')) library(scico) q <- list() for(i in 1:6){ figname=ifelse(names(NoOverlap_all)[[i]]=="FT_MCWD",'Full term MCWD (mm)', ifelse(names(NoOverlap_all)[[i]]=="FT_VPD",'Full termVPD (kPa)', ifelse(names(NoOverlap_all)[[i]]=="Ch_MCWD",'Delta MCWD (mm)', ifelse(names(NoOverlap_all)[[i]]=="Ch_VPD",'Delta VPD (kPa)', ifelse(names(NoOverlap_all)[[i]]=="cec_0to30cm",'CEC (mmol (c)/Kg','Clay (%)'))))) rastmap<-NoOverlap_all[[i]] rastmap<-rastmap%>%as.data.frame(xy=T)%>%na.omit() names(rastmap)<-c("x","y","Val") rastmap<-rastmap%>%subset(Val==0) plotx<-ggplot(rastmap, aes(x = x, y = y)) + ggtitle(figname)+ geom_sf(data=world1,fill='gray80',col='gray90',inherit.aes = FALSE)+ geom_raster(aes(fill = Val)) + scale_colour_manual(values = 'darkblue')+ theme(text = element_text(size = 10, colour = "black"),panel.background = element_rect(fill = "white")) + borders(colour = "gray60", size = 0.05) + coord_sf(expand = FALSE, xlim = c(-120, 181), ylim = c(-23.5, 23.5))+ theme(legend.position = "right", legend.title = element_blank(), plot.background = element_blank(), strip.text = element_text(size = 12, colour = "black"), panel.border = element_blank(), axis.text.y = element_text(angle = 90, hjust = 0.5), axis.text = element_text(size = 9, colour = "black"), axis.title = element_text(size = 9, colour = "black"), axis.title.y = element_text(angle = 90, size = 9, colour = "black")) + labs(x = "Longitude", y = "Latitude")+ scale_y_continuous(breaks = c(-20, 0, 20)) q[[i]]<-plotx } do.call(grid.arrange,c(q, ncol=2)) ggarrange(q[[1]] + theme(axis.title.x = element_blank(),legend.position = "none") + ggtitle('MCWD (mm)'), q[[2]] + theme(axis.title.x = element_blank(),legend.position = "none") + ggtitle('VPD (kPa)'), q[[3]] + theme(axis.title.x = element_blank(),legend.position = "none") + ggtitle('ΔMCWD (mm)'), q[[4]] + theme(axis.title.x = element_blank(),legend.position = "none") + ggtitle('ΔVPD (kPa)'), q[[5]] + theme(legend.position = "none") + ggtitle('CEC (mmol (c)/Kg)'), q[[6]] + theme(legend.position = "none") + ggtitle('Clay (%)'), ncol=2,nrow = 3, align = c( "hv")) ############################################## # Save the data ############################ save.image("FD_FRed_analysis.RData")