#--------------------------------------------------------------------------------------------------# #------CODE developed for analyses presented in Schlaepfer et al.: Future regeneration of big sagebrush support predicted changes in habitat suitability at trailing and leading edges #------CODE developed and written by # Daniel R Schlaepfer (dschlaep@uwyo.edu): 2014 # for contact and further information see also: sites.google.com/site/drschlaepfer #------DISCLAIMER: This program is distributed in the hope that it will be useful, #but WITHOUT ANY WARRANTY; without even the implied warranty of #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. #--------------------------------------------------------------------------------------------------# #------------------------------- #---R packages libraries <- c("RSQLite", "MASS", "seqinr", "energy") l <- lapply(libraries, FUN=function(lib) stopifnot(require(lib, character.only=TRUE, quietly=TRUE))) #---Directories dir.prj <- "/Users/drschlaep/Documents/drschlaepfer/2_Research/200907_UofWyoming_PostDoc/Projects_My/Product_NCCSC/BigSagebrushRegeneration_LeadingTrailingEdge/5_Manuscript/Code_Data" dir.dat <- file.path(dir.prj, "Data") dir.create(dir.out <- file.path(dir.prj, "Figures&Tables")) #------------------------------- #---Helper functions #Convert DOY of gregorian calendar to water year GregToWaterYear <- function(dat, vars){ dat[, vars, ] <- ifelse(dat[, vars, ] > 273, dat[, vars, ] - 273, dat[, vars, ] + 91) return(dat) } #Draws contours at alpha*100% of points (x, y) based on code by Pascal Haenggi 2008 (https://stat.ethz.ch/pipermail/r-help/2008-June/166079.html) draw.contourCIpoints <- function(x, y, alpha=0.95, col="black", lwd=1, fill=FALSE){ require(MASS) if(sd(x) > 0 & sd(y) > 0){ xlimP <- c(min(x), max(x)) ylimP <- c(min(y), max(y)) xlimKd <- c(xlimP[1]*(1-sign(xlimP[1])*1), xlimP[2]*(1+sign(xlimP[2])*1)) ylimKd <- c(ylimP[1]*(1-sign(ylimP[1])*1), ylimP[2]*(1+sign(ylimP[2])*1)) h <- c(bandwidth.nrd(x), bandwidth.nrd(y)) if(any(temp <- h == 0)) h[temp] <- sqrt(.Machine$double.eps) kerneld <- kde2d(x, y, h=h, n=100, lims = c(xlimKd, ylimKd)) pp <- array() for (i in 1:length(x)){ z.x <- max(which(kerneld$x < x[i])) z.y <- max(which(kerneld$y < y[i])) pp[i] <- kerneld$z[z.x, z.y] } confidencebound <- quantile(pp, 1-alpha, na.rm = TRUE) if(fill) { cl <- contourLines(kerneld, levels = confidencebound) require(seqinr) for(cli in seq(along=cl)){ polygon(x=cl[[cli]]$x, y=cl[[cli]]$y, border=NA, col=col2alpha(col, 0.3) ) } } else { contour(kerneld, levels=confidencebound, col=col, lwd=lwd, drawlabels=FALSE, xlim=xlimP, ylim=ylimP, add = TRUE) } } } #------------------------------- #---Database functionality name.dbScen <- "dbTables_Schlaepfer_BigSagebrushRegeneration_LeadingTrailingEdges.sqlite3" drv <- dbDriver("SQLite") addHeaderToWhereClause <- function(whereClause, headers=headerFields){ temp1 <- res <- strsplit(whereClause, split=" ", fixed=TRUE)[[1]] #Locate all "Label='x'" temp1F <- strsplit(temp1, split="=", fixed=TRUE) temp2 <- temp1[ielem <- grepl("=", temp1) & !grepl("header.", temp1, fixed=TRUE) & sapply(temp1F, FUN=function(ch) ch[1] %in% headers)] if(length(temp2) > 0) res[ielem] <- paste0("header.", temp2) #add 'header.' return(paste(res, collapse=" ")) } #Get data of variables in the overall aggregation table for one of the scenarios get.SeveralOverallVariables_Scenario <- function(responseName, MeanOrSD="Mean", scenario="Current", whereClause=NULL){ dat <- outOrder <- iColumns.iTable <- iColumns.header <- NULL if(length(responseName) > 0){ con <- dbConnect(drv, file.path(dir.dat, name.dbScen)) iTable <- (temp <- dbListTables(con))[grepl(pattern=paste0("Overall_", MeanOrSD), x=temp, ignore.case=T, fixed=FALSE)] if(length(iTable) == 1){ fields.header <- dbListFields(con,"header") fields.iTable <- dbListFields(con, iTable) responseNameC <- responseName if("P_id" %in% responseName) { addPid <- TRUE responseName<-responseName[!(responseName %in% "P_id")] } else { addPid <- FALSE } if(length(responseName) > 0) { responseName <- gsub(".", "_", responseName, fixed=TRUE) for(i in seq_along(responseName)){ iColumns.iTable <- c(iColumns.iTable, fields.iTable[grepl(pattern=responseName[i], x=gsub(".", "_", fields.iTable, fixed=TRUE), fixed=FALSE)]) iColumns.header <- c(iColumns.header, fields.header[grepl(pattern=responseName[i], x=gsub(".", "_", fields.header, fixed=TRUE), fixed=FALSE)]) } fields.iTable <- fields.iTable[!(fields.iTable %in% "P_id")] for(i in seq_along(responseNameC)) { outOrder <- c(outOrder, fields.iTable[grepl(pattern=responseNameC[i], x=gsub(".", "_", fields.iTable, fixed=TRUE), fixed=FALSE)], fields.header[grepl(pattern=responseNameC[i], x=gsub(".", "_", fields.header, fixed=TRUE), fixed=FALSE)]) } } if(length(iColumns.header) > 0 | length(iColumns.iTable) > 0 | addPid){ if(length(whereClause) > 0){ sql <- paste0("SELECT ", if(addPid) paste("header.P_id AS P_id", if(length(iColumns.header)>0 | length(iColumns.iTable)>0) ", ",sep=""), if(length(iColumns.header)>0) paste(paste0(paste0("\"", iColumns.header, "\"",sep=""), collapse=", "), if(length(iColumns.iTable)>0) ", ",sep=""), if(length(iColumns.iTable)>0) paste0(paste0("\"", iColumns.iTable, "\"",sep=""), collapse=", "), " FROM ", iTable, " INNER JOIN header ON ",iTable,".P_id=header.P_id WHERE header.Scenario=", shQuote(scenario), " AND ", addHeaderToWhereClause(whereClause), " ORDER BY header.P_id;") } else { sql <- paste0("SELECT ", if(addPid) paste("header.P_id AS P_id", if(length(iColumns.header)>0 | length(iColumns.iTable)>0) ", ",sep=""), if(length(iColumns.header)>0) paste(paste0(paste0("\"", iColumns.header, "\"",sep=""), collapse=", "), if(length(iColumns.iTable)>0) ", ",sep=""), if(length(iColumns.iTable)>0) paste0(paste0("\"", iColumns.iTable, "\"",sep=""), collapse=", "), " FROM ", iTable, " INNER JOIN header ON ",iTable,".P_id=header.P_id WHERE header.Scenario=", shQuote(scenario), " ORDER BY header.P_id;") } dat <- dbGetQuery(con, sql) } } dbDisconnect(con) } dat <- dat[,outOrder] return(dat) } #Get header and data for an entire table for one of the scenarios get.Table_Scenario <- function(responseName, MeanOrSD="Mean", scenario="Current", whereClause=NULL, header=FALSE){ dat <- NULL if(length(responseName) > 0){ con <- dbConnect(drv, file.path(dir.dat, name.dbScen)) iTable <- (temp <- dbListTables(con))[grepl(pattern=paste0(responseName, "_", MeanOrSD), x=temp, ignore.case=T, fixed=FALSE)] if(length(iTable) == 1){ fields <- dbListFields(con, iTable)[-1] if(length(whereClause) > 0){ sql <- paste0("SELECT ", if(header) "header.*, ", paste0(paste0("\"", fields, "\"",sep=""), collapse=", ") ," FROM ", iTable, " INNER JOIN header ON ",iTable,".P_id=header.P_id WHERE header.Scenario=", shQuote(scenario), " AND ", whereClause, " ORDER BY header.P_id;") } else { sql <- paste0("SELECT ", if(header) "header.*, ", paste0(paste0("\"", fields, "\"",sep=""), collapse=", ") ," FROM ", iTable, " INNER JOIN header ON ",iTable,".P_id=header.P_id WHERE header.Scenario=", shQuote(scenario), " ORDER BY header.P_id;") } dat <- dbGetQuery(con, sql) } dbDisconnect(con) } return(dat) } #Extract data for climate conditions and variables extractData <- function(var.Names, length=NULL, dimnames2=NULL){ if(is.null(length)) length <- length(var.Names) if(is.null(dimnames2)) dimnames2 <- var.Names array(unlist(lapply(climate.conditions, FUN=function(sc) get.SeveralOverallVariables_Scenario(responseName=var.Names, MeanOrSD="Mean", scenario=sc))), dim=c(nrow(locations), length, length(climate.conditions)), dimnames=list(NULL, dimnames2, climate.conditions)) } #------------------------------- #---Settings mo <- 1:12 #Months of year doy <- 1:365 # Days of year doysAtFirstOfMonths <- as.POSIXlt(seq(from=as.Date("2014-01-01"), to=as.Date("2014-12-31"), by="1 month"))$yday + 1 layerDepths_cm <- c(0, 10, 20, 60, 150) climate.conditions <- c("Current", "hybrid-delta.RCP45.CanESM2", "hybrid-delta.RCP45.CESM1-CAM5", "hybrid-delta.RCP45.CSIRO-Mk3-6-0", "hybrid-delta.RCP45.EC-EARTH", "hybrid-delta.RCP45.FGOALS-g2", "hybrid-delta.RCP45.FGOALS-s2", "hybrid-delta.RCP45.GFDL-CM3", "hybrid-delta.RCP45.GISS-E2-R", "hybrid-delta.RCP45.HadGEM2-CC", "hybrid-delta.RCP45.HadGEM2-ES", "hybrid-delta.RCP45.inmcm4", "hybrid-delta.RCP45.IPSL-CM5A-MR", "hybrid-delta.RCP45.MIROC-ESM", "hybrid-delta.RCP45.MIROC5", "hybrid-delta.RCP45.MPI-ESM-MR", "hybrid-delta.RCP45.MRI-CGCM3", "hybrid-delta.RCP85.CanESM2", "hybrid-delta.RCP85.CESM1-CAM5", "hybrid-delta.RCP85.CSIRO-Mk3-6-0", "hybrid-delta.RCP85.EC-EARTH", "hybrid-delta.RCP85.FGOALS-g2", "hybrid-delta.RCP85.FGOALS-s2", "hybrid-delta.RCP85.GFDL-CM3", "hybrid-delta.RCP85.GISS-E2-R", "hybrid-delta.RCP85.HadGEM2-CC", "hybrid-delta.RCP85.HadGEM2-ES", "hybrid-delta.RCP85.inmcm4", "hybrid-delta.RCP85.IPSL-CM5A-MR", "hybrid-delta.RCP85.MIROC-ESM", "hybrid-delta.RCP85.MIROC5", "hybrid-delta.RCP85.MPI-ESM-MR", "hybrid-delta.RCP85.MRI-CGCM3") rcps <- c("RCP45", "RCP85") #Color schemes colorLines <- c("red", rep("black", times=length(climate.conditions)-1)) lwdLines <- c(2, rep(1, times=length(climate.conditions)-1)) pchPoints <- c(16, rep(46, times=length(climate.conditions)-1)) pchPoints2 <- c(16, rep(4, times=length(climate.conditions)-1)) colorLinesByEdge <- matrix(c(c("dodgerblue3", rep("blue", times=length(climate.conditions)-1)), c("darkred", rep("red", times=length(climate.conditions)-1))), nrow=2, byrow=TRUE) colorLinesByEdge2 <- matrix(c(c("darkblue", rep("deepskyblue", times=length(climate.conditions)-1)), c("darkred", rep("orange", times=length(climate.conditions)-1))), nrow=2, byrow=TRUE) colorPointsByEdge <- matrix(c(c("darkblue", rep("deepskyblue", times=length(climate.conditions)-1)), c("darkred", rep("orange", times=length(climate.conditions)-1))), nrow=2, byrow=TRUE) colorPointsByEdge2 <- matrix(c(c("dodgerblue3", rep("blue", times=length(climate.conditions)-1)), c("darkred", rep("red", times=length(climate.conditions)-1))), nrow=2, byrow=TRUE) #---Variables and labels to analyse var.Edges <- c("Include_LeadingEdge", "Include_TrailingEdge") label.Edges <- c("Leading edge", "Trailing edge") var.ClimChange <- c("SWinput_ClimatePerturbations_PrcpMultiplier", "SWinput_ClimatePerturbations_TmaxAddand", "SWinput_ClimatePerturbations_TminAddand") label.ClimChange <- c("Scenario PPT/Current PPT (-)", "Scenario Tmax - Current Tmax (C)", "Scenario Tmin - Current Tmin (C)") var.Regeneration <- c( "ArtemisiaTridentataNotSpecifiedSubspecies.Germination.SuitableYears_fraction_mean", "ArtemisiaTridentataNotSpecifiedSubspecies.Seedlings1stSeason.SuitableYears_fraction_mean") label.Regeneration <- c("Germination success\n(Fraction of years)", "Seedling success\n(Fraction of years)") var.Explanatory <- c( "SWinput_Soil_topLayers_Sand_fraction", "SWinput_Soil_bottomLayers_Sand_fraction", "SWinput_Soil_topLayers_Clay_fraction", "SWinput_Soil_bottomLayers_Clay_fraction", "MAT_C_mean", "MAP_mm_mean", "PET_mm_mean", "Seasonality_monthlyTandPPT_PearsonCor_mean", "UNAridityIndex_Normals_none_mean", "SnowOfPPT_fraction_mean", "TempAir_m1_C_mean", "TminBelowNeg9degCwithoutSnowpack_days_mean", "TmaxAbovePos34degC_days_mean", "RelRecharge_topLayers_DailyMax_Fraction_mean", "RelRecharge_bottomLayers_DailyMax_Fraction_mean", "RelRecharge_topLayers_DailyMax_doy_mean", "RelRecharge_bottomLayers_DailyMax_doy_mean", "ThermalSnowfreeDryPeriods_SWPcrit2300kPa_topLayers_DrySpellsAllLayers_maxDuration_days_mean", "ThermalSnowfreeDryPeriods_SWPcrit2300kPa_bottomLayers_DrySpellsAllLayers_maxDuration_days_mean", "ThermalSnowfreeDryPeriods_SWPcrit2300kPa_topLayers_DrySpellsAtLeast10DaysAllLayers_Start_doy_mean", "ThermalSnowfreeDryPeriods_SWPcrit2300kPa_bottomLayers_DrySpellsAtLeast10DaysAllLayers_Start_doy_mean", "ThermalSnowfreeWetPeriods_SWPcrit2300kPa_topLayers_AvailableWater_mm_mean", "ThermalSnowfreeWetPeriods_SWPcrit2300kPa_bottomLayers_AvailableWater_mm_mean", "ThermalSnowfreeWetPeriods_SWPcrit2300kPa_topLayers_Duration_days_quantile0_5", "ThermalSnowfreeWetPeriods_SWPcrit2300kPa_bottomLayers_Duration_days_quantile0_5", "WetSoilPeriods_SWPcrit100kPa_NSadj_topLayers_AllLayersWet_Duration_LongestContinuous_days_mean", "WetSoilPeriods_SWPcrit100kPa_NSadj_bottomLayers_AllLayersWet_Duration_LongestContinuous_days_mean" ) label.Explanatory <- c( "Sand 0-20 cm (-)", "Sand 20-150 cm (-)", "Clay 0-20 cm (-)", "Clay 20-150 cm (-)", "MAT (C)", "MAP (mm)", "PET (mm)", "Seasonality =\ncor(monthlyT, monthly PPT)", "Aridity index = MAP / PET", "Snowfall/PPT (-)", "Mean January Temperature (C)", "Snowfree days with Tmin < -9 C", "Days with Tmax > 34 C", "Maximum recharge at 0-20 cm (-)", "Maximum recharge at 20-150 cm (-)", "Day of water year of\nmaximum recharge in 0-20 cm (-)", "Day of water year of\nmaximum recharge in 20-150 cm (-)", "Longest snowfree frost-free drought with\nSWP < -2.3 MPa at 0-20 cm (days)", "Longest snowfree frost-free drought with\nSWP < -2.3 MPa at 20-150 cm (days)", "ThermalSnowfreeDryPeriods_SWPcrit2300kPa_topLayers_DrySpellsAtLeast10DaysAllLayers_Start_doy_mean", "ThermalSnowfreeDryPeriods_SWPcrit2300kPa_bottomLayers_DrySpellsAtLeast10DaysAllLayers_Start_doy_mean", "Soil water with SWP > 2.3 Mpa during snowfree frost-free periods at 0-20 cm (mm)", "Soil water with SWP > 2.3 Mpa during snowfree frost-free periods at 20-150 cm (mm)", "Median suitable period with SWP > 2.3 MPa during snowfree frost-free periods at 0-20 cm (days)", "Median suitable period with SWP > 2.3 MPa during snowfree frost-free periods at 20-150 cm (days)", "Longest saturated period with SWP > 0.1 MPa at 0-20 cm (days)", "Longest saturated period with SWP > 0.1 MPa at 20-150 cm (days)" ) iExp2a <- c(1, 3, 5, 6, 9, 10) var.Explanatory2a <- c("ELEV_m", var.Explanatory[iExp2a]) label.Explanatory2a <- c("Elevation (m)", label.Explanatory[iExp2a]) iExp2b <- c(12, 13, 14, 16, 18, 19) var.Explanatory2b <- var.Explanatory[iExp2b] label.Explanatory2b <- label.Explanatory[iExp2b] var.SoilSpace <- c( "SWinput_Soil_topLayers_Sand_fraction", "SWinput_Soil_bottomLayers_Sand_fraction", "SWinput_Soil_topLayers_Clay_fraction", "SWinput_Soil_bottomLayers_Clay_fraction") label.SoilSpace <- c( "Sand 0-20 cm (-)", "Sand 20-150 cm (-)", "Clay 0-20 cm (-)", "Clay 20-150 cm (-)") var.ClimateSpace <- c("MAP_mm_mean", "MAT_C_mean") label.ClimateSpace <- c("MAP (mm)", "MAT (C)") #------------------------------- #---Extract data from database locations <- get.SeveralOverallVariables_Scenario(responseName=c("X_WGS84", "Y_WGS84"), scenario="Current") edge <- get.SeveralOverallVariables_Scenario(responseName=var.Edges, scenario="Current") climateDelta <- extractData(var.ClimChange, length=length(var.ClimChange)*length(mo), dimnames2=paste0(rep(var.ClimChange, each=length(mo)), "_m", mo)) climateSpace <- extractData(var.ClimateSpace) soilSpace <- extractData(var.SoilSpace) explanatories <- extractData(var.Explanatory) explanatories <- GregToWaterYear(dat=explanatories, vars=c("RelRecharge_topLayers_DailyMax_doy_mean", "RelRecharge_bottomLayers_DailyMax_doy_mean")) explanatories2a <- extractData(var.Explanatory2a) explanatories2b <- extractData(var.Explanatory2b) explanatories2b <- GregToWaterYear(dat=explanatories2b, vars=c("RelRecharge_topLayers_DailyMax_doy_mean")) regenerationResponse <- extractData(var.Regeneration) #------------------------ #---Analysis and figures #Climate space fig.2dVarSpaceByRCP <- function(datExp, labelsExp, xlim=NULL, ylim=NULL, points=TRUE, median=FALSE, jitter=FALSE, rcp, xlegend="bottomright", dir.fig, figname){ have.labelsCR <- grepl("\n", labelsExp) h.panel <- 3; w.panel <- 3; w.edge <- 0.5 + if(have.labelsCR[1]) 0.2 else 0; h.edge <- 0.45 + if(have.labelsCR[2]) 0.2 else 0 pdf(height=h.edge+h.panel*length(rcp), width=w.edge+w.panel*1, file=file.path(dir.fig, figname)) layout(matrix(1:((1+length(rcp)) * (1+1)), ncol=1+1, byrow=TRUE), heights=c(rep(h.panel, times=length(rcp)), h.edge), widths=c(w.edge, rep(w.panel, times=1))) cex <- 1 op <- par(mar=c(0.5, 0, 1, 0), mgp=c(1.6, 0.6, 0), las=0, xpd=NA, cex=cex) if(is.null(xlim)) xlim <- c(min(datExp[, 1, ]), max(datExp[, 1, ])) if(is.null(ylim)) ylim <- c(0, max(datExp[, 2, ])) for(ip in seq_along(rcp)){ iclims <- c(1, grep(rcp[ip], climate.conditions)) for(ix in 1:2){ for(isc in rev(iclims)){ for(ie in 1:ncol(edge)){ iedge <- which(edge[, ie] > 0) x <- datExp[iedge, 1, isc] y <- datExp[iedge, 2, isc] if(jitter){ x <- jitter(x, amount=0.001) y <- jitter(y, amount=0.001) } if(ix == 1 && ie == 1 && isc == tail(iclims, n=1)){ plot.new() plot(x, y, type="n", xlim=xlim, ylim=ylim, bty="l", xlab="", ylab="", axes=FALSE) if(ip == length(rcp)) mtext(side=1, line=1, text=labelsExp[1], padj=1, cex=cex) mtext(side=2, line=1.5, text=labelsExp[2], cex=cex) axis(side=1, labels=(ip == length(rcp))) axis(side=2, labels=TRUE) if(length(rcp) > 1) mtext(side=3, line=-1, adj=1, text=paste("Current vs", rcp[ip], " "), cex=cex) } if(median && ix == 1 || !median && ix == 2){ draw.contourCIpoints(x, y, alpha=0.9, col=if(median) colorLinesByEdge2[ie, isc] else colorLinesByEdge2[ie, isc], lwd=lwdLines[isc], fill=FALSE) } else { if(points) points(x, y, pch=pchPoints[isc], col=colorPointsByEdge2[ie, isc], cex=if(isc == 1) 0.5 else 1) if(median) points(median(x), median(y), pch=pchPoints2[isc], col=colorPointsByEdge2[ie, isc], cex=1) } } } } if(ip == 1){ lpos <- legend(x=xlegend, bty="n", cex=0.7*cex, legend=c("Leading edge: current", "Leading edge: future", "Trailing edge: current", "Trailing edge: future"), lwd=c(2, 1, 2, 1), col=if(median) c(t(colorLinesByEdge2[1:2, 1:2])) else c(t(colorLinesByEdge2[1:2, 1:2]))) pchs <- ifelse((temp <- if(median) pchPoints2 else pchPoints) == 16 | temp == 46, 21, temp)[rep(1:2, 2)] pcexs <- if(median) rep(0.7, 4) else rep(c(0.7, 0.3), 2) pcols <- if(median) c(t(colorPointsByEdge2[1:2, 1:2])) else c(t(colorPointsByEdge2[1:2, 1:2])) pbgs <- pcols; pbgs[c(1, 3)] <- "white" for(il in 1:4) points(x=lpos$text$x[il] - 2*0.7*par("cex")*xinch(par("cin")[1L], warn.log = FALSE), y=lpos$text$y[il], pch=pchs[il], col=pbgs[il], bg=pcols[il], cex=pcexs[il]*cex) } } par(op) dev.off() } #Figure 2 fig.2dVarSpaceByRCP(datExp=climateSpace, labelsExp=label.ClimateSpace, rcp="RCP85", dir.fig=dir.out, figname="Fig_ClimateSpaceRCP85.pdf") #Figure S3 fig.2dVarSpaceByRCP(datExp=climateSpace, labelsExp=label.ClimateSpace, rcp="RCP45", dir.fig=dir.out, figname="Fig_ClimateSpaceRCP45.pdf") tab.ClimateSpaceByRCP <- function(rcp, dir.tab, tabname){ nclims <- max(sapply(seq_along(rcp), FUN=function(ip) length(c(1, grep(rcp[ip], climate.conditions))))) resp <- c("Mean", "Quantile_5prc", "Quantile_95prc", "Spread_ScenToCur", "Delta_ScenToCur", "Delta_ScenMinusCur") vars <- dimnames(climateSpace)[2][[1]] exps <- c("Edge", "RCP", "Scenario") res <- data.frame(matrix(NA, ncol=length(exps) + length(vars) * length(resp), nrow=nclims, dimnames=list(NULL, c(exps, paste(rep(vars, each=length(resp)), rep(resp, times=length(vars)), sep="_"))))) for(ip in seq_along(rcp)){ iclims <- c(1, grep(rcp[ip], climate.conditions)) res[, 2] <- rcp[ip] res[, 3] <- climate.conditions[iclims] for(ie in 1:ncol(edge)){ iedge <- which(edge[, ie] > 0) res[, 1] <- colnames(edge)[ie] for(isc in seq_along(iclims)){ for(iv in seq_along(vars)){ x <- climateSpace[iedge, iv, iclims[isc]] res[isc, length(exps) + 1 + (iv - 1) * length(resp)] <- mean(x) res[isc, length(exps) + 2:3 + (iv - 1) * length(resp)] <- quantile(x, probs=c(0.05, 0.95)) res[isc, length(exps) + 4 + (iv - 1) * length(resp)] <- (res[isc, length(exps) + 3 + (iv - 1) * length(resp)] - res[isc, length(exps) + 2 + (iv - 1) * length(resp)]) / (res[1, length(exps) + 3 + (iv - 1) * length(resp)] - res[1, length(exps) + 2 + (iv - 1) * length(resp)]) res[isc, length(exps) + 5 + (iv - 1) * length(resp)] <- res[isc, length(exps) + 1 + (iv - 1) * length(resp)] / res[1, length(exps) + 1 + (iv - 1) * length(resp)] res[isc, length(exps) + 6 + (iv - 1) * length(resp)] <- res[isc, length(exps) + 1 + (iv - 1) * length(resp)] - res[1, length(exps) + 1 + (iv - 1) * length(resp)] } } res[length(iclims) + 1, -(1:3)] <- apply(res[-1, -(1:3)], 2, mean) if(ip == 1 && ie == 1){ write.csv(res, row.names=FALSE, file=file.path(dir.tab, tabname)) } else { write.table(res, append=TRUE, sep=",", dec=".", qmethod="double", col.names=FALSE, row.names=FALSE, file=file.path(dir.tab, tabname)) } res <- res[1:length(iclims), ] } } return(0) } #Table 1 tab.ClimateSpaceByRCP(rcp=rcps, dir.out, tabname="Tab_ClimateSpace.csv") #------------------------ #Climate vs soil space fig.ClimateSoilSpace <- function(dir.fig, figname){ h.panel <- 2; w.panel <- 2; w.edge <- 0.4; h.edge <- 0.4 cex <- 0.7 pdf(height=h.edge+h.panel*length(var.SoilSpace), width=w.edge+w.panel*length(var.ClimateSpace), file=file.path(dir.fig, figname)) layout(matrix(1:((1+length(var.SoilSpace)) * (1+length(var.ClimateSpace))), ncol=1+length(var.ClimateSpace), byrow=TRUE), heights=c(rep(h.panel, times=length(var.SoilSpace)), h.edge), widths=c(w.edge, rep(w.panel, times=length(var.ClimateSpace)))) op <- par(mar=c(0.5, 0, 0.5, 0), mgp=c(1.7, 0.75, 0), las=0, xpd=NA, cex=cex) ylim <- c(0, 1) for(iv2 in seq_along(var.SoilSpace)){ for(iv1 in seq_along(var.ClimateSpace)){ xlim <- c(min(climateSpace[, iv1, ]), max(climateSpace[, iv1, ])) for(ix in 1:2){ for(ie in 1:ncol(edge)){ iedge <- which(edge[, ie] > 0) x <- climateSpace[iedge, iv1, 1] y <- soilSpace[iedge, iv2, 1] if(ix == 1){ if(ie == 1){ if(iv1 == 1 && ie == 1) plot.new() plot(x, y, type="p", pch=16, cex=cex, col=colorPointsByEdge[ie, 1], xlim=xlim, ylim=ylim, bty="l", xlab=ifelse(iv2 == length(label.SoilSpace), label.ClimateSpace[iv1], ""), ylab=ifelse(iv1 == 1 && ie == 1, label.SoilSpace[iv2], ""), axes=FALSE) axis(side=1, labels=(iv2 == length(label.SoilSpace))) axis(side=2, labels=(iv1 == 1 && ie == 1)) } else { points(x, y, pch=16, cex=cex, col=colorPointsByEdge[ie, 1]) } } else { if(any(ok <- is.finite(x) & is.finite(y))) draw.contourCIpoints(x[ok], y[ok], alpha=0.9, col=colorLinesByEdge[ie, 1], lwd=lwdLines[1], fill=FALSE) } } } if(iv1 == 1 && iv2 == 1){ temp <- legend(x="bottomright", bty="n", cex=cex, inset=0.01, legend=c("Leading edge", "Trailing edge"), lwd=c(2, 2), col=c("dodgerblue3", "darkred")) for(il in 1:2) points(x=temp$text$x[il] - 2*0.7*par("cex")*xinch(par("cin")[1L], warn.log = FALSE), y=temp$text$y[il], pch=21, col="white", bg=c("darkblue", "darkred")[il], cex=c(0.7, 0.7)[il]) } } } par(op) dev.off() } #Figure S2 fig.ClimateSoilSpace(dir.out, figname="Fig_ClimateSoilSpace.pdf") #------------------------ #Climate scenario deltas fig.MonthlyClimateScenarioDelta <- function(dir.fig, figname){ ymin <- c(0, -3, -3) ymax <- c(5, 12, 12) pdf(height=3*length(var.ClimChange), width=5*length(rcps)*ncol(edge), file=file.path(dir.fig, figname)) op <- par(mfcol=c(length(var.ClimChange), length(rcps)*ncol(edge)), mar=c(3, 3, 1, 0), mgp=c(1.5, 1, 0), las=1, cex=1) for(ip in seq_along(rcps)){ iclims <- c(1, grep(rcps[ip], climate.conditions)) for(ie in 1:ncol(edge)){ iedge <- which(edge[, ie] > 0) for(iv in seq_along(var.ClimChange)){ for(ix in 1:2){ for(isc in iclims){ y <- climateDelta[iedge, (iv - 1)*length(mo) + mo, isc] if(ix == 1){ #draw polygons yMin <- apply(y, 2, min) yMax <- apply(y, 2, max) if(isc == 1){ plot(mo, yMin, type="n", ylim=c(ymin[iv], ymax[iv]), xlab=ifelse(iv == length(var.ClimChange), "Months", ""), ylab=ifelse(ip == 1 && ie == 1, label.ClimChange[iv], ""), lwd=2, col=colorLines[isc], bty="l", axes=FALSE) axis(side=1, pos=ymin[iv], at=mo, labels=(iv == length(var.ClimChange))) axis(side=2, pos=1, at=sort(unique(c(axTicks(2), ymin[iv], ymax[iv]))), labels=(ip == 1 && ie == 1)) if(iv == 1) title(main=paste(label.Edges[ie], rcps[ip])) } polygon(x=c(mo, rev(mo)), y=c(yMin, rev(yMax)), col=col2alpha("gray", 0.2), border=NA) } else { #draw lines on top of polygons yMean <- apply(y, 2, mean) lines(mo, yMean, lwd=1, col=colorLines[isc]) } } } } } } par(op) dev.off() } #Figure S4 fig.MonthlyClimateScenarioDelta(dir.out, figname="Fig_ClimateScenariosDeltas.pdf") #------------------------ #PCA of explanatory variables panel.hist <- function(x, ...) { usr <- par("usr"); on.exit(par(usr)) par(usr = c(usr[1:2], 0, 1.5) ) h <- hist(x, plot = FALSE) breaks <- h$breaks; nB <- length(breaks) y <- h$counts; y <- y/max(y) rect(breaks[-nB], 0, breaks[-1], y, col="grey", ...) } panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...) { usr <- par("usr"); on.exit(par(usr)) par(usr = c(0, 1, 0, 1)) r <- abs(cor(x, y)) txt <- format(c(r, 0.123456789), digits=digits)[1] txt <- paste(prefix, txt, sep="") if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) text(0.5, 0.5, txt, cex = cex.cor * r) } fig.doPCA_Current <- function(datExp, labelsExp, chosenVars1=NULL, chosenVars2=NULL, transform=FALSE, title=FALSE, dir.fig, figname){ dat <- datExp[, , 1] temp <- apply(dat, 2, sd, na.rm=FALSE) dat <- dat[, !is.na(temp) & temp > 0] pdf(width=1*dim(dat)[2], height=1*dim(dat)[2], file=file.path(dir.fig, sub(".pdf", "_PairedScatter.pdf", figname))) op <- par(mfcol=c(1, 2)) for(ie in 1:ncol(edge)) pairs(dat[edge[, ie] == 1, ], diag.panel=panel.hist, lower.panel=panel.smooth, upper.panel=panel.cor) #pairwise-scatterplot & correlation coefficients par(op) dev.off() if(transform) { #transform variables if very skewed transform.log <- NULL transform.sqrt <- NULL transform.sqr <- NULL transform.asin <- NULL transform.exp <- NULL if(names(list.final[vl])=="Snow"){ transform.asin <- 1 transform.log <- 2 } if(names(list.final[vl])=="SoilWaterFlow"){ transform.log <- 1:3 } if(names(list.final[vl])=="SoilWaterQuantityAB"){ transform.exp <- c(1, 3) transform.log <- c(12, 14) transform.sqrt <- c(5, 7, 8, 10) } if(names(list.final[vl])=="SoilWaterTimingAB"){ transform.sqrt <- c(1:4) transform.sqr <- c(10) } if(names(list.final[vl])=="Vegetation"){ transform.asin <- 2:3 } if(length(transform.log)>0){ dt[, transform.log] <- log1p(dt[, transform.log]) colnames(dt)[transform.log] <- paste("ln(1+", colnames(dt)[transform.log], ")", sep="") } if(length(transform.asin)>0){ dt[, transform.asin] <- asin(sqrt(dt[, transform.asin])) colnames(dt)[transform.asin] <- paste("asin(sqrt(", colnames(dt)[transform.asin], "))", sep="") } if(length(transform.sqrt)>0){ dt[, transform.sqrt] <- sqrt(dt[, transform.sqrt]+0.5) colnames(dt)[transform.sqrt] <- paste("sqrt(0.5+", colnames(dt)[transform.sqrt], ")", sep="") } if(length(transform.sqr)>0){ dt[, transform.sqr] <- dt[, transform.sqr]^2 colnames(dt)[transform.sqr] <- paste("sqr(", colnames(dt)[transform.sqr], ")", sep="") } if(length(transform.exp)>0){ dt[, transform.exp] <- exp(dt[, transform.exp]) colnames(dt)[transform.exp] <- paste("exp(", colnames(dt)[transform.exp], ")", sep="") } #correlations jpeg(width=10, height=10, units="in", res=300, type="quartz", file=paste(dir.out, "/", dir.out.folder, "/", dir.out.subfolder, "/1_Distributions/", vl, "_Scatterplot_", names(list.final[vl]), "_Transformed.jpg", sep="")) pairs(dt, diag.panel = panel.hist, lower.panel = panel.smooth, upper.panel = panel.cor) #pairwise-scatterplot & correlation coefficients dev.off() } #PCA by a singular value decomposition of the (centered and scaled) data matrix Correlation.biplot <- function(pca){ plot(pca$x[,1]/pca$sdev[1], pca$x[,2]/pca$sdev[2], xlab="PCA 1", ylab="PCA 2", type="p", pch=19, col=c(rep("royalblue", 100), rep("red", 100))[if(ie > 2) rep(TRUE, 100) else edge[, ie] == 1], cex=0.6) arrows(0, 0, pca$rotation[, 1] * pca$sdev[1], pca$rotation[, 2] * pca$sdev[2], length=0.1, angle=20, col={ ctemp <- rep("gray", length(res.pca[[ie]]$center)); if(length(chosenVars1) > 0) { chosenIds1 <- match(chosenVars1, table=names(res.pca[[ie]]$center)); ctemp[chosenIds1] <- "darkgreen"; ctemp }; if(length(chosenVars2) > 0) { chosenIds2 <- match(chosenVars2, table=names(res.pca[[ie]]$center)); ctemp[chosenIds2] <- "purple"; ctemp }; ctemp} ) text(pca$rotation[, 1] * pca$sdev[1] * 1.2, pca$rotation[, 2] * pca$sdev[2] * 1.2, 1:length(res.pca[[ie]]$center), col=ctemp, adj=0, cex=0.6) legend(x="topright", bty="n", legend=c("Leading edge", "Trailing edge", "Overall climate/soil", "Selected variables"), pch=c(19, 19, -1, -1), cex=0.5, lty=c(-1, -1, 1, 1), col=c("royalblue", "red", "darkgreen", "purple")) } pdf(width=3*3, height=3*5, file=file.path(dir.fig, figname)) res.pca <- sum.pca <- list() op <- par(mfcol=c(4, 3)) for(ie in 1:(1+ncol(edge))){ if(ie > ncol(edge)){ edat <- dat } else { edat <- dat[edge[, ie] == 1, ] } res.pca[[ie]] <- prcomp(edat, retx=TRUE, center=TRUE, scale.=TRUE, na.action=na.omit) biplot(res.pca[[ie]], pc.biplot=FALSE) #correlation biplot Correlation.biplot(pca=res.pca[[ie]]) #checks max(res.pca[[ie]]$sdev - sd(res.pca[[ie]]$x)) # sd(principal components) - sd(sample scores): should be close to zero [?] all.equal(sum(res.pca[[ie]]$sdev^2), dim(edat)[2]) # should be equal #summed residuals across variables = q-values q <- apply(scale(edat) - predict(res.pca[[ie]]), MARGIN=1, FUN=sum) plot(sort(q)) #broken-stick method to determin which PC to retain: Jackson, D. A. 1993. Stopping Rules in Principal Components Analysis: A Comparison of Heuristical and Statistical Approaches. Ecology 74:2204-2214. p <- dim(edat)[2] eigen.crit <- rev(cumsum(1/(p:1))) pca.retain <- data.frame(eigen <- res.pca[[ie]]$sdev^2, eigen.crit, retained <- eigen>eigen.crit) #screeplot(res.pca[[ie]]) xtemp <- barplot(eigen[1:(temp1 <- min(10, length(eigen)))], col={cols <- rep("gray", temp1); cols[retained[1:temp1]] <- "red"; cols}) lines(xtemp, eigen.crit[1:temp1], lwd=2) #summary pca.importance <- t(as.matrix(summary(res.pca[[ie]])$importance)) pca.correlation <- t(res.pca[[ie]]$rotation)*res.pca[[ie]]$sdev pca.loading <- t(res.pca[[ie]]$rotation) sum.pca[[ie]] <- data.frame(rownames(pca.importance), pca.importance, pca.retain, pca.correlation, pca.correlation) colnames(sum.pca[[ie]]) <- c("PC", colnames(pca.importance), "Eigenvalues", "Ecrit (broken-stick)", "PCretain (broken-stick)", paste("Cor", colnames(pca.correlation), sep="."), paste("Load", colnames(pca.loading), sep=".")) write.csv(sum.pca[[ie]], file=file.path(dir.fig, sub(".pdf", paste0("_summary", ie, ".csv"), figname))) } par(op) dev.off() pdf(width=3, height=3, file=file.path(dir.fig, sub(".pdf", paste0(ie, "Biplot.pdf"), figname))) op <- par(mar=c(2.5, 2.5, 0.1, 0)+0.1, mgp=c(1.5, 0.6, 0), las=0, xpd=FALSE, cex=1) Correlation.biplot(pca=res.pca[[ie]]) par(op) dev.off() return(list(res.pca, sum.pca)) } #Figure S1, Table S2 pca <- fig.doPCA_Current(datExp=explanatories, labelsExp=label.Explanatory, chosenVars1=var.Explanatory2a, chosenVars2=var.Explanatory2b, title=TRUE, dir.fig=dir.out, figname="Fig_PCA_ExplanatoryVariables.pdf") #------------------------ #Regeneration fig.ResponseBoxplot <- function(rcp, do.change=FALSE, title=FALSE, dir.fig, figname){ h.panel <- 2.5; w.panel <- 5; w.edge <- 0.75; h.edge <- 1.5 pdf(height=h.edge+h.panel*length(rcp)*length(var.Regeneration), width=w.edge+w.panel*ncol(edge), file=file.path(dir.fig, figname)) layout(matrix(1:((1+length(rcp)*length(var.Regeneration))*(1+ncol(edge))), ncol=1+ncol(edge), byrow=TRUE), widths=c(w.edge, rep(w.panel, times=ncol(edge))), heights=c(rep(h.panel, times=length(rcp)*length(var.Regeneration)), h.edge)) op <- par(mar=c(0.5, 0, 1, 0), mgp=c(1.5, 1, 0), las=2, xpd=FALSE, cex=1) if(do.change){ ylim <- c(-1, 1) ylabs <- gsub("\n", " change\n", label.Regeneration) } else { ylim <- c(0, 1) ylabs <- label.Regeneration } do.title <- FALSE letter <- 1 for(ip in seq_along(rcp)){ iclims <- grep(rcp[ip], climate.conditions) if(!do.change) iclims <- c(1, iclims) if(title) do.title <- TRUE for(ir in seq_along(var.Regeneration)){ plot.new() for(ie in rev(1:ncol(edge))){ iedge <- which(edge[, ie] > 0) y <- regenerationResponse[iedge, ir, iclims] if(do.change) y <- y - regenerationResponse[iedge, ir, 1] bdat <- boxplot(y, range=0, notch=TRUE, ylim=ylim, col=colorPointsByEdge[ie, iclims], medcol=colorPointsByEdge2[ie, iclims], names=rep("", length(iclims)), xlab="", ylab="", axes=FALSE) segments(x0=1, y0=bdat$stats[3, 1], x1=length(iclims), y1=bdat$stats[3, 1], col=colorPointsByEdge[ie, 1], lwd=2) axis(side=1, pos=ylim[1], at=seq_along(iclims), labels=if(ip*ir == length(rcp)*length(var.Regeneration)) sapply(strsplit(climate.conditions[iclims], split=".", fixed=TRUE), FUN=function(x) tail(x, n=1)) else FALSE) axis(side=2, pos=0.5, labels=(ie == 2)) if(do.change) lines(x=c(0.5, length(iclims)+0.5), y=c(0, 0), col="red") if(ie == 2) mtext(side=2, line=1.5, text=ylabs[ir], las=3) if(do.title) title(main=paste0(label.Edges[ie], ": ", rcp[ip], ifelse(do.change, " - Current", ""))) mtext(side=3, line=-0.25, text=paste0("(", letters[letter], ")"), las=0, adj=0) letter <- letter + 1 } do.title <- FALSE } } par(op) dev.off() } #Figure S7 fig.ResponseBoxplot(rcp="RCP45", do.change=FALSE, title=FALSE, dir.out, figname="Fig_Boxplot_Regeneration_RCP45.pdf") #Figure 4 fig.ResponseBoxplot(rcp="RCP85", do.change=FALSE, title=FALSE, dir.out, figname="Fig_Boxplot_Regeneration_RCP85.pdf") #2nd sub-section of result section tab.ResponseBoxplot <- function(rcp, dir.tab, tabname){ nclims <- max(sapply(seq_along(rcp), FUN=function(ip) length(c(1, grep(rcp[ip], climate.conditions))))) resp <- paste(rep(temp <- c("Mean", "SD", "Median", "MAD"), times=2), rep(c("Absolute", "Change"), each=length(temp)), sep=".") vars <- dimnames(regenerationResponse)[2][[1]] exps <- c("Edge", "RCP", "Scenario") res <- data.frame(matrix(NA, ncol=length(exps) + length(vars) * length(resp), nrow=nclims, dimnames=list(NULL, c(exps, paste(rep(vars, each=length(resp)), rep(resp, times=length(vars)), sep="_"))))) for(ip in seq_along(rcp)){ iclims <- c(1, grep(rcp[ip], climate.conditions)) res[, 2] <- rcp[ip] res[, 3] <- climate.conditions[iclims] for(ie in 1:ncol(edge)){ iedge <- which(edge[, ie] > 0) res[, 1] <- colnames(edge)[ie] for(isc in seq_along(iclims)){ for(iv in seq_along(vars)){ for(ix in 1:2){# 1 is absolute, 2 is change if(ix == 1){ x <- regenerationResponse[iedge, iv, iclims[isc]] icoff <- 0 } else { x <- x - regenerationResponse[iedge, iv, 1] icoff <- 4 } res[isc, length(exps) + 1 + icoff + (iv - 1) * length(resp)] <- mean(x) res[isc, length(exps) + 2 + icoff + (iv - 1) * length(resp)] <- sd(x) res[isc, length(exps) + 3 + icoff + (iv - 1) * length(resp)] <- median(x) res[isc, length(exps) + 4 + icoff + (iv - 1) * length(resp)] <- mad(x) } } } #Overall meansĀ±SD imean <- grep("Mean.", colnames(res)) res[length(iclims) + 1, imean] <- apply(res[2:nclims, imean], 2, mean) res[length(iclims) + 1, imean + 1] <- apply(res[2:nclims, imean], 2, sd) #Overall medianĀ±MAD imed <- grep("Median.", colnames(res)) res[length(iclims) + 1, imed] <- apply(res[2:nclims, imed], 2, median) res[length(iclims) + 1, imed + 1] <- apply(res[2:nclims, imed], 2, mad) if(ip == 1 && ie == 1){ write.csv(res, row.names=FALSE, file=file.path(dir.tab, tabname)) } else { write.table(res, append=TRUE, sep=",", dec=".", qmethod="double", col.names=FALSE, row.names=FALSE, file=file.path(dir.tab, tabname)) } res <- res[1:length(iclims), ] } } return(0) } tab.ResponseBoxplot(rcp="RCP85", dir.tab=dir.out, tabname="Tab_Regeneration_RCP85.csv") #------------------------ #Regeneration vs gradients plots (leading and trailing edge in same panels) with absolute and change panels in same plot fig.ResponseVSgradients3 <- function(datExp, labelsExp, rcp, line.method="lm", alpha=0.05, title=FALSE, dir.fig, figname){ line.method <- match.arg(line.method, choices=c("lm", "dcor")) varsExp <- dimnames(datExp)[2][[1]] h.panel <- 1.4; w.panel <- 2; w.edge <- 0.5; h.edge <- 0.5 pdf(height=h.edge+h.panel*2*length(rcp)*length(var.Regeneration), width=w.edge+w.panel*length(varsExp), file=file.path(dir.fig, figname)) layout(matrix(1:((1+length(varsExp))*(1+2*length(rcp)*length(var.Regeneration))), nrow=1+2*length(rcp)*length(var.Regeneration), byrow=FALSE), heights=c(rep(h.panel, times=2*length(rcp)*length(var.Regeneration)), h.edge), widths=c(w.edge, rep(w.panel, times=length(varsExp)))) cex <- 0.6 op <- par(mar=c(0.25, 0.5, 1, 0.1), mgp=c(1.5, 0.6, 0), las=0, xpd=FALSE, cex=cex) ylim <- list(c(0, 1), c(-1, 1)) ylabs <- list(label.Regeneration, gsub("\n", " change\n", label.Regeneration)) lcount <- rep(0, 2) letter <- 1 lettersExt <- c(letters, sapply(letters, FUN=function(x) paste0(x, x))) for(i in 1:(1+2*length(rcp)*length(var.Regeneration))) plot.new() for(iv in seq_along(varsExp)){ xlim <- c(min(datExp[, iv, ], na.rm=TRUE), max(datExp[, iv, ], na.rm=TRUE)) #xlimCur <- sapply(1:ncol(edge), FUN=function(ie) c(min(explanatories[which(edge[, ie] > 0), iv, 1], na.rm=TRUE), max(explanatories[which(edge[, ie] > 0), iv, 1], na.rm=TRUE))) xlimCur <- sapply(1:ncol(edge), FUN=function(ie) quantile(datExp[which(edge[, ie] > 0), iv, 1], probs=c(0.05, 0.95), na.rm=TRUE)) #xlimScen <- sapply(1:ncol(edge), FUN=function(ie) c(min(explanatories[which(edge[, ie] > 0), iv, -1], na.rm=TRUE), max(explanatories[which(edge[, ie] > 0), iv, -1], na.rm=TRUE))) #xlimScen <- sapply(1:ncol(edge), FUN=function(ie) quantile(explanatories[which(edge[, ie] > 0), iv, -1], probs=c(0.05, 0.95), na.rm=TRUE)) xlimScen <- sapply(1:ncol(edge), FUN=function(ie) c(min((temp <- sapply(2:dim(datExp)[3], FUN=function(isc) quantile(datExp[which(edge[, ie] > 0), iv, isc], probs=c(0.05, 0.95), na.rm=TRUE)))[1, ]), max(temp[2, ]))) for(ic in seq_along(rcp)){ iclims <- list(c(1, temp <- grep(rcp[ic], climate.conditions)), temp) for(ir in seq_along(var.Regeneration)){ for(ip in 1:2){ #first plot absolute values, second plot changes for(ix in 1:2){ for(isc in rev(iclims[[ip]])){ for(ie in 1:ncol(edge)){ iedge <- which(edge[, ie] > 0) x <- datExp[iedge, iv, isc] y <- list(temp <- regenerationResponse[iedge, ir, isc], temp - regenerationResponse[iedge, ir, 1]) if(ix == 1){#first draw points if(isc == tail(iclims[[ip]], n=1) && ie == 1){#set up plot plot(x, y[[ip]], xlim=xlim, ylim=ylim[[ip]], xlab="", ylab="", pch=19, cex=ifelse(isc == 1, 0.5, 0.1), col=colorPointsByEdge[ie, isc], axes=FALSE, type="n") axis(side=1, pos=ylim[[ip]][1], labels=(ic == length(rcp) && ir == length(var.Regeneration) && ip == 2)) axis(side=2, labels=(iv == 1)) if(ic == length(rcp) && ir == length(var.Regeneration) && ip == 2) mtext(side=1, line=1, text=labelsExp[iv], padj=1, cex=cex) if(iv == 1) mtext(side=2, line=1.5, text=ylabs[[ip]][ir], cex=cex) if(ip == 2) lines(x=par()$usr[1:2], y=c(0, 0), col="black") if(title && ir == 1) title(main=paste0(strsplit(label.Regeneration[ir], split="\n", fixed=TRUE)[[1]][1], ": ", rcp[ic])) } if(isc == tail(iclims[[ip]], n=1)){#add range polygons polygon(x=c(xlimCur[, ie], rev(xlimCur[, ie])), y=c(rep(ylim[[ip]][2] - 0.2*(ylim[[ip]][2] - ylim[[ip]][1]), 2), rev(rep(ylim[[ip]][2], 2))), border=NA, col=col2alpha(colorLinesByEdge[ie, 1], 0.1)) polygon(x=c(xlimScen[, ie], rev(xlimScen[, ie])), y=c(rep(ylim[[ip]][1], 2), rev(rep(ylim[[ip]][1] + 0.2*(ylim[[ip]][2] - ylim[[ip]][1]), 2))), border=NA, col=col2alpha(colorLinesByEdge[ie, 1], 0.1)) } if(isc == tail(iclims[[ip]], n=1)){#set up counter of lines lcount[ie] <- 0 } points(x, y[[ip]], pch=19, cex=ifelse(isc == 1, 0.5, 0.1), col=colorPointsByEdge[ie, isc]) } else {#second draw lines if(sd(y[[ip]]) > 0 && any(ok <- is.finite(x) & is.finite(y[[ip]]))){ if(line.method == "lm"){ m <- lm(y[[ip]] ~ x) p.value <- anova(m)[1, "Pr(>F)"] if(is.finite(p.value) && p.value < alpha / length(iclims[[ip]])){ lines(sort(x[ok]), predict(m, newdata=data.frame(x=sort(x[ok]))), lwd=lwdLines[isc], col=colorLinesByEdge[ie, isc]) if(isc != tail(iclims[[ip]], n=1)) lcount[ie] <- lcount[ie] + 1 } } else if(line.method == "dcor"){ dc <- energy::dcov.test(x[ok], y[[ip]][ok], R=2000) if(dc$p.value < alpha / length(iclims[[ip]])){ lines(stats::lowess(x[ok], y[[ip]][ok], f=2/3, iter=3), lwd=lwdLines[isc], col=colorLinesByEdge[ie, isc]) #text(x=par("xaxp")[2], y=par("yaxp")[1], label=paste0("dcor = ", formatC(dc$estimates[2], format="f", digits=3), "\np < ", formatC(dc$p.value, format="f", digits=3)), adj=c(1,0), cex=0.8) if(isc != tail(iclims[[ip]], n=1)) lcount[ie] <- lcount[ie] + 1 } } } } } } } mtext(side=3, line=-0.25, text=paste0("(", lettersExt[letter], ")"), las=0, adj=0, cex=cex) letter <- letter + 1 mtext(side=3, line=-0.4, text=expression("GCMs: " * phantom("16 (trailing)") * ", " * phantom("16 (leading)")), cex=cex) mtext(side=3, line=-0.4, text=bquote(phantom("GCMs: ") * .(lcount[2]) * " (trailing)" * phantom(", 16 (leading)")), cex=cex, col="red") mtext(side=3, line=-0.4, text=bquote(phantom("GCMs: 16 (trailing), ") * .(lcount[1]) * " (leading)"), cex=cex, col="blue") } } } plot.new() } par(op) dev.off() } #Figure S8 fig.ResponseVSgradients3(datExp=explanatories2a, labelsExp=label.Explanatory2a, rcp="RCP45", line.method="dcor", title=FALSE, dir.fig=dir.out, figname="Fig_Scatter_RegenerationChangeVSgradients3a_RCP45.pdf") fig.ResponseVSgradients3(datExp=explanatories2a, labelsExp=label.Explanatory2a, rcp="RCP85", line.method="dcor", title=FALSE, dir.fig=dir.out, figname="Fig_Scatter_RegenerationChangeVSgradients3a_RCP85.pdf") #Figure S9 fig.ResponseVSgradients3(datExp=explanatories2b, labelsExp=label.Explanatory2b, rcp="RCP45", line.method="dcor", title=FALSE, dir.fig=dir.out, figname="Fig_Scatter_RegenerationChangeVSgradients3b_RCP45.pdf") #Figure 5 fig.ResponseVSgradients3(datExp=explanatories2b, labelsExp=label.Explanatory2b, rcp="RCP85", line.method="dcor", title=FALSE, dir.fig=dir.out, figname="Fig_Scatter_RegenerationChangeVSgradients3b_RCP85.pdf") #------------------------ #Median daily soil water fig.medianDailySoilWater <- function(rcp, soilWater="SWP", do.change=FALSE, title=TRUE, dir.fig, figname){ soilWater <- match.arg(soilWater, c("SWP", "VWC", "SWC")) swUnit <- switch(soilWater, SWP="MPa", VWC="'m'^3 %.% 'm'^-3", SWC="mm") layerN <- length(layerDepths_cm) - 1 h.panel <- 1.4; w.panel <- 2.5; w.edge <- 0.5; h.edge <- 0.4 pdf(height=h.edge+h.panel*layerN, width=w.edge+w.panel*ncol(edge)*length(rcp), file=file.path(dir.fig, figname)) layout(matrix(1:((1+layerN)*(1+ncol(edge)*length(rcp))), ncol=1+ncol(edge)*length(rcp), byrow=TRUE), heights=c(rep(h.panel, times=layerN), h.edge), widths=c(w.edge, rep(w.panel, times=ncol(edge)*length(rcp)))) cex <- 0.7 op <- par(mar=c(0.25, 0, 1, 0), mgp=c(1.5, 0.75, 0), las=0, xpd=FALSE, cex=cex) tblSoilWaterdoy <- paste0("aggregation_doy_", soilWater) if(do.change){ ylim <- switch(soilWater, SWP=c(-4.5, 4.5), VWC=c(-0.1, 0.1), SWC=c(-50, 50)) ylabs <- paste0(bquote(.(soilWater)), "*~change~(", swUnit, ")") } else { ylim <- switch(soilWater, SWP=c(-8, 0), VWC=c(0, 0.35), SWC=c(0, 100)) ylabs <- paste0(bquote(.(soilWater)), "*~(", swUnit, ")") } do.title <- FALSE letter <- 1 for(il in 1:layerN){ plot.new() for(ie in rev(1:ncol(edge))){ if(do.change){ y_current <- get.Table_Scenario(responseName=tblSoilWaterdoy, MeanOrSD="Mean", scenario="Current", whereClause=paste0("Soil_Layer=", shQuote(il), " AND ", var.Edges[ie], "=1")) if(length(itemp <- grep("Soil_Layer", colnames(y_current))) > 0) y_current <- y_current[, -itemp] y_current <- y_current[, doy] } if(title) do.title <- TRUE for(ip in seq_along(rcp)){ iclims <- grep(rcp[ip], climate.conditions) if(!do.change) iclims <- c(1, iclims) for(ix in 1:2){ for(isc in rev(iclims)){ y <- get.Table_Scenario(responseName=tblSoilWaterdoy, MeanOrSD="Mean", scenario=climate.conditions[isc], whereClause=paste0("Soil_Layer=", shQuote(il), " AND ", var.Edges[ie], "=1")) if(length(itemp <- grep("Soil_Layer", colnames(y))) > 0) y <- y[, -itemp] y <- y[, doy] if(do.change) y <- y - y_current if(nrow(y) > 0){ if(ix == 1){ #draw polygons yMin <- apply(y, 2, min) yMax <- apply(y, 2, max) if(isc == max(iclims)){ plot(doy, yMin, type="n", xlim=c(min(doy), max(doy)), ylim=ylim, xlab="", ylab="", col=colorLines[isc], bty="l", axes=FALSE) axis(side=1, pos=par()$usr[3], at=c(doysAtFirstOfMonths, 365), labels=if(il == layerN) c(month.abb, "") else FALSE) axis(side=2, pos=1, at=sort(unique(c(axTicks(2), ylim))), labels=(ip == 1 && ie == 2)) if(il == layerN) mtext(side=1, line=2, text="Day of year", cex=cex) #if(ip == 1 && ie == 2) mtext(side=2, line=1.5, text=paste0("Soil depth: ", layerDepths_cm[il], "-", layerDepths_cm[il+1], " cm\n", ylabs), las=3, cex=cex) if(ip == 1 && ie == 2){ mtext(side=2, line=2, text=parse(text=ylabs), las=3, cex=cex) mtext(side=2, line=1, text=paste0("at ", layerDepths_cm[il], "-", layerDepths_cm[il+1], " cm"), las=3, cex=cex) } #if(ip == 1 && ie == 2) mtext(side=2, line=1.5, text=paste0(ylabs, "\nat ", layerDepths_cm[il], "-", layerDepths_cm[il+1], " cm"), las=3, cex=cex) if(il == 1 && do.title) title(main=paste0(label.Edges[ie], ": ", rcp[ip], ifelse(do.change, " - Current", "")), cex=cex) } polygon(x=c(doy, rev(doy)), y=c(yMin, rev(yMax)), col=col2alpha(if(length(rcp) > 1) colorPointsByEdge[ie, isc] else "gray", 1/(0.6*length(iclims))), border=NA) } else { #draw lines on top of polygons yMean <- apply(y, 2, median) lines(doy, yMean, lwd=lwdLines[isc], col=colorLines[isc]) } } } } if(do.change) lines(x=c(min(doy), max(doy)), y=c(0, 0), col="red") mtext(side=3, line=-0.5, text=paste0("(", letters[letter], ")"), las=0, adj=0.05, cex=cex) letter <- letter + 1 } do.title <- FALSE } } par(op) dev.off() } sw <- "SWP" #Figure S5 fig.medianDailySoilWater(rcp=rcps, soilWater=sw, do.change=FALSE, title=TRUE, dir.fig=dir.out, figname=paste0("Fig_medianDaily", sw, ".pdf")) sw <- "VWC" #Figure S6 fig.medianDailySoilWater(rcp="RCP45", soilWater=sw, do.change=TRUE, title=FALSE, dir.fig=dir.out, figname=paste0("Fig_medianDaily", sw, "ChangeRCP45.pdf")) #Figure 3 fig.medianDailySoilWater(rcp="RCP85", soilWater=sw, do.change=TRUE, title=FALSE, dir.fig=dir.out, figname=paste0("Fig_medianDaily", sw, "ChangeRCP85.pdf"))