This R Markdown document presents the functions and the steps used to estimate species’ Area of Habitat (AOH) as described by Dória & Dobrovolski (2020) in the paper “Improving post-2020 conservation of terrestrial vertebrates in Caatinga”.
install.packages ('sp')
install.packages ('raster')
install.packages ('maptools')
install.packages ('maps')
install.packages ('rworldmap')
install.packages ('colorRamps')
install.packages ('rgdal')
install.packages ('here')library (sp)
library (raster)
library (maptools)## Warning: package 'maptools' was built under R version 3.6.1## Checking rgeos availability: TRUElibrary (maps)
library (rworldmap)## ### Welcome to rworldmap ##### For a short introduction type : vignette('rworldmap')library (colorRamps)
library (rgdal)## rgdal: version: 1.3-4, (SVN revision 766)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/Thais/OneDrive/Documentos/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/Thais/OneDrive/Documentos/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-1library (here)## Warning: package 'here' was built under R version 3.6.1## here() starts at F:/DOC/PAPERS/Paper1/MATERIAL SUPLEMENTAR/RPubsf.hab.redlist <- function(z){
hab.spp <- rl_habitats(z, key="38e913256ac81b0159bc5421b633ad467b94be67e9c886a3ca86ccc78710ef50")
return(hab.spp$result$code)
}
#Notes:
#'z' in this function is matrix with specie's names
#in rl_habitats function, the argument 'key' must be requested from IUCN.f.pref.hab <- function(y, x){
for (i in 1:nrow(y))
y[i,"artif"] <- ifelse(sum(x[i,]>14 & x[i,]<15),1,0)
yart <- y
{
for (i in 1:nrow(yart))
yart[i,"opena"] <- ifelse(sum(x[i,]>2 & x[i,]<5 | x[i,]>=12.2 & x[i,]<=12.3 | x[i,]==13.3),1,0)
yopen <- yart
{
for (i in 1:nrow(yopen))
yopen[i,"fores"] <- ifelse(sum(x[i,]<=1.9),1,0)
return(yopen)}}
}
# Notes:
# 'y' in this function is a matrix with categories of habitat that will be considered as preference to each species
# 'x' in this function is a matrix with species code's indicating habitat preference according to IUCN databasef.habelev.ref <- function(a, b, c) {
nsp <- nlayers(a)
for (i in 1:nsp) {
tab2 <- cbind(categs.pref, (as.numeric(b[i,3:7])))
landcover.reclass2 <- reclassify(landcover.reclass, tab2, progress="text")
if(c[i,2]!=0){
elev.reclass <- elev.caat>=c[i,1]&elev.caat<=c[i,2]
landcover.reclass2 <- overlay(landcover.reclass2,elev.reclass,fun=function (x,y){x*y}, progress="text")
}
landelev.reclass.res <- aggregate(landcover.reclass2, fun=sum, fact=172, progress="text") #0.05 resolution
landelev.abfrac <- landelev.reclass.res/29584
spp <- a[[i]]
proj <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
proj4string(spp) <- proj
resamplesp <- resample(landelev.abfrac,spp,method="ngb", progress="text")
mult <- overlay(resamplesp,spp,fun=function (x,y){x*y})
nome <- unique[i] #species list created in STEP 1
writeRaster(mult,paste(nome,"0.05_HSMs.asc",sep = ""),format="ascii")
}}
#Notes:
# 'a' in this function is the original raster of species
# 'b' in this function is matrix of habitat's preference of species
# 'c' in this function is the matrix of elevation's preference of speciesadata <- "https://github.com/thaisdoria/TD_Thesis_Chapter1_RData/blob/master/DoriaTAF_Amphibians.RData?raw=true"
load(url(adata)) # RData file created on 2017/2018 containing the files listed bellow
unifiedshp.amphibians <- AMPHCAAT #amphibians species distribution (IUCN data)
shape.caat <- caatinga #shapefile of the region
landcover.caat <- mpbio.ori #raster of the region landcover
tab.reclass <- tab.reclass #table for reclassify landcover
pref.hab.spp <- habnull.a #table containing only preferences data on water bodies (manually assigned following specific criteria) that will be filled with habitat preferences from IUCN
elev.caat <- elev.resample #raster of the region elevation
pref.elev.spp <- elevspp.a #table with the interval of elevation in which each species occurs (data compiled manually from IUCN database)ATTENTION: polygons from IUCN were downloaded at 2017/2018 and after that the genus of species “Hypsiboas albomarginatus”, “Hypsiboas albopunctatus”, “Hypsiboas atlanticus”, was modified for “Boana”. Thus, to run this script adequately and obtain correct correspondences between the species names from shapefiles and the most recent habitat data from IUCN database it is recommended to download the up-dated polygons and not use those available in this RData.
##### STEP 1 - Creating an individual shapefile for each species distribution ####
data <- unifiedshp.amphibians
names(data) #used to verify the column that contained the species names
unique <- sort(unique(data@data$binomial)) #accessing 'binomial', column that contain the species names, and saving these names removing the duplicate ones
for (i in 1:length(unique)) {
tmp <- data[data$binomial== unique[i], ]
writeOGR(tmp, dsn=getwd(), unique[i], driver="ESRI Shapefile",
overwrite_layer=TRUE)
}
#### STEP 2 - Creating an individual raster for each species distribution ####
extent(landcover.caat)
r0.05 <- raster(res=0.05, xmn=-44.50838, xmx=-35.17174, ymn=-16.08856, ymx=-2.783879)
r0.05[] <- 1:ncell(r0.05)
r0.05 # 266, 187, 49742 (nrow, ncol, ncell)
lista.spp.a <- list.files(pattern = ".shp$", full.names = F) #list of amphibians shapefiles
# Presence/Absence
for (i in 1:length(lista.spp.a)) {
spp <- shapefile(lista.spp.a[i])
spp.r <- rasterize(spp, r0.05, background = 0, getCover=T)
spp.r[spp.r > 1] <- 1
spp.c <- crop(spp.r, shape.caat)
spp.m <- mask(spp.c, shape.caat)
e <- extent(-44.8, -34.8, -16.6, -2.2)
spp.e <- extend(spp.m,e)
nome <- unique[i]
writeRaster(spp.e, paste(nome, "0.05.asc", sep = ""), format = "ascii", overwrite = T)
}
#### STEP 3 - Obtaining information of habitat species preferences data from IUCN ####
spp.nomes <- as.matrix(unique) #remember: the 'unique' object was created in the STEP 1
f.hab.redlist <- function(spp.nomes){
hab.spp <- rl_habitats(spp.nomes, key="99ad94015eeccd35abe9f869c633e1ebad67ee390ab444902eead0c309c800d5")
return(hab.spp$result$code)
}
lista.hab <- lapply(1:length(spp.nomes),function(f0) f.hab.redlist(spp.nomes[f0])) #creating an object with all information about species habitat preference running 'hab_redlist' function which returns only the information about preference, while 'rl_habitats' function returns in addition the name of the species
spp.codes <- do.call(rbind, lista.hab) # transforming the list with habitat codes of species in a matrix of habitat codes
#### STEP 4 - Editing the habitat preference to obatin correspondence with the landcover categories ####
## Preparing the matrix
pref.hab.spp # table already loaded with file .RData
hab.pref <- as.matrix(pref.hab.spp)
# Remove lines for which there is no available IUCN data (verify this on objetc 'lista.hab')
hab.pref1 <- hab.pref[-109,] # the number here corresponds to lines corresponding to a species for which the habitat preference information was empty, that is, there wasn't habitat preference information for this species in IUCN database
hab.pref1 <- replace(hab.pref1, list = is.na(hab.pref1), values = 0) #replace NA for 0
## Using the function (to fill the matrix with IUCN data for 128 spp)
hab.pref2 <- f.pref.hab (hab.pref1, spp.codes) #creating an matrix with habitat preferences of each species with function f.pref.hab
## Editing the matrix of spp. without IUCN data (1 spp)
hab.pref3 <- hab.pref[109,]
hab.pref3 <- replace(hab.pref3, list = is.na(hab.pref3), values = 1) # replace NA data for 1 (to maintain original distribution on artificial, open areas and forest areas for species without IUCN data)
# Creating a final (total) matrix
pref.hab.spp.a <- rbind(hab.pref2, hab.pref3) # Join the two matrix (spp with and without IUCN data)
pref.hab.spp.a <- pref.hab.spp.a[order(as.numeric(as.factor(pref.hab.spp.a[,2]))),] # reorder in alphabetical order
## Saving the habitat matrix of species
write.table(pref.hab.spp.a, file="HabMat_Amphs_FILLED.csv", sep=";", header=T) #saving the matrix with habitat preferences of each species
#### STEP 5 - Reclassifying the land cover to obtain equivalence with the IUCN habitat categories ####
landcover.reclass <- reclassify(landcover.caat, tab.reclass)
writeRaster(landcover.reclass, filename="landcoverReclass.tif")
#### STEP 6 - Refining the individual raster according to the land cover and habitat preferences #####
## Preparing data
categs.pref <- c(0,1,2,3,4)
landcover.reclass <- raster("landcoverReclass.tif") # uploading the reclassified land cover
elev.caat # raster already loaded with file .RData
pref.hab.spp.a # table created in the STEP 4
pref.elev.spp # table already loaded with file .RData
pref.elev.spp.a <- as.matrix(pref.elev.spp)
pref.elev.spp.a <- matrix(as.numeric(pref.elev.spp.a[, 2:3]), ncol = 2 )
range.orig.a <- raster::stack(list.files(path=" ", pattern = "asc$", full.names = F)) #raster files created in the STEP 2
## Using the function
f.habelev.ref(range.orig.a, pref.hab.spp.a, pref.elev.spp.a)rdata <- "https://github.com/thaisdoria/TD_Thesis_Chapter1_RData/blob/master/DoriaTAF_Reptilesg.RData?raw=true"
load(url(rdata)) #RData file created on 2019 containing the files listed bellow:
unifiedshp.reptiles.gard <- REPCAATgard #reptiles species distribution (GARD data)
shape.caat <- caatinga #shapefile of the region
landcover.caat <- mpbio.ori #raster of the region landcover
tab.reclass <- tab.reclass #table for reclassify landcover function
pref.hab.spp <- habnull.rg #table containing only preferences data on water bodies (manually assigned following specific criteria) that will be filled with habitat preferences from IUCN
elev.caat <- elev.resample #raster of the region elevation
pref.elev.spp <- elevspp.rg #table with the interval of elevation in which each species occurs (data compiled manually from IUCN database)##### STEP 1 - Creating an individual shapefile for each species distribution ####
data <- unifiedshp.reptiles.gard
names(data) #used to verify the column that contained the species names
unique.gard <- sort(unique(data@data$Binomial)) #accessing 'binomial', column that contain the species names, and saving these names removing the duplicate ones
for (i in 1:length(unique.gard)) {
tmp <- data[data$Binomial== unique.gard[i], ]
writeOGR(tmp, dsn=getwd(), unique.gard[i], driver="ESRI Shapefile",
overwrite_layer=TRUE)
}
#### STEP 2 - Creating an individual raster for each species distribution ####
extent(landcover.caat) # this information was used to create a reference raster, the chosen resolution was 300m (0.002777778 x 0.002777778 degree)
r0.05 <- raster(res=0.05, xmn=-44.50838, xmx=-35.17174, ymn=-16.08856, ymx=-2.783879)
r0.05[] <- 1:ncell(r0.05)
r0.05 # 266, 187, 49742 (nrow, ncol, ncell)
lista.spp.r.gard <- list.files(pattern = ".shp$", full.names = F) #list of reptiles shapefiles
# Presence/Absence
for (i in 1:length(lista.spp.r.gard)) {
spp <- shapefile(lista.spp.r.gard[i])
spp.r <- rasterize(spp, r0.05, background = 0, getCover=T)
spp.r[spp.r > 1] <- 1
spp.c <- crop(spp.r, shape.caat)
spp.m <- mask(spp.c, shape.caat)
e <- extent(-44.8, -34.8, -16.6, -2.2)
spp.e <- extend(spp.m,e)
nome <- unique.gard[i]
writeRaster(spp.e, paste(nome, "0.05g.asc", sep = ""), format = "ascii", overwrite = T)
}
#### STEP 3 - Obtaining information of habitat species preferences data from IUCN ####
spp.nomes.gard <- as.matrix(unique.gard) #remember: the 'unique.gard' object was created in the STEP 1
f.hab.redlist <- function(spp.nomes.gard){
hab.spp <- rl_habitats(spp.nomes.gard, key="99ad94015eeccd35abe9f869c633e1ebad67ee390ab444902eead0c309c800d5")
return(hab.spp$result$code)
}
lista.hab.gard <- lapply(1:length(spp.nomes.gard),function(f0) f.hab.redlist(spp.nomes.gard[f0])) #creating an object with all information about species habitat preference running 'hab_redlist' function which returns only the information about preference, while 'rl_habitats' function returns in addition the name of the species
spp.codes.gard <- do.call("rbind", lista.hab.gard) # transforming the list with habitat codes of species in a matrix of habitat codes
spp.codes.gard
#### STEP 4 - Editing the habitat preference to obatin correspondence with the landcover categories ####
## Preparing the matrix
pref.hab.spp # c table already loaded with file .RData
hab.pref.gard <- as.matrix(pref.hab.spp)
# There are IUCN data only for 82 spp.
# Remove lines for which there is no available IUCN data (verify this on objetc 'lista.hab.gard')
hab.pref.gard1 <- hab.pref.gard[-c(1,3:6,8:10,14:16, 18, 20:30, 32:33, 36:49, 51:61, 63:65, 68, 70, 72:74, 78, 80, 82:85, 87:100, 103:104, 107, 109, 113:115, 117:122, 127:140, 142, 145, 148, 150:152, 155:160, 162:174, 177:188, 190:192, 194:196, 198:208, 211, 213:215, 217:219, 222, 226:229, 232:236, 240:241, 247:251, 253, 255:258, 260:261, 264, 266:268, 270, 276:281, 283, 285),] # the number here corresponds to lines corresponding to a species for which the habitat preference information was empty, that is, there wasn't habitat preference information for this species in IUCN database
hab.pref.gard1 <- replace(hab.pref.gard1, list = is.na(hab.pref.gard1), values = 0) #replace NA for 0
## Using the function (to fill the matrix with IUCN data for 82 spp)
hab.pref.gard2 <- f.pref.hab (hab.pref.gard1, spp.codes.gard) #creating an matrix with habitat preferences of each species with function f.pref.hab
## Identifing the spp. which IUCN data have no equivalence with LULC class (i.e. with data that not match with artificial, open and forest areas from MAPBIOMAS)
hab.pref.gard3 <- hab.pref.gard2[hab.pref.gard2[,4] == 0 & hab.pref.gard2[,5] == 0 & hab.pref.gard2[,6] == 0, ]
hab.pref.gard3[,4:6] <- 1 # replace NA data for 1 (to maintain original distribution on artificial, open areas and forest areas for species without IUCN data)
filtro1=match(as.character(hab.pref.gard3[,2]),as.character(hab.pref.gard2[,2]))
hab.pref.gard2=hab.pref.gard2[-c(filtro1),]
hab.pref.gard2
## Editing the matrix of spp. without IUCN data (203 spp)
hab.pref.gard4 <- hab.pref.gard[c(1,3:6,8:10,14:16, 18, 20:30, 32:33, 36:49, 51:61, 63:65, 68, 70, 72:74, 78, 80, 82:85, 87:100, 103:104, 107, 109, 113:115, 117:122, 127:140, 142, 145, 148, 150:152, 155:160, 162:174, 177:188, 190:192, 194:196, 198:208, 211, 213:215, 217:219, 222, 226:229, 232:236, 240:241, 247:251, 253, 255:258, 260:261, 264, 266:268, 270, 276:281, 283, 285),]
hab.pref.gard4 <- replace(hab.pref.gard4, list = is.na(hab.pref.gard4), values = 1) # replace NA data for 1 (to maintain original distribution on artificial, open areas and forest areas for species without IUCN data)
# Creating a final (total) matrix
pref.hab.spp.rgard <- rbind(hab.pref.gard2, hab.pref.gard3, hab.pref.gard4) # Join the two matrix (spp with and without IUCN data)
pref.hab.spp.rgard <- pref.hab.spp.rgard[order(as.numeric(as.factor(pref.hab.spp.rgard[,2]))),] # reorder in alphabetical order
## Saving the habitat matrix of species
write.table(pref.hab.spp.rgard, file="HabMat_RepsGard_FILLED.csv", sep=";") #saving the matrix with habitat preferences of each species
#### STEP 5 - Reclassifying the land cover to obtain equivalence with the IUCN habitat categories ####
landcover.reclass <- reclassify(landcover.caat, tab.reclass)
writeRaster(landcover.reclass, filename="landcoverReclass.tif")
plot(landcover.reclass)
#### STEP 6 - Refining the individual raster according to the land cover and habitat preferences #####
## Preparing data
categs.pref <- c(0,1,2,3,4)
landcover.reclass <- raster("landcoverReclass.tif") # uploading the reclassified land cover
elev.caat # raster already loaded with file .RData
elev.caat <- raster("TopoCaat_resample.asc")
pref.hab.spp.rgard # table created in the STEP 4
pref.elev.spp # renamed object from table already loaded with file .RData
pref.elev.spp.rgard <- as.matrix(pref.elev.spp)
pref.elev.spp.rgard <- matrix(as.numeric(pref.elev.spp.rgard[,3:4]), ncol = 2 )
range.orig.rgard <- raster::stack(list.files(path=" ", pattern = "asc$", full.names = F)) #raster files created in the STEP 2
## Using the function
f.habelev.refgard(range.orig.rgard, pref.hab.spp.rgard, pref.elev.spp.rgard)avdata <- "https://github.com/thaisdoria/TD_Thesis_Chapter1_RData/blob/master/DoriaTAF_Aves.RData?raw=true"
load(url(avdata)) # RData file created on 2017/2018 containing the files listed bellow:
#### Running steps ####
unifiedshp.aves <- AVESCAAT #aves species distribution (IUCN data)
shape.caat <- caatinga #shapefile of the region
landcover.caat <- mpbio.ori #raster of the region landcover
tab.reclass <- tab.reclass #table for reclassify landcover function
pref.hab.spp <- habnull.av #table containing only preferences data on water bodies (manually assigned following specific criteria) that will be filled with habitat preferences from IUCN
elev.caat <- elev.resample #raster of the region elevation
pref.elev.spp <- elevspp.av #table with the interval of elevation in which each species occurs (data compiled manually from IUCN database)ATTENTION: polygons from IUCN were downloaded at 2017/2018 and after that the genus of specie “Phalacrocorax brasilianus” was modified for “Nannopterum”. Thus, to run this script adequately and obtain correct correspondences between the species names from shapefiles and the most recent habitat data from IUCN database it is recommended to download the up-dated polygons and not use those available in this RData.
#### STEP 1 - Creating an individual shapefile for each species distribution ####
data <- unifiedshp.aves # repeat this step to all groups (amphibians, reptiles, mammals, aves)
names(data) #used to verify the column that contained the species names
unique <- sort(unique(data@data$SCINAME)) #accessing 'SCINAME', column that contain the species names, and saving these names removing the duplicate ones
for (i in 1:length(unique)) {
tmp <- data[data$SCINAME== unique[i], ]
writeOGR(tmp, dsn=getwd(), unique[i], driver="ESRI Shapefile",
overwrite_layer=TRUE)
}
#### STEP 2 - Creating an individual raster for each species distribution ####
extent(landcover.caat) # this information was used to create a reference raster, the chosen resolution was 300m (0.002777778 x 0.002777778 degree)
r0.05 <- raster(res=0.05, xmn=-44.50838, xmx=-35.17174, ymn=-16.08856, ymx=-2.783879)
r0.05[] <- 1:ncell(r0.05)
r0.05 # 266, 187, 49742 (nrow, ncol, ncell)
lista.spp.av <- list.files(pattern = ".shp$", full.names = F) # list of aves shapefiles
# Presence/Absence
for (i in 1:length(lista.spp.av)) {
spp <- shapefile(lista.spp.av[i])
spp.r <- rasterize(spp, r0.05, background = 0, getCover=T)
spp.r[spp.r > 1] <- 1
spp.c <- crop(spp.r, shape.caat)
spp.m <- mask(spp.c, shape.caat)
e <- extent(-44.8, -34.8, -16.6, -2.2)
spp.e <- extend(spp.m,e)
nome <- unique[i]
writeRaster(spp.e, paste(nome, "0.05.asc", sep = ""), format = "ascii", overwrite = T)
}
#### STEP 3 - Obtaining information of habitat species preferences data from IUCN ####
spp.nomes <- as.matrix(unique) #remember: the 'unique' object was created in the STEP 1
f.hab.redlist <- function(spp.nomes){
hab.spp <- rl_habitats(spp.nomes, key="99ad94015eeccd35abe9f869c633e1ebad67ee390ab444902eead0c309c800d5")
return(hab.spp$result$code)
}
lista.hab <- lapply(1:length(spp.nomes),function(f0) f.hab.redlist(spp.nomes[f0])) #creating an object with all information about species habitat preference running 'hab_redlist' function which returns only the information about preference, while 'rl_habitats' function returns in addition the name of the species
spp.codes <- do.call(rbind, lista.hab) # transforming the list with habitat codes of species in a matrix of habitat codes
#### STEP 4 - Editing the habitat preference to obatin correspondence with the landcover categories ####
## Preparing the matrix
pref.hab.spp # table already loaded with file .RData
hab.pref <- as.matrix(pref.hab.spp)
# There are IUCN data for all aves species
hab.pref1 <- replace(hab.pref, list = is.na(hab.pref), values = 0) #replace NA for 0
## Using the function
hab.pref2 <- f.pref.hab (hab.pref1, spp.codes) #creating an matrix with habitat preferences of each species with function f.pref.hab
## Identifing the spp. which IUCN data have no equivalence with LULC class (i.e. with data that not match with artificial, open and forest areas from MAPBIOMAS)
hab.pref3 <- hab.pref2[hab.pref2[,4] == 0 & hab.pref2[,5] == 0 & hab.pref2[,6] == 0, ]
hab.pref3[,4:6] <- 1 # replace NA data for 1 (to maintain original distribution on artificial, open areas and forest areas for species without IUCN data)
filtro1=match(as.character(hab.pref3[,2]),as.character(hab.pref2[,2]))
hab.pref2=hab.pref2[-c(filtro1),]
# Creating a final (total) matrix
pref.hab.spp.av <- rbind(hab.pref2, hab.pref3) # Join the two matrix (spp with and without IUCN data)
pref.hab.spp.av <- pref.hab.spp.av[order(as.numeric(as.factor(pref.hab.spp.av[,2]))),] # reorder in alphabetical order
## Saving the habitat matrix of species
write.table(pref.hab.spp.m, file="HabMat_Aves_FILLED.csv", sep=";") #saving the matrix with habitat preferences of each species
#### STEP 5 - Reclassifying the land cover to obtain equivalence with the IUCN habitat categories ####
landcover.reclass <- reclassify(landcover.caat, tab.reclass)
writeRaster(landcover.reclass, filename="landcoverReclass.tif")
#### STEP 6 - Refining the individual raster according to the land cover and habitat preferences #####
## Preparing data
categs.pref <- c(0,1,2,3,4)
landcover.reclass <- raster("landcoverReclass.tif") # uploading the reclassified land cover
elev.caat #raster already loaded with file .RData
pref.hab.spp.av # table created in the STEP 4
pref.elev.spp # renamed object from table already loaded with file .RData
pref.elev.spp.av <- as.matrix(pref.elev.spp)
pref.elev.spp.av <- matrix(as.numeric(pref.elev.spp.av[, 2:3]), ncol = 2 )
range.orig.av <- raster::stack(list.files(path=" ", pattern = "asc$", full.names = F))
## Using the function
f.habelev.ref(range.orig.av, pref.hab.spp.av, pref.elev.spp.av)mdata <- "https://github.com/thaisdoria/TD_Thesis_Chapter1_RData/blob/master/DoriaTAF_Mammals.RData?raw=true"
load(url(mdata)) # RData file created on 2017/2018 containing the files listed bellow:
unifiedshp.mammals <- MAMICAAT #mammals species distribution (IUCN data)
shape.caat <- caatinga #shapefile of the region
landcover.caat <- mpbio.ori #raster of the region landcover
tab.reclass <- tab.reclass #table for reclassify landcover function
pref.hab.spp <- habnull.m #table containing only preferences data on water bodies (manually assigned following specific criteria) that will be filled with habitat preferences from IUCN
elev.caat <- elev.resample #raster of the region elevation
pref.elev.spp <- elevspp.m #table with the interval of elevation in which each species occurs (data compiled manually from IUCN database)ATTENTION: after the polygons from IUCN were downloaded at 2017/2018 and after that the genus of specie “Mimon crenulatum” was modified for “Gardnerycteris”. Thus, to run this script adequately and obtain correct correspondences between the species names from shapefiles and the most recent habitat data from IUCN database it is recommended to download the up-dated polygons and not use those available in this RData.
#### STEP 1 - Creating an individual shapefile for each species distribution ####
data <- unifiedshp.mammals
names(data) #used to verify the column that contained the species names
unique <- sort(unique(data@data$binomial)) #accessing 'binomial', column that contain the species names, and saving these names removing the duplicate ones
for (i in 1:length(unique)) {
tmp <- data[data$binomial== unique[i], ]
writeOGR(tmp, dsn=getwd(), unique[i], driver="ESRI Shapefile",
overwrite_layer=TRUE)
}
#### STEP 2 - Creating an individual raster for each species distribution ####
extent(landcover.caat) # this information was used to create a reference raster, the chosen resolution was 300m (0.002777778 x 0.002777778 degree)
r0.05 <- raster(res=0.05, xmn=-44.50838, xmx=-35.17174, ymn=-16.08856, ymx=-2.783879)
r0.05[] <- 1:ncell(r0.05)
r0.05 # 266, 187, 49742 (nrow, ncol, ncell)
lista.spp.m <- list.files(pattern = ".shp$", full.names = F) #list of mammals shapefiles
# Presence/Absence
for (i in 1:length(lista.spp.m)) {
spp <- shapefile(lista.spp.m[i])
spp.r <- rasterize(spp, r0.05, background = 0, getCover=T)
spp.r[spp.r > 1] <- 1
spp.c <- crop(spp.r, shape.caat)
spp.m <- mask(spp.c, shape.caat)
e <- extent(-44.8, -34.8, -16.6, -2.2)
spp.e <- extend(spp.m,e)
nome <- unique[i]
writeRaster(spp.e, paste(nome, "0.05.asc", sep = ""), format = "ascii", overwrite = T)
}
#### STEP 3 - Obtaining information of habitat species preferences data from IUCN ####
spp.nomes <- as.matrix(unique) #remember: the 'unique' object was created in the STEP 1
f.hab.redlist <- function(spp.nomes){
hab.spp <- rl_habitats(spp.nomes, key="99ad94015eeccd35abe9f869c633e1ebad67ee390ab444902eead0c309c800d5")
return(hab.spp$result$code)
}
lista.hab <- lapply(1:length(spp.nomes),function(f0) f.hab.redlist(spp.nomes[f0])) #creating an object with all information about species habitat preference running 'hab_redlist' function which returns only the information about preference, while 'rl_habitats' function returns in addition the name of the species
spp.codes <- do.call(rbind, lista.hab) # transforming the list with habitat codes of species in a matrix of habitat codes
#### STEP 4 - Editing the habitat preference to obatin correspondence with the landcover categories ####
## Preparing the matrix
pref.hab.spp # table already loaded with file .RData
hab.pref <- as.matrix(pref.hab.spp)
# There are IUCN data for all 208 mammals
hab.pref1 <- replace(hab.pref, list = is.na(hab.pref), values = 0) #replace NA for 0
## Using the function
hab.pref2 <- f.pref.hab (hab.pref1, spp.codes) #creating an matrix with habitat preferences of each species with function f.pref.hab
## Identifing the spp. which IUCN data have no equivalence with LULC class (i.e. with data that not match with artificial, open and forest areas from MAPBIOMAS)
hab.pref3 <- hab.pref2[hab.pref2[,4] == 0 & hab.pref2[,5] == 0 & hab.pref2[,6] == 0, ]
hab.pref3[,4:6] <- 1 # replace NA data for 1 (to maintain original distribution on artificial, open areas and forest areas for species without IUCN data)
filtro1=match(as.character(hab.pref3[,2]),as.character(hab.pref2[,2]))
hab.pref2=hab.pref2[-c(filtro1),]
# Creating a final (total) matrix
pref.hab.spp.m <- rbind(hab.pref2, hab.pref3) # Join the two matrix (spp with and without IUCN data)
pref.hab.spp.m <- pref.hab.spp.m[order(as.numeric(as.factor(pref.hab.spp.m[,2]))),] # reorder in alphabetical order
## Saving the habitat matrix of species
write.table(pref.hab.spp.m, file="HabMat_Mammals_FILLED.csv", sep=";") #saving the matrix with habitat preferences of each species
#### STEP 5 - Reclassifying the land cover to obtain equivalence with the IUCN habitat categories ####
landcover.reclass <- reclassify(landcover.caat, tab.reclass)
writeRaster(landcover.reclass, filename="landcoverReclass.tif")
#### STEP 6 - Refining the individual raster according to the land cover and habitat preferences #####
## Preparing data
categs.pref <- c(0,1,2,3,4)
landcover.reclass <- raster("landcoverReclass.tif") # uploading the reclassified land cover
elev.caat # raster already loaded with file .RData
pref.hab.spp.m # table created in the STEP 4
pref.elev.spp # renamed object from table already loaded with file .RData
pref.elev.spp.m <- as.matrix(pref.elev.spp)
pref.elev.spp.m <- matrix(as.numeric(pref.elev.spp.m[, 2:3]), ncol = 2 )
range.orig.m <- raster::stack(list.files(path=" ", pattern = "asc$", full.names = F)) #raster files created in the STEP 2
# Using the function
f.habelev.ref(range.orig.m, pref.hab.spp.m, pref.elev.spp.m)