######### Read before starting analysis #############
##### We do not own the source data, and datasets used to support the findings of this study are available for download by request from their respective providers. ####
##### Species assessments (habitat preference, elevational range) and range map polygons can be accessed and requested from the IUCN Red List of Threatened Species (www.iucnredlist.org/) and BirdLife International #####
##### Protected area maps can be requested from the World Database on Protected Areas (www.protectedplanet.net). ####### 
##### Records of protected area downgrading, downsizing or degazettement events were obtained from PADDDtracker (www.PADDDtracker.org) Version 2.1. ######


######### Step 1: Loading required packages #########
library(sf)
library(terra)
######### Step 2: Loading data from IUCN #############

######### Step 2a: Species level assessment data #######
######### Can be extracted by searching for species in https://www.iucnredlist.org/ ######
######### For example: Balebreviceps hillmani #####
######### Download "Search Results" and "Range data - Polygons" #####
######### Open folder from "Search Results" #########

setwd("SearchResults_Folder") #### load search results folder 
assessment <- read.csv("assessments.csv") #### load assessment ###
elevation <- read.csv("all_other_fields.csv")[,c(2,10:11)] #### load elevation data ####
elevation[is.na(elevation[,2]),2] <- -999 #### just in case data are coded as NA #####
elevation[is.na(elevation[,3]),3] <- 999999 #### just in case data are coded as NA #####
habitats <- read.csv("habitats.csv") ### loading habtiat preferences ####

legend <- read.csv("IUCN_mapping_legend.csv") #### We coded a translator to match habitat preference data from "habitat" object to the map from Jung et al. 2020;  "IUCN_mapping_legend.csv" is uploaded together with this code ###

######### Step 2b: Species level polygon data ############
sppShp <- st_read("Balebreviceps_hillmani.shp") #### Load downloaded shapefile

######### Step 3: Load habitat map from Jung et al. 2020, elevation map from Fick et al. 2017, and land use change map from Buchhorn et al. 2020 #### 
habitatMap <- rast("habitatclassification.tif")
elevationMap <- rast("elevation.tif")
changeMap <- rast("change.tif") ### map has been reclassified such that all unchanged areas are zeros and all anthropogenic expansion areas are NA ###

######### Step 4: Converting sf file to vector  #######
x1 <- vect(sppShp)

######### Step 5: Habitat mapping ##
hab1 <- crop(habitatMap,ext(x1)) ### cropping to area of focus (speeds up process)
hx1 <- mask(hab1,x1) #### removing all areas outside range polygon

spphab <- habitats[,"code"] #### which habitats do the species use ? ####
rcl <- as.matrix(data.frame(from=legendConv[legendConv [,1] %in% spphab,2],to=0))  #### Creating a matrix of suitable habitat for reclassification ####
hx2 <- classify(hx1,rcl=rcl,othersNA=T) ##### reclassification: all suitable habitats become zero and unsuitable becomes NA ######
	
######### Step 6: Elevation mapping ###
ele1 <- crop(elevationMap,ext(x1))
ele2 <- mask(ele1,x1)

elee <-  elevation[,2:3]
elee1 <- elee[order(elee)]
ele3 <- classify(ele2,rcl=as.matrix(cbind(elee1,0)),othersNA=T) ### reclassifying all values within suitable elevation range to become zero and outside the range to be NA ####

######### Step 7: AOH = habitat + elevation ####
aoh <- sum(ele3,hx2) ### combining elevational maps and habitat maps ##

######### Step 8: Updating to 2019 ####
change <- crop(changeMap,ext(aoh))
aoh1 <- sum(aoh,change)

######### Step 9: Calculating area of habitat #####
size <- cellSize(aoh1,unit="ha")
global(size,fun="sum",na.rm=T) #### area of habitat value in hectares 

######### Step 10: Loading protected area and PADDD data ####
pa <- rast("protectedAreasConsolidated.tif") ### we rasterize the consolidated WDPA shapefile data using the fasterize package ####
paddd <- rast("PADDDConsolidated.tif") ### we rasterize the PADDD data from PADDDtracker.org after the data was matched to data from WDPA. Rasterized data reflects PAs that are afftected by 1) downgrading, 2) downsizing or 3) degazettement events alone, 4) multiple types of events that have not been reversed, and 5) reversed downgrading events, 6) reversed downsizing events, 7) reversed degazettement events, and 8) reversed combination of event types  ###
pa1 <- crop(pa,ext(aoh1))
paddd1 <- crop(paddd,ext(aoh1))
combinedPAandPADDD <- mosaic(pa1,paddd1) ### Each category (including unaffected PAs) are coded separately ####
######### Step 11: Calculating protected area and PADDD area per species ####
extracted <- na.omit(data.frame(codes=values(combinedPAandPADDD),size=values(size)))
as.data.frame(xtabs(extracted[,2]~extracted[,1])) ### breakdown of habitat area in each category

######### Step 12: Loading future human land-use change expansion ######
### Data is obtained from the LUH2 project ###
### We isolate maps of cropland and urban expansion (to 2050) and overlay these maps with existing protected areas ###
### These maps were resampled to a 1km resolution matching the rest of our data ##
### SSP maps directly indicate cropland and urban expansion within existing PAs with zeros indicating areas with expansion and NAs indicating areas with no expansion ###
ssp1 <- rast("ssp1.tif") 
ssp2 <- rast("ssp2.tif")
ssp3 <- rast("ssp3.tif")
ssp5 <- rast("ssp5.tif")

######### Step 13: Calculating species area of habitat affected by future human land-use change expansion ######
ssp1a <- sum(ssp1,size)
global(ssp1a,fun="sum",na.rm=T) 
ssp2a <- sum(ssp2,size)
global(ssp2a,fun="sum",na.rm=T) 
ssp3a <- sum(ssp3,size)
global(ssp3a,fun="sum",na.rm=T) 
ssp5a <- sum(ssp5,size)
global(ssp5a,fun="sum",na.rm=T) 

####### Results are then tabulated accordingly ###