############################################################################################################################
# SUPPORTING INFORMATION TO

# Title: Pervasive gaps in Amazonian ecological research
# Authors: Raquel L. Carvalho1,2, Angelica Resende1,2, Jos Barlow3, Felipe França4, Mario R. Moura5,6, et al.
# Journal: Current Biology
# 1 Universidade de São Paulo, São Paulo, SP, Brazil
# 2 Empresa Brasileira de Pesquisa Agropecuária, Amazônia Oriental, Belém, PA, Brazil
# 3 Lancaster University, Lancaster, UK
# 4 University of Bristol, Bristol, UK
# 5 Departamento de Biologia Animal, Universidade Estadual de Campinas, Campinas, Brazil
# 6 Departamento de Ciências Biológicas, Universidade Federal da Paraíba, Areia, Brazil
# *Corresponding author for this script: Angelica Resende, gel.florestal@gmail.com

## Open packages
install.packages("pacman")
pacman::p_load(raster, sp, rgdal, rlang, ENMTools, scales,  
               ggfortify, FactoMineR, factoextra, classInt, dismo, XML,
               maps, scico, dplyr, tidyverse, stars, ggExtra, bivariatemaps)

## Open data
## Open Study area shapefile
SA = raster::shapefile("Shapefiles/Bullock_mask_polygon.shp")
plot(SA)

## CLIMATE
## Open IPCC layers
dir  = dir(path = "GlobalChangeLayers/IPCC_layers/")
climate = lapply(paste0("GlobalChangeLayers/IPCC_layers/",dir), raster)
clim_stack = raster::brick(unlist(sapply(climate, crop, y = extent(SA)+1)))
cstack = raster::as.data.frame(clim_stack)

## Data preparation
scaled = cstack # not centre, absolute values only, 0-1 scale
for(i in 1:ncol(scaled)){
  scaled[,i] = scales::rescale(abs(cstack[,i]))}

#### PCA
pca = prcomp(scaled, center = F, scale. = F)
autoplot(pca, variance_percentage = T, xlim = c(-0.06,0.01))

# Preparing Axis to use only the two first
idx = which(!is.na(scaled))
ncomp <- ncol(pca$rotation) # All principal components
r.pca <- clim_stack[[1:ncomp]]
names(r.pca) = c("PC1", "PC2","PC3","PC4","PC5","PC6","PC7","PC8",
                 "PC9","PC10","PC11","PC12","PC13")
for(i in 1:ncomp) {r.pca[[i]][idx] <- pca$x[,i] } 
plot(mask(r.pca,SA))

# First and second components
PCA12 = r.pca[[1]] + r.pca[[2]]
PCA12@data@values = scales::rescale(abs(PCA12@data@values))
plot(PCA12)

## Resample the climate layer to match Research Prob
PCAClim = raster::mask(raster::resample(PCA12, AVGResProb, method = 'bilinear'),SA)
plot(PCAClim)
plot(SA, add=T)

# Saving the layer
writeRaster(PCAClim, "AvgOutputs/PCA12_DeltaClimate_1km.tiff", overwrite=T)