library(TPD) library(shape) library(mgcv) library(colortools) library(RColorBrewer) library(scico) library(viridis) source("Aux_Functions.R") source("Aux_PlotVariables.R") PCAImputedALL <- readRDS("dataPrepared/PCAImputedALL.rds") targetGroupName <- names(PCAImputedALL) colGradient <- c("white", "yellow", "red") gradientColorsF <- colorRampPalette(colGradient, space = "Lab") ncolors <- 1000 ColorRamp <- rev(gradientColorsF(ncolors)) contourLevels <- c(0.5, 0.99) summaryResults <- readRDS("dataPrepared/NullModelsFRichnessSummaryResults.rds") RawResultsNull <- readRDS("dataPrepared/NullModelsFRichnessRawResults.rds") MAEResults <- readRDS("dataPrepared/observedMAE.rds") pdf("Figures/Figure_FunctionalShifts.pdf", width = 33, height = 20, pointsize=0.85*ptSize) layout(matrix(c(1, 2, 3, 7, 4, 5, 6, 7), nrow=2, byrow=T), widths=c(0.3, 0.3, 0.3, 0.1)) par(mar = c(5, 5, 2, 0.3), oma = c(1, 1, 1, 1), mgp = c(2.2, 0.1, 0), cex.axis = 1.25, cex.lab = 1.5) for(gr in 1:length(PCAImputedALL)){ cat(paste("\n Plotting: ", targetGroupName[gr], "\n")) data <- PCAImputedALL[[gr]] dimensions <- PCAImputedALL[[gr]]$dimensions traitsUSE <- PCAImputedALL[[gr]]$traitsUse IUCNStatus <- which(!is.na(PCAImputedALL[[gr]]$Threat)) dataAnalysis <- data.frame(traitsUSE[IUCNStatus,], Threat = PCAImputedALL[[gr]]$Threat[IUCNStatus]) commNames <- c("ALL", "IUCN", "WithoutThreatened") commMatrix <- matrix(0, nrow = length(commNames), ncol = nrow(dataAnalysis), dimnames = list(commNames, rownames(dataAnalysis))) commMatrix["ALL", ] <- 1 commMatrix["IUCN", rownames(dataAnalysis)] <- 1 commMatrix["WithoutThreatened", rownames(subset(dataAnalysis, dataAnalysis$Threat == 0))] <- 1 TPDsAux <- readRDS(paste0("dataPrepared/TPDs/TPDs_2D_highDef",targetGroupName[gr], ".rds")) TPDcAux <- TPDc(TPDs = TPDsAux, sampUnit = commMatrix) xlab <- paste0("PC1 (", round(100 * PCAImputedALL[[gr]]$variance[1], 2), "%)") ylab <- paste0("PC2 (", round(100 * PCAImputedALL[[gr]]$variance[2], 2), "%)") imageMatALL <-imageTPD(TPDcAux, thresholdPlot = 0.99)[,,"ALL"] imageMatIUCN <-imageTPD(TPDcAux, thresholdPlot = 0.99)[,,"IUCN"] imageMatRemain <-imageTPD(TPDcAux, thresholdPlot = 0.99)[,,"WithoutThreatened"] imageMatIUCN1 <-imageTPD(TPDcAux, thresholdPlot = 1)[,,"IUCN"] imageMatRemain1 <-imageTPD(TPDcAux, thresholdPlot = 1)[,,"WithoutThreatened"] imageMatChange <- imageMatIUCN - imageMatRemain imageMatLostSpace <- replace(imageMatIUCN, imageMatIUCN>0, NA) imageMatLostSpace <- replace(imageMatLostSpace, imageMatIUCN>0 & is.na(imageMatRemain), 1) Comp.1Vec <- unique(TPDcAux$data$evaluation_grid[,1]) Comp.2Vec <- unique(TPDcAux$data$evaluation_grid[,2]) limX <- limXlist[[gr]] limY <- limYlist[[gr]] ncol <- 1000 colorsPick <- "vik" ColorRamp <- rev(scico(n=ncol, palette=colorsPick)) nHalf <- sum(!is.na(imageMatChange))/2 Min <- -0.36 Max <- 0.29 Thresh <- 0 ## Make vector of colors for values below threshold rc1 <- colorRampPalette(colors = c(ColorRamp[1:500]), space="Lab")(nHalf) ## Make vector of colors for values above threshold rc2 <- colorRampPalette(colors = c(ColorRamp[501:1000]), space="Lab")(nHalf) rampcols <- c(rc1, rc2) rb1 <- seq(Min, Thresh, length.out=nHalf+1) rb2 <- seq(Thresh, Max, length.out=nHalf+1)[-1] rampbreaks <- c(rb1, rb2) image(x=Comp.1Vec, y=Comp.2Vec, z=imageMatChange, xlim=limX, ylim=limY, col=rampcols, breaks=rampbreaks, xaxs="r", yaxs="r", axes=F, asp=1, xlab=paste0("PC1 (", round(100 * PCAImputedALL[[gr]]$variance[1], 2), "%)"), ylab=paste0("PC2 (", round(100 * PCAImputedALL[[gr]]$variance[2], 2), "%)")) contIUCN99<-contourLines(x = Comp.1Vec, y = Comp.2Vec, z = imageMatIUCN1, levels= c(0.99)) contRemain99<-contourLines(x = Comp.1Vec, y = Comp.2Vec, z = imageMatRemain1, levels= c(0.99)) for(fig in 1:length(contIUCN99)){ polygon(x=c(contIUCN99[[fig]]$x), y=c(contIUCN99[[fig]]$y), col=1) } image(x=Comp.1Vec, y=Comp.2Vec, z=imageMatChange, xlim=limX, ylim=limY, col=rampcols, breaks=rampbreaks, add=T) box(which="plot") axis(1,tcl=0.3,lwd=0.8) axis(2, las=1, tcl=0.3,lwd=0.8) mtext(text=targetGroupName[gr],side=3, cex = cexMain, line=0.5) mtext(text=letters[gr], font=2, side=3, cex = cexMain, line=0.5, adj=0) percLoss <- round(100*summaryResults[gr, "%lossWithoutThreat"],2) MAE <- round(100*MAEResults[[gr]]["ObsMAE"], 2) percLossSignificance <- c("ns", "***", ".", "ns", "ns", "*") percMAESignificance <- c("***", "***", "***", "***", "***", "ns") legendPos <- ifelse(gr<6, "bottomleft", "topleft") legend(legendPos, legend = c(paste0("Lost space = ", percLoss, "% (", percLossSignificance[gr],")"), paste0("Mean absolute quantile change = ", MAE, "% (", percMAESignificance[gr], ")")), bty="n", cex=1.5) } ###GRADIENT LEGEND: maxRange <- max(abs(c(Min, Max))) par(mar=c(1,1.5,2,1)) ##Para cada grafica: A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot legend_image <- as.raster(matrix(rev(ColorRamp), ncol=1)) plot(c(0,2),c(-maxRange, maxRange),type = 'n', axes = F,xlab = '', ylab = '') mtext('Quantile\nchange', side=3, cex = cexMain, line=-0.5) text(x=1.5, y = c(-maxRange, 0, maxRange), labels = c(-100*maxRange, 0, 100*maxRange), cex=1.25) rasterImage(legend_image, xleft=0, ybottom=-maxRange, xright=1, ytop=maxRange) heights <- seq(0, 1, by=0.1) lines(x=c(0,1), y=c(0,0), lwd=1.5) dev.off()