############################################################################################################################
# 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, ggplot2, ENMTools, scales, plotly, 
               ggfortify, FactoMineR, factoextra, classInt, dismo, XML,
               maps, scico, dplyr, tidyverse, stars, ggExtra, bivariatemaps)

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

AVGResProb = raster::mask(raster::raster("AvgOutputs/AvgResearchProb_Carvalho_et_al_2023.tif"),SA)
plot(AVGResProb, col = c("lightyellow","#BEFDA5","#84B44C","#5A7C2A", "#354819"), axes = F)
plot(SA,interior=T,add=T)

PCAClim = raster("GlobalChangeLayers/PCA12_DeltaClimate_1km.tiff") # set diretory

## Open Lapola et al 2023's BAU carbon layer
Carbon = raster("GlobalChangeLayers/Degradation/bau_carbon_final_WGS84.tif")

# Resample to match other layers and mask to the Study Area
Carbon = raster::mask(raster::resample(Carbon, AVGResProb, method = "ngb"), SA)
npixCarbon = nrow(na.omit(as.data.frame(Carbon$bau_carbon_final_WGS84)))
plot(Carbon)

## Splitting degrad/defor in two layers > Degradation and Deforestation
## For more information refer to Lapola et al 2023
### Degradation
Degrad = Carbon
length(Degrad@data@values[Degrad@data@values == 0])
Degrad@data@values[Degrad@data@values == -3] = 0
Degrad@data@values[Degrad@data@values == -1] = NA # Current defor
Degrad@data@values[Degrad@data@values == -2] = NA # Future defor
max(na.omit(Degrad@data@values)) # 600

## Reescaling to 0-1
Degrad@data@values = scales::rescale(abs(Degrad@data@values)) # Reescale
max(na.omit(Degrad@data@values)) # 1
min(na.omit(Degrad@data@values)) # 0
# New resample to match other layers and mask to the Study Area
Degrad = raster::mask(raster::resample(Degrad, AVGResProb, method = "ngb"),SA)
plot(is.na(Degrad), add=T, alpha=0.4)
plot(is.na(Degrad), color = "grey", add = T)

## Saving the final Degradation layer
writeRaster(Degrad, "GlobalChangeLayers/Loss_Degrad_BAU_1km.tiff", overwrite=T)

# ## Deforestation
Defor = Carbon
Defor@data@values[Defor@data@values == -3] = NA # No data or no loss?
Defor@data@values[Defor@data@values >= 0] = NA # No data or no loss?
Defor@data@values[Defor@data@values == -1] = 1 # Current  
Defor@data@values[Defor@data@values == -2] = 2 # Future
#Defor@data@values = rescale(abs(Defor@data@values)) # Reescale
max(na.omit(Defor@data@values)) # 2
min(na.omit(Defor@data@values)) # 1
Defor = raster::mask(raster::resample(Defor, AVGResProb, method = "ngb"),SA)
plot(Defor, col = c("grey","black"))
writeRaster(Defor, "GlobalChangeLayers/Current_Future_Defor_BAU_1km.tiff", overwrite=T)

npixCarbon == nrow(na.omit(as.data.frame(Degrad)))+nrow(na.omit(as.data.frame(Defor)))

### Defining Deforestation's NA mask 
DefNA = na.omit(Defor) 
plot(DefNA)
## masking "Research probability" layer for places with data
### this step is necessary to match layers for bivarplot
AVGRP_deg = raster::mask(AVGResProb, DefNA, maskvalue = T)
plot(AVGRP_deg)

## Listing quartiles for plotting
# Research Prob
RPqt = data.frame(quantile(AVGResProb@data@values, na.rm = T)) #Clima
RPDeg_qt = data.frame(quantile(AVGRP_deg@data@values, na.rm = T)) #Degrad
#Climate
Climqt = data.frame(quantile(PCAClim, na.rm = T))
# Degrad
Degqt = data.frame(quantile(Degrad@data@values, na.rm = T))
# Defor
Defqt = data.frame(quantile(Defor@data@values, na.rm = T))

## Applying bivarmap palette and saving legends
jpeg(filename = "Figures/Figure4_Climate_legend.jpeg")
col.matrix.clim <- colmat(nquantiles= 4, 
                          upperleft= "#5B244E", # Vanimo 
                          upperright= "#1A1D10", 
                          bottomleft= "lightyellow", 
                          bottomright= "#5A7C2A",
                          xlab="Research Probability",
                          ylab="Delta Climate Change")
axis(1, at = c(0,.25, .5, .75,1), labels = c("",round(RPqt[2:4,1],2),""))#X
axis(2, at = c(0,.25, .5, .75,1), labels = c("",round(Climqt[2:4,1],2),""))#Y
dev.off()

jpeg(filename = "Figures/Figure4_Carbon_legend.jpeg")
col.matrix.carb <- colmat(nquantiles= 4, 
                          upperleft= "#5B244E", # Vanimo 
                          upperright= "#1A1D10", 
                          bottomleft= "lightyellow", 
                          bottomright= "#5A7C2A",
                          xlab="Research Probability",
                          ylab="Carbon Loss by Degradation")
axis(1, at = c(0,.25, .5, .75,1), labels = c("",round(RPDeg_qt[2:4,1],2),""))#X
axis(2, at = c(0,.25, .5, .75,1), labels = c("",round(Degqt[2:4,1],2),""))#Y
dev.off()

### Calculating plotting values for bivarplot
## Climate
bivmap     <- bivariate.map(AVGResProb, PCAClim, colormatrix=col.matrix.clim, 
                            nquantiles=4)
### Carbon
bivmapCarb <- bivariate.map(AVGResProb, Degrad, colormatrix=col.matrix.carb, 
                            nquantiles=4)

## Plotting BIVARPLOT
### Climate (Fig 4 - panel A)
jpeg(filename = "Figures/Figure4_BivarPlot_Climate.jpeg", width = 480, height = 480,
     units = "px", pointsize = 12)
raster::plot(raster::mask(bivmap,SA),frame.plot=F,axes=F,box=F,add=F,legend=F,
             col=as.vector(col.matrix.clim))
raster::plot(SA, interior=F, add=T)
dev.off()

### Carbon loss (Degradation and Deforestation) (Fig 4 - panel B)
jpeg(filename = "Figures/Figure4_BivarPlot_CarbonLoss.jpeg", width = 480, height = 480,
     units = "px", pointsize = 12)
raster::plot(raster::mask(bivmapCarb,SA),frame.plot=F,axes=F,box=F,add=F,legend=F,
             col=as.vector(col.matrix.carb))
raster::plot(raster::mask(Defor,SA), col = c("grey", "black"),add=T, legend=F)
raster::plot(SA,interior=T,add=T)
dev.off()

## Joining previous layers for scatterplots with histograms and boxplot (Fig 4 - panels C to E)
data = raster::as.data.frame(mask(raster::brick(AVGResProb,PCAClim,Degrad,bivmap,
                                                bivmapCarb,Defor),SA)) 
names(data)= c("Research_prob", "Delta_Clim", "Degrad", "colclim","colcarb","defor")

## Taking a sample to plot
sample = data[sample(nrow(data), 35000), ]
sample = na.omit(sample)
sample = sample[order(sample$Research_prob, decreasing=TRUE), ]


## Scatterplotting (panels C and D)

## Climate (panel C)
scplot = 
  ggplot(sample, aes(Research_prob, Delta_Clim, col = as.factor(colclim)))+
  geom_point(alpha = 0.5, show.legend = F)+
  scale_color_manual(values = col.matrix.clim[2:5,2:5])+
  xlab("Research probability")+
  ylab("Delta Climate Change")+
  xlim(c(0,1))+ylim(c(0,1))+
  geom_smooth(method="loess", alpha=0.5,  
              linetype = "dashed", color="black", fill=NA, se = F)
ggplotly(scplot)

## Interactive plot 
scplot2 = scplot+
  theme(axis.text.x = element_text(color = "black", size = 18), 
        axis.title.x = element_text(color = "black", size = 18)) +
  theme(axis.text.y = element_text(color = "black", size = 18), 
        axis.title.y = element_text(color = "black", size = 18))+
  theme_classic()
scplot2

# Adding lateral histograms
scplotextra = ggExtra::ggMarginal(scplot2, 
                                  type="density", margins="both", fill="gray85", size=10)
scplotextra
ggsave(scplotextra, filename = 'Figures/Figure4_climate_scatterplot.jpeg', device = "jpeg") # you can set size and other parameters

## Degradation (Panel D)
scplot_degrad = 
  ggplot(sample, aes(Research_prob, Degrad, col = as.factor(colcarb)))+
  geom_point(alpha = 0.5, show.legend = F)+
  scale_color_manual(values = col.matrix.carb[2:5,2:5])+
  xlab("Research probability")+
  ylab("Carbon Loss by Degradation")+
  xlim(c(0,1))+ylim(c(0,1))+
  geom_smooth(method="loess", alpha=0.5,  
              linetype = "dashed", color="black", fill=NA, se = F)
ggplotly(scplot_degrad)

## Interactive plot 
scplot_degrad2 = scplot_degrad+
  theme(axis.text.x = element_text(color = "black", size = 18), 
        axis.title.x = element_text(color = "black", size = 18)) +
  theme(axis.text.y = element_text(color = "black", size = 18), 
        axis.title.y = element_text(color = "black", size = 18))+
  theme_classic()
scplot_degrad2

# Adding lateral histograms
scplotextra_d = ggExtra::ggMarginal(scplot_degrad2, 
                                    type="density", margins="both", fill="gray85", size=10)
scplotextra_d
ggsave(scplotextra_d, filename = 'Figures/Figure4_carbon_loss_scatterplot.jpeg', device = "jpeg") # you can set size and other parameters

## Deforestation boxplots (Panel E)
dt_defor = na.omit(as.data.frame(brick(AVGResProb,Defor)))
names(data)
dt_defor = na.omit(data[,c(1,6)])
names(dt_defor)= c("Research_prob", "Defor")
sample_defor = dt_defor[sample(nrow(dt_defor), 100000), ]
sample_defor$Defor[sample_defor$Defor == 1] = "Current"
sample_defor$Defor[sample_defor$Defor == 2] = "Future"

bplot_defor = 
ggplot(sample_defor, aes(y = Research_prob, 
                         x = Defor, fill=Defor))+
  geom_violin(alpha = 0.8, show.legend = F)+
  geom_boxplot(show.legend = F, width=0.08, col = c("black", "lightgrey"))+
  ylab("Research probability")+
  xlab("Deforestation")+
  scale_fill_manual(values = c("grey","black"))+
  theme(axis.text.x = element_text(color = "black", size = 18), 
        axis.title.x = element_text(color = "black", size = 18)) +
  theme(axis.text.y = element_text(color = "black", size = 18), 
        axis.title.y = element_text(color = "black", size = 18))+
  theme_classic()
bplot_defor
ggsave(bplot_defor, filename = 'Figures/Figure4_boxplot_deforestation.jpeg', device = "jpeg") # you can set size and other parameters