In this code we will account for changes in land cover in the range of each mammal species with distribution area available in the IUNC database.

We will use The Long term Global Land Change of: Song, X.P., Hansen, M.C., Stehman, S.V., Potapov, P.V., Tyukavina, A., Vermote, E.F., & Townshend, J.R. (2018). Global land change 1982-2016. Nature 560, 639–643.

The data can be downloaded here

These data consist of three different characterizations of land cover changes: Tree canopy cover, short vegetation cover and bare ground cover over the past 35 years (1982-2016). We will analyze the three data separately. Pixel values (-100 to 100) represent net percent change over the 35-year period. Negative values represent loss; positive values represent gain; 0 represents no change.

library(raster)
library(rgdal)
library(rgeos)
library(tidyverse)
library(dplyr)

setwd("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_data")

First, let’s work with Short vegetation raster

ras<-raster("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_data/Short_vegetation_change_1982-2016.tif")


# Open mammals ranges:
fold<- list.files("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/IUCN_range_data")
orders<-unlist(str_split(fold, "_"))
orders<- split(orders, 1:2)$`1`

With a loop we select one of the 23 orders (with the polygons contained in a .shp file) and calculate all cover changes for each species.

for (j in 1:23) {
  o<- orders[j]
  
  setwd("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/IUCN_range_data")
  sha<-readOGR(paste(fold[j], "data_0.shp", sep="/", collapse=NULL))
  
  #A list with the species names of the j order:
  spp<-unique(sha$BINOMIAL)
  
  #A table to store the results:
  change<-data.frame(matrix(nrow=length(spp), ncol=4))
  colnames(change)<-c("species", "Short_vegetation_change", "n_pixels", "change_perc")
  
  
  for(i in 1:length(spp)){
    try({
      x <- sha[sha$BINOMIAL == spp[i], ] #select the specie
      ras_crop<-crop(ras,x)              #crop the raster less accurately
      Env_shape1<-mask(ras_crop, x)      #cut the raster more accurately
      
      s<-sum(as.matrix(Env_shape1), na.rm = T) #the sum of all values in the raster
      change$Short_vegetation_change[i]<-s
      nc<-ncell(Env_shape1) - freq(Env_shape1, value=NA) #counting the number of cells with data
      change$species[i]<- spp[i] #store the species name
      change$n_pixels[i]<-nc
      change$change_perc[i]<- (s*(100/nc))/100}) #store the percentage of change
    
    n<-paste(orders[j],"Short_vegetation_change.csv", sep="_") 
    write.csv(change, paste("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_change_tables", n, sep="/"))}} #saving the results.

Now, we do the same thing with other two cover data type: Tree canopy cover and bare ground cover.

#### Bare ground change raster
rm(list=ls())

ras<-raster("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_data/Bare_ground_change_1982-2016.tif")
fold<- list.files("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/IUCN_range_data")
orders<-unlist(str_split(fold, "_"))
orders<- split(orders, 1:2)$`1`

for (j in 1:23) {
o<- orders[j]

setwd("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/IUCN_range_data")
sha<-readOGR(paste(fold[j], "data_0.shp", sep="/", collapse=NULL))

spp<-unique(sha$BINOMIAL)

change<-data.frame(matrix(nrow=length(spp), ncol=4))
colnames(change)<-c("species", "Bare_ground_change", "n_pixels", "change_perc")


for(i in 1:length(spp)){
try({
x <- sha[sha$BINOMIAL == spp[i], ]
ras_crop<-crop(ras,x)               
Env_shape1<-mask(ras_crop, x)

s<-sum(as.matrix(Env_shape1), na.rm = T)
change$Bare_ground_change[i]<-s
nc<-ncell(Env_shape1) - freq(Env_shape1, value=NA)
change$species[i]<- spp[i]
change$n_pixels[i]<-nc
change$change_perc[i]<- (s*(100/nc))/100})

n<-paste(orders[j],"Bare_ground_change.csv", sep="_")
write.csv(change, paste("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_change_tables", n, sep="/"))}}
#### Tree canopy change raster
rm(list=ls())

ras<-raster("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_data/Tree_canopy_change_1982-2016.tif")


fold<- list.files("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/IUCN_range_data")
orders<-unlist(str_split(fold, "_"))
orders<- split(orders, 1:2)$`1`


for (j in 1:23) {
  o<- orders[j]
  
  setwd("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/IUCN_range_data")
  sha<-readOGR(paste(fold[j], "data_0.shp", sep="/", collapse=NULL))
  
  spp<-unique(sha$BINOMIAL)
  
  change<-data.frame(matrix(nrow=length(spp), ncol=4))
  colnames(change)<-c("species", "Tree_canopy_change", "n_pixels", "change_perc")
  
  
  for(i in 1:length(spp)){
    try({
      x <- sha[sha$BINOMIAL == spp[i], ]
      ras_crop<-crop(ras,x)               
      Env_shape1<-mask(ras_crop, x)
      
      s<-sum(as.matrix(Env_shape1), na.rm = T)
      change$Tree_canopy_change[i]<-s
      nc<-ncell(Env_shape1) - freq(Env_shape1, value=NA)
      change$species[i]<- spp[i]
      change$n_pixels[i]<-nc
      change$change_perc[i]<- (s*(100/nc))/100})
    
    n<-paste(orders[j],"Tree_canopy_change.csv", sep="_")
    write.csv(change, paste("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_change_tables", n, sep="/"))}}

Combining the tables

To facilitate the manipulation of this data, we need to combine the data:

rm(list=ls())

setwd("C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/Land_cover_proportions")

bare<-list.files(".", pattern = "Bare_ground_change")
all_bare <- Reduce(rbind, lapply(bare, read.csv))
all_bare<-all_bare[-(3:8)]
all_bare$data<-"Bare_ground"

tree<-list.files(".", pattern = "Tree_canopy_change")
all_tree <- Reduce(rbind, lapply(tree, read.csv))
# all_tree<-all_tree[-(3:8)]
all_tree$data<-"Tree_canopy"

short<-list.files(".", pattern = "Short_vegetation")
all_short <- Reduce(rbind, lapply(short, read.csv))
all_short<-all_short[-(3:8)]
all_short$data<-"Short_vegetation"

colnames(all_bare)
colnames(all_tree)
colnames(all_short)

tab<-rbind(all_bare, all_tree, all_short)
write.csv(tab, "C:/Users/raque/OneDrive - ufpr.br/Doutorado/3ºCap/prop_change_all_spp.csv")