This R Markdown document contains all of the code accompanying the
publication Over-the-horizon threat evaluation predicts rapidly
shifting geographic and taxonomic priorities for the conservation of
Australian tetrapods. Annotation is included to explain the
pipelines in detail. Please note these analyses are computationally
extensive and extremely time-consuming.
In order to run the code properly, make sure you have the latest Java
Development Environment (JDE; https://www.oracle.com/java/technologies/downloads/)
installed on your computer, and that your working directory contains
maxent.jar (https://biodiversityinformatics.amnh.org/open_source/maxent/)
and the following folders, which can be downloaded from the associated
data repository:
occurrences - contains all occurrence records downloaded from ALA as individual CSV files per species.
invasive - contains the above subfolder but only for the invasive species.
RAW - contains all raw data used in the subsequent data preparation steps. These are not available in the data repository, and can be sourced directly using the citations provided in the text.
SDM - contains all of the environmental layers used in the SDM.
AA - contains the output from the Automated Assessment procedure, including results of analyses done on the AA output.
SDM - contains the output from the Species Distribution Modelling. Subdivided into one folder with plots of individual response curves per species (Response Curves), one folder for SDMs on trained and projected to current climatic conditions (Current), and one folder for SDMs projected to future climatic conditions (Future). Each of the latter two contains three sub-folders:
Maps - contains pdf files with individual maps for each species, with the suffix _pa denoting presence-absence based on an optimal threshold, and the suffix _prob denoting predicted probability. In the Current folder _pa maps also include black dots denoting current occurrence record.
Rasters - contains individual rasters for each species in .geotiff format, with the suffix _pa denoting presence-absence based on an optimal threshold, and the suffix _pr denoting predicted probability. In the Future folder each raster is a stack of projections under different climate change scenarios.
Shapefiles - contains individual shapefiles of the presence-absence based on an optimal threshold for each species. In the Future folder each shapefile is a multipolygon of projections under different climate change scenarios.
The first step is data preparation.
The raw distribution data are sourced directly from the IUCN (mammals
and amphibians), BirdLife (birds), and the GARD Initiative
(reptiles).
The code in this chunk filters the raw shapefiles to only include native
Australian species, downloads missing IUCN Red List Assessments,
downloads occurrence data from ALA, and downloads and prepares
environmental predictors. Please note that to run this code you will
require an API token for the IUCN Red List (https://apiv3.iucnredlist.org/api/v3/token), and an
email associated with a registered account in Atlas of Living Australia
(https://auth.ala.org.au/userdetails/registration/createAccount).
library(tidyverse) # install.packages("tidyverse")
library(sf) # install.packages("sf")
library(rnaturalearth) # install.packages("rnaturalearth")
library(terra) # install.packages("terra")
library(galah) # install.packages("galah")
library(taxize) # install.packages("taxize")
library(foreach) # install.packages("foreach")
library(doParallel) # install.packages("doParallel")
sf_use_s2(FALSE)
## Create shapefiles of Australian tetrapods
# import shapefiles of global verteberate distributions
mam_shp <-
bind_rows(
read_sf("data/RAW/MAMMALS_TERRESTRIAL_ONLY.shp"),
read_sf("data/RAW/MAMMALS_FRESHWATER.shp")
)
bird_shp <- readRDS("data/RAW/BirdMaps_JetzNames_filtered.rds") %>%
filter(!Jetz_name %in% as.character(read_csv("data/birds_to_drop.csv") %>% pull(Binomial))) %>% # drop seabirds
mutate(
Jetz_name = case_when(
Jetz_name == "Philemon novaeguineae" ~ "Philemon buceroides",
Jetz_name == "Cacatua roseicapilla" ~ "Eolophus roseicapilla",
Jetz_name == "Casmerodius albus" ~ "Ardea alba",
Jetz_name == "Eopsaltria pulverulenta" ~ "Peneoenanthe pulverulenta",
Jetz_name == "Hirundo nigricans" ~ "Petrochelidon nigricans",
Jetz_name == "Ixobrychus minutus" ~ "Ixobrychus dubius",
Jetz_name == "Limicola falcinellus" ~ "Calidris flacinellus",
Jetz_name == "Tadorna radjah" ~ "Radjah radjah",
Jetz_name == "Sterna nereis" ~ "Sternula nereis",
Jetz_name == "Monarcha trivirgatus" ~ "Symposiachrus trivirgatus",
Jetz_name == "Megalurus timoriensis" ~ "Cincloramphus timoriensis",
Jetz_name == "Manucodia keraudrenii" ~ "Phonygammus keraudrenii",
Jetz_name == "Hirundo daurica" ~ "Cecropis daurica",
Jetz_name == "Alcedo pusilla" ~ "Ceyx pusillus",
TRUE ~ Jetz_name
)
) # make some fixes to taxonomy
rep_shp <- read_sf("data/RAW/Gard_1_7_ranges.shp") %>%
mutate(
Jetz_name = case_when(
binomial == "Ctenotus aphrodite" ~ "Ctenotus septenarius",
binomial == "Lerista talpina" ~ "Lerista petersoni",
binomial == "Pogona mitchelli" ~ "Pogona minor",
binomial == "Varanus similis" ~ "Varanus scalaris",
TRUE ~ binomial
)
) # make some fixes to taxonomy
amph_shp <- read_sf("data/RAW/ANURA.shp")
# create map of Australia
oz_map <-
ne_countries(country = "Australia",
scale = 10,
returnclass = "sf") %>%
st_transform(st_crs(mam_shp)) %>%
st_simplify(preserveTopology = FALSE, dTolerance = 0.01)
# generate a vector of polygons that intersect Australia
sel_oz = st_intersects(x = mam_shp, y = oz_map)
sel_logical = lengths(sel_oz) > 0
# filter out polygons that don't intersect Australia, non-native species, then combine all polygons to create a single polygon per species
mam_dist <- mam_shp[sel_logical,] %>%
st_intersection(oz_map) %>%
filter(!Binomial == "Bos javanicus") %>% # introduced in Australia
filter(presence == 1,
origin %in% c(1, 2)) %>%
group_by(binomial, category) %>%
summarise() %>%
mutate(binomial = str_replace(binomial, " ", "_")) %>%
rename(Binomial = binomial) %>%
# fix to taxonomy and distribution of Dobsonia moluccensis
filter(!Binomial == "Dobsonia_magna") %>%
mutate(Binomial = case_when(Binomial == "Dobsonia_moluccensis" ~ "Dobsonia_magna",
T ~ Binomial))
# generate a vector of polygons that intersect Australia
sel_oz = st_intersects(x = bird_shp, y = oz_map)
sel_logical = lengths(sel_oz) > 0
# filter out polygons that don't intersect Australia, non-native species, then combine all polygons to create a single polygon per species
bird_dist <- bird_shp[sel_logical,] %>%
st_intersection(oz_map) %>%
filter(PRESENC == 1,
ORIGIN %in% c(1, 2)) %>%
group_by(Jetz_name) %>%
summarise() %>%
mutate(Jetz_name = str_replace(Jetz_name, " ", "_")) %>%
rename(Binomial = Jetz_name) %>%
filter(!Binomial %in% c("Dacelo_tyro", "Elanus_caeruleus", "Hieraaetus_weiskei", "Malurus_cyanocephalus", "Oceanites_maorianus", "Oedistoma_pygmaeum", "Sipodotus_wallacii")) # drop some species that do not actually occur in Australia
# generate a vector of polygons that intersect Australia
sel_oz = st_intersects(x = rep_shp, y = oz_map)
sel_logical = lengths(sel_oz) > 0
# filter out polygons that don't intersect Australia, non-native species, then combine all polygons to create a single polygon per species
rep_dist <- rep_shp[sel_logical,] %>%
st_intersection(oz_map) %>%
group_by(binomial) %>%
summarise() %>%
mutate(binomial = str_replace(binomial, " ", "_")) %>%
rename(Binomial = binomial)
# generate a vector of polygons that intersect Australia
sel_oz = st_intersects(x = amph_shp, y = oz_map)
sel_logical = lengths(sel_oz) > 0
# filter out polygons that don't intersect Australia, non-native species, then combine all polygons to create a single polygon per species
amph_dist <- amph_shp[sel_logical,] %>%
st_intersection(oz_map) %>%
filter(presence == 1,
origin %in% c(1, 2)) %>%
group_by(binomial, category) %>%
summarise() %>%
mutate(binomial = str_replace(binomial, " ", "_")) %>%
rename(Binomial = binomial)
## Download IUCN Red List Assessments
# set API for IUCN RedList
API = "YOUR_TOKEN_HERE" ## replace with your IUCN Red List API token
# query IUCN categories (for birds and reptiles only, already included in IUCN shapefiles for mammals and amphibians) and add to shapefiles
IUCN.bird <-
iucn_summary(str_replace(bird_dist$Binomial, "_", " "), key = API)
IUCN.bird <- iucn_status(IUCN.bird) %>%
as.data.frame() %>%
rownames_to_column(var = "Binomial") %>%
as_tibble() %>%
mutate(Binomial = str_replace(Binomial, " ", "_")) %>%
set_names(c("Binomial", "category"))
bird_dist <- bird_dist %>%
left_join(IUCN.bird) %>%
arrange(Binomial)
IUCN.rep <-
iucn_summary(str_replace(rep_dist$Binomial, "_", " "), key = API)
IUCN.rep <- iucn_status(IUCN.rep) %>%
as.data.frame() %>%
rownames_to_column(var = "Binomial") %>%
as_tibble() %>%
mutate(Binomial = str_replace(Binomial, " ", "_")) %>%
set_names(c("Binomial", "category"))
# manually update some turtle assessments
rep_dist <- rep_dist %>%
left_join(IUCN.rep) %>%
mutate(
category = case_when(
Binomial == "Chelodina_longicollis" ~ "LC",
Binomial == "Chelodina_steindachneri" ~ "LC",
Binomial == "Elseya_dentata" ~ "LC",
Binomial == "Emydura_macquarii" ~ "LC",
Binomial == "Emydura_subglobosa" ~ "LC",
Binomial == "Emydura_victoriae" ~ "LC",
Binomial == "Myuchelys_latisternum" ~ "LC",
category == "LR/lc" ~ "LC",
category == "LR/nt" ~ "NT",
TRUE ~ category
)
) %>% # manually update evaluated but not listed (TTWG)
arrange(Binomial)
# export Australian tetrapod shapefiles
write_rds(mam_dist, "data/mam_dist.rds")
write_rds(bird_dist, "data/bird_dist.rds")
write_rds(rep_dist, "data/rep_dist.rds")
write_rds(amph_dist, "data/amph_dist.rds")
## generate layers of environmental predictors for SDMs
# generate vector of current WorldClim2 layers
wc_2 <- dir("data/RAW/WC2/25m", full.names = T)
# import and stack rasters
currentEnv <- rast(wc_2)
# clip rasters to Australia
currentEnv <- crop(currentEnv, vect(oz_map))
currentEnv <- mask(currentEnv, vect(oz_map))
# select variables to include in SDM as predictors
currentEnv <- currentEnv[[c(1, 14, 2, 3, 4, 7, 8, 9, 20)]]
# import soil layers
soils <- rast(dir("data/RAW/soils/", full.names = T))
# resample soils to be in the same extent and resolution as other layers (aggregating into lower resolution by avergaing cell values)
soils <-
resample(soils, currentEnv[[1]], method = "bilinear")
# add Australian soil layers
currentEnv <- c(currentEnv, soils)
names(currentEnv) <-
c(
"bio1",
"bio4",
"bio10",
"bio11",
"bio12",
"bio15",
"bio16",
"bio17",
"elev",
"AWC",
"BDW",
"CLY",
"pHc",
"SLT",
"SND",
"SOC"
)
# export raster stack for SDM
writeRaster(
currentEnv,
filename = "data/SDM/Current",
format = "GTiff",
overwrite = T
)
# generate vector of future WorldClim2 layers
future <- dir("data/RAW/WC2/HadGEM3-GC31-LL/", full.names = T)
# generate raster stacks for SDMs
for (i in 1:length(future)) {
# import and stack future layers (along with elevation layer from current) and select variables
futureEnv <-
c(rast(future[[i]]), rast(wc_2[[20]]))[[c(1, 4, 10, 11, 12, 15, 16, 17, 20)]]
# clip rasters to Australia
futureEnv <- crop(futureEnv, vect(oz_map))
futureEnv <- mask(futureEnv, vect(oz_map))
# add soil layer
futureEnv <- c(futureEnv, soils)
names(futureEnv) <-
c(
"bio1",
"bio4",
"bio10",
"bio11",
"bio12",
"bio15",
"bio16",
"bio17",
"elev",
"AWC",
"BDW",
"CLY",
"pHc",
"SLT",
"SND",
"SOC"
)
# export raster stack for SDM
writeRaster(
futureEnv,
filename = paste(
"data/SDM/",
str_match(future[[i]], "bioc_\\s*(.*?)\\s*.tif")[2],
sep = ""
),
format = "GTiff",
overwrite = T
)
}
## create dataset of occurrence points in Australia
# generate vector of taxa
oz_tax <-
c(mam_dist$Binomial,
bird_dist$Binomial,
rep_dist$Binomial,
amph_dist$Binomial)
# download occurrences per species
cores = detectCores()
cl <-
makeCluster(cores[1] - 2) ## set to use all but one thread - replace if necessary
registerDoParallel(cl)
tibble(Binomial = character(),
n_pre = numeric(),
n_post = numeric(),
date = character(),
doi = character()) %>%
write_csv("data/sample_sizes.csv")
foreach(i = 1:length(oz_tax),
.packages = c("tidyverse", "galah", "sf", "dismo", "terra")) %dopar% {
sf_use_s2(FALSE)
# test if already downloaded (if so, skip to next iteration)
if (!oz_tax[[i]] %in% str_remove(dir("data/occurrences/"), ".csv")) {
# set email registered to ALA
galah_config(email = "YOUR_EMAIL_HERE") ## replace with your registered email in ALA
result <- tryCatch(
galah_call() |>
galah_identify(str_replace(oz_tax[[i]], "_", " ")) |>
galah_filter(profile = "ALA", year > 1990) |> # filter after 1990
galah_select(basisOfRecord, group = "basic") |>
atlas_occurrences(mint_doi = T),
error = function(e)
e
)
if (!inherits(result, "error")) {
# if there are occurrence points left after filtering
# convert to spatial object
result <-
st_as_sf(
result %>% drop_na(decimalLatitude, decimalLongitude, scientificName),
coords = c("decimalLongitude", "decimalLatitude"),
crs = st_crs(oz_map)
)
# filter shapefile
result_shp <-
bind_rows(mam_dist, bird_dist, rep_dist, amph_dist) %>% filter(Binomial == oz_tax[[i]]) %>%
st_transform(crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=km +no_defs") %>%
st_buffer(dist = 100) %>% # numeric buffer of 100 km to capture points just outside of ranges
st_transform(crs = st_crs(oz_map)) %>%
st_make_valid() %>% # fix geometry errors due to buffer
st_intersection(oz_map)
# filter out points which do not fall within IUCN shapefile
result_points <-
st_join(result, result_shp, join = st_within) %>%
filter(!is.na(Binomial))
if (dim(result_points)[1] > 0) {
# if there are occurrence points left after filtering
# import current environmental layers
currentEnv <- rast("data/SDM/Current.tif")
names(currentEnv) <-
c(
"bio1",
"bio4",
"bio10",
"bio11",
"bio12",
"bio15",
"bio16",
"bio17",
"elev",
"AWC",
"BDW",
"CLY",
"pHc",
"SLT",
"SND",
"SOC"
)
# perform spatial thinning on points
set.seed(42)
if(dim(result_points)[1] > 1){
points_occ <-
dismo::gridSample(as.matrix(st_coordinates(result_points)), currentEnv, n = 1)
} else {
points_occ <- st_coordinates(result_points)
}
# convert to tibble
result_points <-
as_tibble(points_occ) %>%
mutate(Binomial = oz_tax[[i]]) %>%
rename(decimalLongitude = X,
decimalLatitude = Y)
citation_text <- atlas_citation(result)
doi <- str_extract(citation_text, "https://doi\\.org/[0-9a-zA-Z./\\-]+")
# save sample sizes
tibble(Binomial = oz_tax[[i]],
n_pre = nrow(result),
n_post = nrow(result_points),
date = as.character(today()),
doi = doi) %>%
write_csv("data/sample_sizes.csv", append = T)
# export points
write_csv(result_points,
paste("data/occurrences/", oz_tax[[i]], ".csv", sep = ""))
}
}
}
}
stopCluster(cl)
# convert to single dataframe
oz_points <-
bind_rows(lapply(dir("data/occurrences/"), function(x)
read_csv(paste("data/occurrences/", x, sep = "")) %>% mutate(Binomial = str_remove(x, ".csv")))) %>%
mutate(
class = case_when(
Binomial %in% mam_dist$Binomial ~ "Mammal",
Binomial %in% bird_dist$Binomial ~ "Bird",
Binomial %in% rep_dist$Binomial ~ "Reptile",
Binomial %in% amph_dist$Binomial ~ "Amphibian"
)
) %>%
filter(decimalLongitude < 154) # drop Lord Howe Island
# export points
write_rds(oz_points %>% filter(class == "Mammal"), "data/mam_points.RDS")
write_rds(oz_points %>% filter(class == "Bird"), "data/bird_points.RDS")
write_rds(oz_points %>% filter(class == "Reptile"), "data/rep_points.RDS")
write_rds(oz_points %>% filter(class == "Amphibian"), "data/amph_points.RDS")Next we run species distribution models (SDMs). We use an ensemble
approach implemented through biomod2. We have run into some
issues running these on RStudio with R v4.2 - if this happens, the
scripts can be run directly through R without using RStudio.
These scripts will only run on recent versions of biomod2
(>= 4.2-5) - these scripts will not work with older versions because
of changes made to how modelling options and tuning are defined!
# install.packages("ENMeval")
# install.packages("dismo")
# install.packages("biomod2")
# install.packages("blockCV")
# install.packages("parallel")
# install.packages("pbapply")
# install.packages("doSNOW")
source("scripts/SDM_amph.R")
source("scripts/SDM_bird.R")
source("scripts/SDM_mam.R")
source("scripts/SDM_rep.R")
bind_rows(
read_csv("output/SDM/amph_summary.csv"),
read_csv("output/SDM/bird_summary.csv"),
read_csv("output/SDM/mam_summary.csv"),
read_csv("output/SDM/rep_summary.csv")
) %>%
write_csv("output/SDM/tet_summary.csv")Now that we have generated SDMs for all species and predicted future
ranges, we can calculate the predictor variables which will be used for
the automated assessment procedure.
The next chunk of code will first source a script that reruns everything
we have done so far but for invasive species in Australia - this is
important as overlap with invasive species ranges will be one of our
predictors. As before, a registered email address needs to be provided
to source data from ALA.
The chunk then runs code which sources the raw human population density
and land use data for different Shared Socioeconomic Pathways (SSPs) and
Global Circulation Models (GCMs) and overlaps these with current and
predicted future ranges to calculate degrees of overlap for each species
in the dataset.
library(raster) # install.packages("raster")
library(stars) # install.packages("stars")
library(ncdf4) # install.packages("ncdf4")
source("scripts/invasive_species.R")
# generate vector of taxa
inv_tax <-
c(
"Rhinella marina",
"Felis catus",
"Vulpes vulpes",
"Camelus dromedarius",
"Oryctolagus cuniculus"
)
# generate vector of taxa
inv_tax <-
c(
"Rhinella marina",
"Felis catus",
"Vulpes vulpes",
"Camelus dromedarius",
"Oryctolagus cuniculus"
)
# generate vector of tetrapod SDMs
tet_sdm <-
unique(gsub("\\..*", "", (
dir("output/SDM/Current/Shapefiles/", full.names = F)
)))
tet_sdm <-
tet_sdm[which(!tet_sdm %in% inv_tax)] # filter out invasive species
# import Australian tetrapod suitable habitat shapefiles
tet_current <-
lapply(tet_sdm, function(x)
read_sf(paste(
"output/SDM/Current/Shapefiles/", x, ".shp", sep = ""
))) %>%
bind_rows() %>%
st_transform(st_crs(oz_map))
tet_future <-
lapply(tet_sdm, function(x)
read_sf(paste(
"output/SDM/Future/Shapefiles/", x, ".shp", sep = ""
)) %>%
mutate(
Year = sapply(str_split(Year, "-"), "[[", 2)
)) %>%
bind_rows() %>%
st_transform(st_crs(oz_map))
# generate vectors with names of models and years of projection
models <-
rep(c(
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP5.85",
"SSP5.85",
"SSP5.85",
"SSP5.85"
), 3)
years <-
rep(c(
"2021-2040",
"2041-2060",
"2061-2080",
"2081-2100",
"2021-2040",
"2041-2060",
"2061-2080",
"2081-2100"
), 3)
disps <-
rep(c("Mean", "Min", "Max"), each = 8)
# import predictors
# import invasive species
inv_current <-
lapply(inv_tax, function(x)
read_sf(paste(
"output/SDM/Current/Shapefiles/", x, ".shp", sep = ""
))) %>%
bind_rows() %>%
st_transform(st_crs(oz_map))
inv_future <-
lapply(inv_tax, function(x)
read_sf(paste(
"output/SDM/Future/Shapefiles/", x, ".shp", sep = ""
)) %>%
mutate(
Year = sapply(str_split(Year, "-"), "[[", 2)
)) %>%
bind_rows() %>%
st_transform(st_crs(oz_map))
inv_current_overlap <- list()
inv_future_overlap <- list()
for (i in 1:length(inv_tax)) {
inv_current_r <-
stars::st_rasterize(inv_current %>% filter(Binomial == inv_tax[[i]]) %>% dplyr::select(geometry))
inv_current_r <-
mask(as(inv_current_r, "Raster"), as(oz_map, "Spatial"))
tet_inv_current <- list()
for (x in 1:nrow(tet_current)) {
if (
nrow(st_as_sf(
stars::st_as_stars(mask(
inv_current_r, as(tet_current[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE
)) == 0
) {
tet_inv_current[[x]] <- tibble(Binomial = tet_current[x,]$Binomial,
area_inv = 0)
} else {
tet_inv_current[[x]] <-
st_as_sf(stars::st_as_stars(mask(
inv_current_r, as(tet_current[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(area_inv = as.numeric(units::set_units(st_area(geometry), "km2")),
Binomial = tet_current[x,]$Binomial) %>%
as_tibble() %>%
dplyr::select(Binomial, area_inv)
}
writeLines(paste0("completed: ", tet_current[x,]$Binomial, " overlap with ", inv_tax[i], "; current"),
"output/progress.txt")
}
inv_current_overlap[[i]] <-
bind_rows(tet_inv_current) %>%
group_by(Binomial) %>%
summarise(area_inv = sum(area_inv)) %>%
mutate(Binomial.1 = inv_tax[[i]])
inv_future_r <-
inv_future %>% filter(Binomial == inv_tax[[i]]) %>% dplyr::select(geometry)
inv_future_r <-
lapply(1:nrow(inv_future_r), function(x)
mask(as(
stars::st_rasterize(inv_future_r[x,]), "Raster"
), as(oz_map, "Spatial")))
tet_inv_future <- list()
for (k in 1:24) {
inv.temp <- inv_future_r[[k]]
tet.temp <- tet_future %>%
filter(
Model == (inv_future %>% filter(Binomial == inv_tax[[i]]))[k,]$Model,
Year == (inv_future %>% filter(Binomial == inv_tax[[i]]))[k,]$Year,
Dispersal == (inv_future %>% filter(Binomial == inv_tax[[i]]))[k,]$Dispersal
)
tet_inv_temp <- list()
for (x in 1:nrow(tet.temp)) {
if (st_is_empty(tet.temp[x,])){
tet_inv_temp[[x]] <- tibble(
Binomial = tet.temp[x,]$Binomial,
area_inv = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
)
} else {
if (nrow(st_as_sf(
stars::st_as_stars(mask(
inv.temp, as(tet.temp[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE
)) == 0){
tet_inv_temp[[x]] <- tibble(
Binomial = tet.temp[x,]$Binomial,
area_inv = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
)
} else{
tet_inv_temp[[x]] <-
st_as_sf(stars::st_as_stars(mask(inv.temp, as(
tet.temp[x,], "Spatial"
))),
as_points = FALSE,
merge = TRUE) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(
Binomial = tet.temp[x,]$Binomial,
area_inv = as.numeric(units::set_units(st_area(geometry), "km2")),
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
) %>%
as_tibble() %>%
dplyr::select(Binomial, area_inv, Model, Year, Dispersal)
}
}
writeLines(paste0("completed: ", tet_current[x,]$Binomial, " overlap with ", inv_tax[i], "; future scenario ", k),
"output/progress.txt")
}
tet_inv_future[[k]] <-
bind_rows(tet_inv_temp) %>%
group_by(Binomial, Model, Year, Dispersal) %>%
summarise(area_inv = sum(area_inv))
}
inv_future_overlap[[i]] <- bind_rows(tet_inv_future) %>%
mutate(Year = as.numeric(Year),
Binomial.1 = inv_tax[[i]])
}
inv_current_overlap <-
inv_current_overlap %>% bind_rows() %>%
mutate(Binomial.1 = str_replace(Binomial.1, " ", "_")) %>%
spread(Binomial.1, area_inv)
inv_future_overlap <-
lapply(inv_future_overlap, bind_rows) %>%
bind_rows() %>%
mutate(Binomial.1 = str_replace(Binomial.1, " ", "_")) %>%
spread(Binomial.1, area_inv) %>%
mutate(Year = as.numeric(Year))
inv_current_overlap %>% write_csv("output/inv_current_overlap.csv")
inv_future_overlap %>% write_csv("output/inv_future_overlap.csv")
# import land-use
lu_126 <- nc_open("data/RAW/AIM-SSPRCP-LUmap-SSP1-26-v.1.0.1.nc")
lu_585 <-
nc_open("data/RAW/AIM-SSPRCP-LUmap-SSP5-Baseline-v.1.0.1.nc")
lon <- ncvar_get(lu_126, "lon")
lat <- ncvar_get(lu_126, "lat", verbose = F)
# store the data in an array
lu_126.array <- ncvar_get(lu_126)
lu_585.array <- ncvar_get(lu_585)
# extract data for 2010 (baseline)
lu_baseline <- lu_126.array[, , , 2]
# extract data for 2040, 2060, 2080, 2100
lu_126_2040 <- lu_126.array[, , , 4]
lu_126_2060 <- lu_126.array[, , , 6]
lu_126_2080 <- lu_126.array[, , , 8]
lu_126_2100 <- lu_126.array[, , , 10]
lu_585_2040 <- lu_585.array[, , , 4]
lu_585_2060 <- lu_585.array[, , , 6]
lu_585_2080 <- lu_585.array[, , , 8]
lu_585_2100 <- lu_585.array[, , , 10]
# sum "cropland_other", "cropland_bioenergy" and "built_up" categories for total human-modified area
lu.baseline <-
raster(
t(lu_baseline[, , 1] + lu_baseline[, , 2] + lu_baseline[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
)
lu.future <-
stack(
raster(
t(lu_126_2040[, , 1] + lu_126_2040[, , 2] + lu_126_2040[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_126_2060[, , 1] + lu_126_2060[, , 2] + lu_126_2060[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_126_2080[, , 1] + lu_126_2080[, , 2] + lu_126_2080[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_126_2100[, , 1] + lu_126_2100[, , 2] + lu_126_2100[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_585_2040[, , 1] + lu_585_2040[, , 2] + lu_585_2040[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_585_2060[, , 1] + lu_585_2060[, , 2] + lu_585_2060[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_585_2080[, , 1] + lu_585_2080[, , 2] + lu_585_2080[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
),
raster(
t(lu_585_2100[, , 1] + lu_585_2100[, , 2] + lu_585_2100[, , 7]),
xmn = min(lon),
xmx = max(lon),
ymn = min(lat),
ymx = max(lat),
crs = raster::crs(oz_map)
)
)
names(lu.future) <-
c(
"SSP1.26_2040",
"SSP1.26_2060",
"SSP1.26_2080",
"SSP1.26_2100",
"SSP5.85_2040",
"SSP5.85_2060",
"SSP5.85_2080",
"SSP5.85_2100"
)
lu.baseline.r <- crop(x = lu.baseline, y = as(oz_map, "Spatial"))
lu.future.r <- crop(x = lu.future, y = as(oz_map, "Spatial"))
lu.baseline.r <- mask(lu.baseline.r, as(oz_map, "Spatial"))
lu.future.r <- mask(lu.future.r, as(oz_map, "Spatial"))
# import human population
hp.baseline <-
raster("data/RAW/BaseYear_1km/baseYr_total_2000.tif") # load human population baseline
hp.future <- stack(
raster("data/RAW/SSP1_1km/ssp1_total_2040.tif"),
raster("data/RAW/SSP1_1km/ssp1_total_2060.tif"),
raster("data/RAW/SSP1_1km/ssp1_total_2080.tif"),
raster("data/RAW/SSP1_1km/ssp1_total_2100.tif"),
raster("data/RAW/SSP5_1km/ssp5_total_2040.tif"),
raster("data/RAW/SSP5_1km/ssp5_total_2060.tif"),
raster("data/RAW/SSP5_1km/ssp5_total_2080.tif"),
raster("data/RAW/SSP5_1km/ssp5_total_2100.tif")
)
names(hp.future) <-
c(
"SSP1.26_2040",
"SSP1.26_2060",
"SSP1.26_2080",
"SSP1.26_2100",
"SSP5.85_2040",
"SSP5.85_2060",
"SSP5.85_2080",
"SSP5.85_2100"
)
hp.baseline.r <- crop(x = hp.baseline, y = as(oz_map, "Spatial"))
hp.future.r <- crop(x = hp.future, y = as(oz_map, "Spatial"))
hp.baseline.r <- aggregate(hp.baseline.r, fact = 60)
hp.future.r <- aggregate(hp.future.r, fact = 60)
hp.baseline.r <- mask(hp.baseline.r, as(oz_map, "Spatial"))
hp.future.r <- mask(hp.future.r, as(oz_map, "Spatial"))
# intersect tetrapod shapefiles with land-use to get area overlapping with land-use
tet_lu_current <- list()
for (x in 1:nrow(tet_current)) {
if (
nrow(st_as_sf(
stars::st_as_stars(mask(
lu.baseline.r >= 0.5, as(tet_current[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE
)) == 0
) {
tet_lu_current[[x]] <- tibble(Binomial = tet_current[x,]$Binomial,
lu = 0)
} else {
tet_lu_current[[x]] <-
st_as_sf(stars::st_as_stars(mask(
lu.baseline.r >= 0.5, as(tet_current[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE) %>%
filter(layer == 1) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(lu = as.numeric(units::set_units(st_area(geometry), "km2")),
Binomial = tet_current[x,]$Binomial) %>%
as_tibble() %>%
dplyr::select(Binomial, lu)
}
}
tet_lu_current <-
bind_rows(tet_lu_current) %>%
group_by(Binomial) %>%
summarise(lu = sum(lu))
tet_lu_future <- list()
for (i in 1:8) {
lu.temp <- lu.future.r[[i]]
tet.temp <-
tet_future %>% filter(Model == str_split(names(lu.temp), "_")[[1]][1],
Year == str_split(names(lu.temp), "_")[[1]][2])
tet_lu_temp <- list()
for (x in 1:nrow(tet.temp)) {
if (st_is_empty(tet.temp[x,])){
tet_lu_temp[[x]] <- tibble(
Binomial = tet.temp[x,]$Binomial,
lu = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
)
} else {
if (
nrow(st_as_sf(
stars::st_as_stars(mask(
lu.temp >= 0.5, as(tet.temp[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE
)) == 0
){
tet_lu_temp[[x]] <- tibble(
Binomial = tet.temp[x,]$Binomial,
lu = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
)
} else {
tet_lu_temp[[x]] <-
st_as_sf(stars::st_as_stars(mask(
lu.temp >= 0.5, as(tet.temp[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE) %>%
filter(layer == 1) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(
Binomial = tet.temp[x,]$Binomial,
lu = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
) %>%
as_tibble() %>%
dplyr::select(Binomial, lu, Model, Year, Dispersal)
}
}
}
tet_lu_future[[i]] <- bind_rows(tet_lu_temp)
}
tet_lu_future <- bind_rows(tet_lu_future) %>%
mutate(Year = as.numeric(Year)) %>%
group_by(Binomial, Model, Year, Dispersal) %>%
summarise(lu = sum(lu))
tet_lu <- bind_rows(tet_lu_current, tet_lu_future)
# intersect tetrapod shapefiles with human population to get area overlapping with human population
tet_hp_current <- list()
for (x in 1:nrow(tet_current)) {
if (
nrow(st_as_sf(
stars::st_as_stars(mask(
hp.baseline.r >= 100, as(tet_current[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE
)) == 0
) {
tet_hp_current[[x]] <- tibble(Binomial = tet_current[x,]$Binomial,
hp = 0)
} else {
tet_hp_current[[x]] <-
st_as_sf(stars::st_as_stars(mask(
hp.baseline.r >= 100, as(tet_current[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE) %>%
filter(layer == 1) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(hp = as.numeric(units::set_units(st_area(geometry), "km2")),
Binomial = tet_current[x,]$Binomial) %>%
as_tibble() %>%
dplyr::select(Binomial, hp)
}
}
tet_hp_current <-
bind_rows(tet_hp_current) %>%
group_by(Binomial) %>%
summarise(hp = sum(hp))
tet_hp_future <- list()
for (i in 1:8) {
hp.temp <- hp.future.r[[i]]
tet.temp <-
tet_future %>% filter(Model == str_split(names(hp.temp), "_")[[1]][1],
Year == str_split(names(hp.temp), "_")[[1]][2])
tet_hp_temp <- list()
for (x in 1:nrow(tet.temp)) {
if (st_is_empty(tet.temp[x,])){
tet_hp_temp[[x]] <- tibble(
Binomial = tet.temp[x,]$Binomial,
hp = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
)
} else {
if (
nrow(st_as_sf(
stars::st_as_stars(mask(
hp.temp >= 100, as(tet.temp[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE
)) == 0
){
tet_hp_temp[[x]] <- tibble(
Binomial = tet.temp[x,]$Binomial,
hp = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
)
} else {
tet_hp_temp[[x]] <-
st_as_sf(stars::st_as_stars(mask(
hp.temp >= 100, as(tet.temp[x,], "Spatial")
)),
as_points = FALSE,
merge = TRUE) %>%
filter(layer == 1) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(
Binomial = tet.temp[x,]$Binomial,
hp = 0,
Model = tet.temp[x,]$Model,
Year = tet.temp[x,]$Year,
Dispersal = tet.temp[x,]$Dispersal
) %>%
as_tibble() %>%
dplyr::select(Binomial, hp, Model, Year, Dispersal)
}
}
}
tet_hp_future[[i]] <- bind_rows(tet_hp_temp)
}
tet_hp_future <- bind_rows(tet_hp_future) %>%
mutate(Year = as.numeric(Year)) %>%
group_by(Binomial, Model, Year, Dispersal) %>%
summarise(hp = sum(hp))
tet_hp <- bind_rows(tet_hp_current, tet_hp_future)
# calculate mean land-use over shapefile
lu.vals <- raster::extract(lu.baseline.r, tet_current)
lu.mean <- sapply(lu.vals, FUN = mean, na.rm = T)
lu.mean[which(is.na(lu.mean))] <- 0
lu.future.means <- list()
for (i in 1:8) {
empty <- NULL
lu.temp <- lu.future.r[[i]]
tet.temp <-
tet_future %>% filter(Model == str_split(names(lu.temp), "_")[[1]][[1]],
Year == str_split(names(lu.temp), "_")[[1]][[2]])
empty <- which(st_is_empty(tet.temp))
if (length(empty) > 0) {
tet.temp <- tet.temp[-empty,]
}
# fix lines
types <- vapply(sf::st_geometry(tet.temp), function(x) {
class(x)[2]
}, "")
tet.temp[which(types == "LINESTRING"),] <-
st_cast(tet.temp[which(types == "LINESTRING"),], "POLYGON")
tet.temp[which(types == "MULTILINESTRING"),] <-
st_cast(tet.temp[which(types == "MULTILINESTRING"),], "MULTIPOLYGON")
lu.vals <- raster::extract(lu.temp, tet.temp)
lu.future.means[[i]] <- sapply(lu.vals, FUN = mean, na.rm = T)
lu.future.means[[i]][which(is.na(lu.future.means[[i]]))] <- 0
if (length(empty) > 0) {
for (k in empty)
lu.future.means[[i]] <-
append(lu.future.means[[i]], 0, after = k - 1)
}
}
# calculate mean human population over shapefile
hp.vals <- raster::extract(hp.baseline.r, tet_current)
hp.mean <- sapply(hp.vals, FUN = mean, na.rm = T)
hp.mean[which(is.na(hp.mean))] <- 0
hp.future.means <- list()
for (i in 1:8) {
empty <- NULL
hp.temp <- hp.future.r[[i]]
tet.temp <-
tet_future %>% filter(Model == str_split(names(hp.temp), "_")[[1]][[1]],
Year == str_split(names(hp.temp), "_")[[1]][[2]])
empty <- which(st_is_empty(tet.temp))
if (length(empty) > 0) {
tet.temp <- tet.temp[-empty,]
}
# fix lines
types <- vapply(sf::st_geometry(tet.temp), function(x) {
class(x)[2]
}, "")
tet.temp[which(types %in% c("LINESTRING", "MULTILINESTRING")),] <-
st_cast(tet.temp[which(types %in% c("LINESTRING", "MULTILINESTRING")),], "MULTIPOLYGON")
hp.vals <- raster::extract(hp.temp, tet.temp)
hp.future.means[[i]] <- sapply(hp.vals, FUN = mean, na.rm = T)
hp.future.means[[i]][which(is.na(hp.future.means[[i]]))] <- 0
if (length(empty) > 0) {
for (k in empty)
hp.future.means[[i]] <-
append(hp.future.means[[i]], 0, after = k - 1)
}
}
# summarise
tet_sum_current <- tet_current %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(area = units::set_units(st_area(geometry), "km2")) %>%
left_join(tet_lu_current) %>%
left_join(tet_hp_current) %>%
mutate(
lu.prop = as.numeric(lu) / as.numeric(area),
hp.prop = as.numeric(hp) / as.numeric(area)
) %>%
mutate(
lu.prop = case_when(is.na(lu.prop) ~ 0, !is.na(lu.prop) ~ lu.prop),
hp.prop = case_when(is.na(hp.prop) ~ 0, !is.na(hp.prop) ~ hp.prop),
Year = 2010
) %>%
as_tibble() %>%
dplyr::select(-c(geometry, lu, hp)) %>%
add_column(lu.mean = lu.mean,
hp.mean = hp.mean)
tet_sum_future <- list()
for (i in 1:8) {
lu.temp <- lu.future[[i]]
tet.temp <-
tet_future %>% filter(Model == str_split(names(lu.temp), "_")[[1]][1],
Year == str_split(names(lu.temp), "_")[[1]][2]) %>%
mutate(Year = as.numeric(Year))
tet_sum_future[[i]] <- tet.temp %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(area = units::set_units(st_area(geometry), "km2")) %>%
left_join(tet_lu_future) %>%
left_join(tet_hp_future) %>%
mutate(
lu.prop = as.numeric(lu) / as.numeric(area),
hp.prop = as.numeric(hp) / as.numeric(area)
) %>%
mutate(
lu.prop = case_when(is.na(lu.prop) ~ 0, !is.na(lu.prop) ~ lu.prop),
hp.prop = case_when(is.na(hp.prop) ~ 0, !is.na(hp.prop) ~ hp.prop)
) %>%
as_tibble() %>%
dplyr::select(-c(geometry, lu, hp)) %>%
add_column(lu.mean = lu.future.means[[i]],
hp.mean = hp.future.means[[i]])
}
tet_sum_future <- bind_rows(tet_sum_future)
tet_sum <- bind_rows(tet_sum_current, tet_sum_future)
tet_sum <- tet_sum %>%
left_join(bind_rows(inv_current_overlap, inv_future_overlap) %>%
mutate(Year = case_when(is.na(Year) ~ 2010, TRUE ~ Year))) %>%
mutate_at(
vars(
Camelus_dromedarius,
Felis_catus,
Oryctolagus_cuniculus,
Rhinella_marina,
Vulpes_vulpes
),
funs(. / area)
)
tet_sum <- tet_sum %>%
mutate(
Camelus_dromedarius = case_when(
is.na(as.numeric(tet_sum$Camelus_dromedarius)) ~ 0,
as.numeric(tet_sum$Camelus_dromedarius) > 1 ~ 1,
TRUE ~ as.numeric(tet_sum$Camelus_dromedarius)
),
Felis_catus = case_when(
is.na(as.numeric(tet_sum$Felis_catus)) ~ 0,
as.numeric(tet_sum$Felis_catus) > 1 ~ 1,
TRUE ~ as.numeric(tet_sum$Felis_catus)),
Oryctolagus_cuniculus = case_when(
is.na(as.numeric(tet_sum$Oryctolagus_cuniculus)) ~ 0,
as.numeric(tet_sum$Oryctolagus_cuniculus) > 1 ~ 1,
TRUE ~ as.numeric(tet_sum$Oryctolagus_cuniculus)
),
Rhinella_marina = case_when(
is.na(as.numeric(tet_sum$Rhinella_marina)) ~ 0,
as.numeric(tet_sum$Rhinella_marina) > 1 ~ 1,
TRUE ~ as.numeric(tet_sum$Rhinella_marina)
),
Vulpes_vulpes = case_when(
is.na(as.numeric(tet_sum$Vulpes_vulpes)) ~ 0,
as.numeric(tet_sum$Vulpes_vulpes) > 1 ~ 1,
TRUE ~ as.numeric(tet_sum$Vulpes_vulpes)
)
)
tet_sum %>% write_csv("output/tet_lu_hp.csv")After we have generated all of the predictors we want, we use the
automated assessment procedure developed by Caetano et al. 2022 PLoS
Biol. 20, e3001544 to predict IUCN Red List assessments for Australian
tetrapods under future climate change scenarios.
First, this code chunk sources several custom-made functions to run the
automated assessment procedure as described in the Methods section of
the paper. The functions in the AA_functions.R script are
annotated.
Next, the code determines the biome for each species in the analyses -
this is because the automated assessment is performed iteratively, with
different hyperparameter tunings and feature selection in each stage and
biome partition. Then, the models are trained on current data using the
custom auto_ass function - note that hyperparameter tuning
can take a bit of time to finish.
After this, the code uses the custom pred_ass function to
generate predictions using the trained model on novel data. First on
current species (both train and test datasets) to have a complete
contemporary dataset as a baseline, and then on future projections using
the sets of projected predictor variables. This stage is really
fast.
library(glue) #install.packages("glue")
library(caret) #install.packages("caret")
library(xgboost) #install.packages("xgboost")
library(pbapply) #install.packages("pdapply")
library(pdp) #install.packages("pdp")
library(cowplot) #install.packages("cowplot")
library(DiagrammeR) #install.packages("DiagrammeR")
library(DiagrammeRsvg) #install.packages("DiagrammeRsvg")
library(rsvg) #install.packages("rsvg")
seed <- 42
# load custom functions
source("scripts/AA_functions.R")
# generate vectors with names of models and years of projection
models <-
rep(c(
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP5.85",
"SSP5.85",
"SSP5.85",
"SSP5.85"
), 3)
years <-
rep(c(
"2040",
"2060",
"2080",
"2100",
"2040",
"2060",
"2080",
"2100"
), 3)
disps <-
rep(c("Min", "Mean", "Max"), each = 8)
# generate vector of invasive taxa
inv_tax <-
c(
"Rhinella marina",
"Felis catus",
"Vulpes vulpes",
"Camelus dromedarius",
"Oryctolagus cuniculus"
)
# generate vector of tetrapod SDMs
tet_sdm <-
unique(gsub("\\..*", "", (
dir("output/SDM/Current/Shapefiles/", full.names = F)
)))
tet_sdm <-
tet_sdm[which(!tet_sdm %in% inv_tax)] # filter out invasive species
# import Australian tetrapod suitable habitat shapefiles
tet_dist <- bind_rows(
readRDS("data/mam_dist.RDS"),
readRDS("data/bird_dist.RDS"),
readRDS("data/rep_dist.RDS"),
readRDS("data/amph_dist.RDS")
) %>%
st_as_sf(crs = st_crs(readRDS("data/mam_dist.RDS")))
oz_map <-
ne_countries(country = "Australia",
scale = 10,
returnclass = "sf") %>%
st_transform(st_crs(tet_dist)) %>%
st_simplify(preserveTopology = FALSE, dTolerance = 0.01)
tet_current <-
lapply(tet_sdm, function(x)
read_sf(paste(
"output/SDM/Current/Shapefiles/", x, ".shp", sep = ""
))) %>%
bind_rows() %>%
st_transform(st_crs(tet_dist))
tet_future <-
lapply(tet_sdm, function(x)
read_sf(paste(
"output/SDM/Future/Shapefiles/", x, ".shp", sep = ""
)) %>%
mutate(
Year = sapply(str_split(Year, "-"), "[[", 2)
)) %>%
bind_rows() %>%
st_transform(st_crs(tet_dist))
# import shapefile of terresterial ecoregions
biomes <- read_sf("data/RAW/wwf_terr_ecos.shp") %>%
group_by(BIOME) %>%
summarise() %>%
st_intersection(oz_map) %>%
dplyr::select(BIOME) %>%
mutate(BIOME = as.factor(BIOME))
# find the biome most of the suitable habitat of each species falls in
tet_current %>%
st_intersection(biomes) %>%
mutate(
Year = 2010,
Model = NA,
area = units::set_units(st_area(geometry), "km2")
) %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
arrange(Binomial, desc(area)) %>%
group_by(Binomial) %>%
filter(row_number() == 1) %>%
ungroup() %>%
write_csv("output/tet_biomes.csv")
tet_current_biomes <- read_csv("output/tet_biomes.csv")
for (i in 1:24) {
tet_future %>%
filter(Year == years[[i]],
Model == models[[i]],
Dispersal == disps[[i]]) %>%
st_intersection(biomes) %>%
mutate(Year = as.numeric(Year),
area = units::set_units(st_area(geometry), "km2")) %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
arrange(Binomial, Year, Model, Dispersal, desc(area)) %>%
group_by(Binomial, Year, Model, Dispersal) %>%
filter(row_number() == 1) %>%
ungroup() %>%
write_csv(paste("output/tet_biomes_", i, ".csv", sep = ""))
}
tet_future_biomes <-
bind_rows(lapply(1:24, function(x)
read_csv(
paste("output/tet_biomes_", x, ".csv", sep = "")
)))
## import species data
# import Australian EPBC listings - these are preferred over IUCN categories when available
EPBC <-
read_csv("data/EPBC.csv") %>% mutate(Binomial = str_replace(Binomial, " ", "_"))
tet_data <- read_csv("output/tet_lu_hp.csv") %>%
left_join(read_csv("data/tet_data.csv") %>% dplyr::select(Binomial, Mass)) %>%
left_join(tet_dist %>% as_tibble() %>% dplyr::select(-geometry)) %>%
left_join(EPBC %>% dplyr::select(-Class)) %>%
left_join(bind_rows(tet_current_biomes, tet_future_biomes) %>% dplyr::select(-area)) %>%
mutate(
category = case_when(is.na(EPBC) ~ category,
TRUE ~ EPBC),
category = case_when(Year == 2010 ~ category),
threatened = case_when(
category %in% c("LC", "NT") ~ "no",
category %in% c("VU", "EN", "CR") ~ "yes"
)
) %>%
dplyr::select(-EPBC)
tet_data[which(is.na(tet_data[, "category"])), "category"] <- "NE"
# drop species with extremely wrong (3 times larger/smaller or more) predicted ranges
tet_data <- tet_data %>%
left_join(
tet_dist %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
mutate(area_iucn = as.numeric(units::set_units(
st_area(geometry), "km2"
))) %>%
as_tibble() %>%
dplyr::select(Binomial, area_iucn)
) %>%
mutate(
pred_real = log(area / area_iucn),
category = case_when(abs(pred_real) >= log(3) &
Year == 2010 ~ "NE",
TRUE ~ category),
threatened = case_when(
category %in% c("LC", "NT") ~ "no",
category %in% c("VU", "EN", "CR") ~ "yes"
)
) %>%
dplyr::select(-c(pred_real, area_iucn)) %>%
drop_na(Mass)
tet_data %>% write_csv("output/AA/tet_data_AA.csv")
# combine all data and omit species without ranges
tet_all <- bind_rows(tet_current, tet_future) %>%
filter(Binomial %in% tet_data$Binomial) %>%
mutate(
id = str_replace(paste(Binomial, Model, Dispersal, Year, sep = "_"), "NA_NA_NA", "2020"),
Model = case_when(is.na(Model) ~ "current",
TRUE ~ Model),
Year = case_when(is.na(Year) ~ "current",
TRUE ~ Year),
Dispersal = case_when(is.na(Dispersal) ~ "current",
TRUE ~ Dispersal)
)
tet_all_sp <- st_centroid(tet_all)
empty <- which(is.na(st_coordinates(tet_all_sp)))
if (length(empty) > 0) {
tet_all_sp <- tet_all_sp[-empty,]
tet_all <- tet_all[-empty,]
}
# format for xgboost
pred <- tet_data %>% filter(Year == 2010) %>% dplyr::select(-c(Binomial, Year, Model, Dispersal)) %>% as.data.frame()
rownames(pred) <- filter(tet_data, Year == 2010)$Binomial
#add class to pred
class <- read.csv("Data/tet_data.csv", row.names = 1)
pred$class <- class[rownames(pred), "Class"]
#group biome into tropical = c(1, 7), temperate = c(8, 4, 10), dry = c(12, 13)
pred$biome <- pred$BIOME
pred$biome[pred$BIOME %in% c(1, 7)] <- "TRO"
pred$biome[pred$BIOME %in% c(8, 4, 10)] <- "TEM"
pred$biome[pred$BIOME %in% c(12, 13)] <- "DRY"
#remove non-assessed species
pred_ne_dd <- pred[!pred$category %in% c("EX", "EW"),]
pred <- pred[!pred$category %in% c("NE", "DD", "EX", "EW"),]
pred_thrt <- pred[pred$threatened == "yes", ]
pred_nthrt <- pred[pred$threatened == "no", ]
#feature names
features <-
c(
"area",
"lu.prop",
"hp.prop",
"lu.mean",
"hp.mean",
"Camelus_dromedarius",
"Felis_catus",
"Oryctolagus_cuniculus",
"Rhinella_marina",
"Vulpes_vulpes",
"Mass"
)
#get contingency table
cat <- "threatened"
table(pred[, c(cat, "class", "biome")])
#create even folds
# set seed for reproducibility
folds <- even_folds(pred, cat, 5, 5, seed)
#define response variable
y <- pred[, cat]
y <- y %>% as.factor %>% {
as.numeric(.) - 1
}
#define predictor
X <- pred[, c(features, "biome", "class")]
data.table::setDT(X)
tr_dummies <- caret::dummyVars(~ ., data = X)
X <- predict(tr_dummies, X)
#run auto assessment
#threatened vs LC
aa_thrt <- auto_ass(X, y, path = "output/AA/", folds, seed)
write_rds(aa_thrt, paste("output/AA/", seed, "threatened.rds", sep = "_"))
# EN_CR vs VU
ind <-
as.logical((pred_thrt$category == "CR") + (pred_thrt$category == "EN"))
pred_thrt$en_cr <- "no"
pred_thrt$en_cr[ind] <- "yes"
#get contingency table
cat <- "en_cr"
table(pred_thrt[, c(cat, "class", "biome")])
#create even folds
# set seed for reproducibility
folds <- even_folds(pred_thrt, cat, 5, 3, seed)
#define response variable
y <- pred_thrt[, cat]
y <- y %>% as.factor %>% {
as.numeric(.) - 1
}
#define predictor
X <- pred_thrt[, c(features, "biome", "class")]
data.table::setDT(X)
tr_dummies <- caret::dummyVars(~ ., data = X)
X <- predict(tr_dummies, X)
#run auto assessment
aa_en_cr <- auto_ass(X, y, path = "output/AA/", folds, seed)
write_rds(aa_en_cr, paste("output/AA/", seed, "en_cr.rds", sep = "_"))
# NT vs LC
nt_ind <-
pred_nthrt[pred_nthrt$threatened == "no", ]$category == "NT"
pred_nthrt$nt <- "no"
pred_nthrt$nt[nt_ind] <- "yes"
#get contingency table
cat <- "nt"
table(pred_nthrt[, c(cat, "class", "biome")])
#create even folds
# set seed for reproducibility
folds <- even_folds(pred_nthrt, cat, 5, 3, seed)
#define response variable
y <- pred_nthrt[, cat]
y <- y %>% as.factor %>% {
as.numeric(.) - 1
}
#define predictor
X <- pred_nthrt[, c(features, "biome", "class")]
data.table::setDT(X)
tr_dummies <- caret::dummyVars(~ ., data = X)
X <- predict(tr_dummies, X)
#run auto assessment
aa_nt <- auto_ass(X, y, path = "output/AA/", folds, seed)
write_rds(aa_nt, paste("output/AA/", seed, "nt.rds", sep = "_"))
#model fit table
1 - mean(aa_thrt$err)
1 - mean(aa_en_cr$err)
1 - mean(aa_nt$err)
#model predictability for table of biome and class with too few samples; This tells if models based on other species can predict the extinction risk of the few species in the biome and class.
1 - colMeans(aa_thrt$err.exclude)
1 - colMeans(aa_en_cr$err.exclude)
1 - colMeans(aa_nt$err.exclude)
#predict current species
X <- pred_ne_dd[, c(features, "biome", "class")]
data.table::setDT(X)
tr_dummies <- caret::dummyVars(~ ., data = X)
X <- predict(tr_dummies, X)
ne_dd <-
pred_ass(
X,
model.threatened = aa_thrt$xgb.final,
model.nt = aa_nt$xgb.final,
model.en_cr = aa_en_cr$xgb.final
)
tet_iucn_current <- tet_data %>% filter(Year == 2010) %>%
mutate(category = ne_dd,
category = case_when(ne_dd == "EN" ~ "EN/CR",
T ~ category)) %>%
dplyr::select(Binomial, Year, Dispersal, Model, BIOME, category) %>%
left_join(read_csv("data/tet_data.csv") %>% dplyr::select(Binomial, Class))
## predict future projections
tet_iucn_future <- list()
for (k in 1:24) {
pred_future <-
tet_data %>% filter(Year == years[[k]], Model == models[[k]], Dispersal == disps[[k]], area > 0) %>% select(-c(Year, Model, Dispersal)) %>% as.data.frame()
rownames(pred_future) <-
pred_future$Binomial
pred_future <- pred_future %>% select(-Binomial)
pred_future$class <- class[rownames(pred_future), "Class"]
#group biome into tropical = c(1, 7), temperate = c(8, 4, 10), dry = c(12, 13)
pred_future$biome <- pred_future$BIOME
pred_future$biome[pred_future$BIOME %in% c(1, 7)] <- "TRO"
pred_future$biome[pred_future$BIOME %in% c(8, 4, 10)] <- "TEM"
pred_future$biome[pred_future$BIOME %in% c(12, 13)] <- "DRY"
X <- pred_future[, c(features, "biome", "class")]
# data.table::setDT(X)
tr_dummies <- caret::dummyVars(~ ., data = X)
X <- predict(tr_dummies, X)
future <-
pred_ass(
X,
model.threatened = aa_thrt$xgb.final,
model.nt = aa_nt$xgb.final,
model.en_cr = aa_en_cr$xgb.final
)
tet_iucn_future[[k]] <-
tet_data %>% filter(Year == years[[k]], Model == models[[k]], Dispersal == disps[[k]]) %>%
dplyr::select(-category) %>%
left_join(tibble(Binomial = row.names(pred_future),
category = future)) %>%
mutate(category = case_when(is.na(category) ~ "EX",
category == "EN" ~ "EN/CR",
T ~ category)) %>%
dplyr::select(Binomial, Year, Dispersal, Model, BIOME, category)
}
tet_iucn_future <-
bind_rows(tet_iucn_future) %>% left_join(read_csv("data/tet_data.csv") %>% dplyr::select(Binomial, Class))
tet_iucn <- bind_rows(tet_iucn_current, tet_iucn_future)
tet_iucn %>% write_csv("output/AA/tet_iucn.csv")
tibble(
accuracy = c(1 - mean(aa_thrt$err), 1 - mean(aa_en_cr$err), 1 - mean(aa_nt$err)),
model = c("Threatened vs. No-threatened", "EN/CR vs VU", "NT vs LC")
) %>% write_csv("output/AA/accuracy.csv")Finally, the last chunks of code perform the analyses described in the paper: examining changing trends in threatened species identity, richness, and spatial distribution. These chunks run analyses, but also export plots which together make up the figures used in the paper.
Before we can run the analyses, we have one final step of data
preparation.
In this step we run some spatial analyses to calculate how ranges are
predicted to change - how much do they contract/expand, do their
centroids shift directionally. We also generate species richness maps,
both of all species and of threatened species only - this will be
important down the line for spatial analyses of threat hotspots.
# generate vectors with names of models and years of projection
models <-
rep(c(
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP5.85",
"SSP5.85",
"SSP5.85",
"SSP5.85"
), 3)
years <-
rep(c(
"2040",
"2060",
"2080",
"2100",
"2040",
"2060",
"2080",
"2100"
), 3)
disps <-
rep(c("Min", "Mean", "Max"), each = 8)
# generate vector of taxa
inv_tax <-
c(
"Rhinella marina",
"Felis catus",
"Vulpes vulpes",
"Camelus dromedarius",
"Oryctolagus cuniculus"
)
# generate vector of tetrapod SDMs
tet_sdm <-
unique(gsub("\\..*", "", (
dir("output/SDM/Current/Shapefiles/", full.names = F)
)))
tet_sdm <-
tet_sdm[which(!tet_sdm %in% inv_tax)] # filter out invasive species
# import Australian tetrapod suitable habitat shapefiles
tet_dist <- bind_rows(
readRDS("data/mam_dist.RDS"),
readRDS("data/bird_dist.RDS"),
readRDS("data/rep_dist.RDS"),
readRDS("data/amph_dist.RDS")
)
# generate spatial objects of tetrapod current and predicted distributions
tet_current <-
lapply(tet_sdm, function(x)
read_sf(paste(
"output/SDM/Current/Shapefiles/", x, ".shp", sep = ""
))) %>%
bind_rows() %>%
st_transform(st_crs(tet_dist))
tet_future <-
lapply(tet_sdm, function(x)
read_sf(paste(
"output/SDM/Future/Shapefiles/", x, ".shp", sep = ""
)) %>%
mutate(
Year = sapply(str_split(Year, "-"), "[[", 2)
)) %>%
bind_rows() %>%
st_transform(st_crs(tet_dist))
# generate a 50*50km equal-area grid of Australia
oz_grid <- tet_dist %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
st_make_grid(cellsize = 50000) %>%
st_intersection(
oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
) %>%
st_cast("MULTIPOLYGON") %>%
st_sf() %>%
mutate(cellid = row_number())
# add Red List categories to ranges and reproject
tet_dist_iucn_current <-
tet_current %>%
left_join(tet_iucn %>%
mutate(threatened = case_when(
category %in% c("LC", "NT") ~ "no",
category %in% c("VU", "EN/CR") ~ "yes"
))) %>%
dplyr::filter(threatened == "yes",
Year == 2010) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
tet_dist_iucn_future <-
tet_future %>%
mutate(Year = as.numeric(Year)) %>%
left_join(tet_iucn %>%
mutate(threatened = case_when(
category %in% c("LC", "NT") ~ "no",
category %in% c("VU", "EN/CR") ~ "yes"
))) %>%
filter(threatened == "yes", !Year == 2010) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
# calculate threatened species richness per cell
thrt_rich_current <-
lapply(c("Mammal", "Bird", "Reptile", "Amphibian"), function(x)
oz_grid %>%
st_join(tet_dist_iucn_current %>% filter(Class == x)) %>%
mutate(overlap = ifelse(!is.na(cellid), 1, 0)) %>%
filter(!is.na(Binomial)) %>%
group_by(cellid, Class) %>%
summarize(num_species = sum(overlap)))
thrt_rich_current <- oz_grid %>%
left_join(bind_rows(thrt_rich_current) %>%
as_tibble() %>%
dplyr::select(-geometry)) %>%
mutate(num_thrt = replace_na(num_species, 0)) %>%
dplyr::select(-num_species)
thrt_rich_future <- lapply(c(1:24), function(x)
lapply(c("Mammal", "Bird", "Reptile", "Amphibian"), function(y)
oz_grid %>%
st_join(
tet_dist_iucn_future %>%
filter(Class == y,
Year == years[[x]],
Model == models[[x]],
Dispersal == disps[[x]])
) %>%
mutate(overlap = ifelse(!is.na(cellid), 1, 0)) %>%
filter(!is.na(Binomial)) %>%
group_by(cellid, Class) %>%
summarize(num_species = sum(overlap)) %>%
mutate(Year = years[[x]],
Model = models[[x]],
Dispersal = disps[[x]])))
thrt_rich_future <- bind_rows(thrt_rich_future)
thrt_rich_future <- oz_grid %>%
left_join(
thrt_rich_future %>%
as_tibble() %>%
mutate(model2 = paste(Model, Year, Dispersal, sep = "_")) %>%
dplyr::select(-c(geometry, Year, Model, Dispersal)) %>%
spread(model2, num_species)
) %>%
gather(model2, num_thrt, -c(cellid, Class, geometry)) %>%
mutate(num_thrt = replace_na(num_thrt, 0)) %>%
mutate(Year = sapply(str_split(model2, "_"), "[[", 2),
Model = sapply(str_split(model2, "_"), "[[", 1),
Dispersal = sapply(str_split(model2, "_"), "[[", 3))
thrt_rich_current %>% write_sf("output/AA/thrt_rich_current.shp")
thrt_rich_future %>% write_sf("output/AA/thrt_rich_future.shp")
# calculate total species richness per cell
rich_current <- oz_grid %>%
st_join(
tet_current %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
) %>%
mutate(overlap = ifelse(!is.na(cellid), 1, 0)) %>%
filter(!is.na(Binomial)) %>%
group_by(cellid) %>%
summarize(num_species = sum(overlap))
rich_future <- lapply(c(1:24), function(x)
oz_grid %>%
st_join(
tet_future %>%
filter(Year == years[[x]],
Model == models[[x]],
Dispersal == disps[[x]]) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
) %>%
mutate(overlap = ifelse(!is.na(cellid), 1, 0)) %>%
filter(!is.na(Binomial)) %>%
group_by(cellid) %>%
summarize(num_species = sum(overlap)) %>%
mutate(Year = years[[x]],
Model = models[[x]],
Dispersal = disps[[x]]))
rich_future <- bind_rows(rich_future)
rich_current %>% write_sf("output/AA/rich_current.shp")
rich_future %>% write_sf("output/AA/rich_future.shp")
# define function to calculate azimuth
st_azimuth2 <- function (x, y)
{
stopifnot(all(st_is(x, "POINT")))
stopifnot(all(st_is(y, "POINT")))
x = st_geometry(x)
y = st_geometry(y)
if (length(x) < length(y)) {
ids = rep(1:length(x), length.out = length(y))
x = x[ids]
}
if (length(y) < length(x)) {
ids = rep(1:length(y), length.out = length(x))
y = y[ids]
}
x_coords = st_coordinates(x)
y_coords = st_coordinates(y)
x1 = x_coords[, 1]
y1 = x_coords[, 2]
x2 = y_coords[, 1]
y2 = y_coords[, 2]
az = (180 / pi) * atan2(x2 - x1, y2 - y1)
names(az) = NULL
az[!is.na(az) & az < 0] = az[!is.na(az) & az < 0] + 360
az[x1 == x2 & y1 == y2] = NA
return(az)
}
# calculate directions of centroid range shift (Direction), distance shifted (Distance), and range contraction/expansion (Magnitude)
models2 <-
c("SSP1.26_Min",
"SSP1.26_Mean",
"SSP1.26_Max",
"SSP5.85_Min",
"SSP5.85_Mean",
"SSP5.85_Max")
tet_directions <- list()
for (i in 1:length(models2)) {
tet_temp <- tet_future %>% mutate(Model = paste(Model, Dispersal, sep = "_")) %>% filter(Model == models2[[i]]) %>% st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
tet_directions_2040 <-
tibble(
Binomial = (tet_temp %>% filter(Year == 2040))$Binomial,
Direction = st_azimuth2(
tet_current %>% st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>% st_centroid(),
tet_temp %>% filter(Year == 2040) %>% st_centroid()
),
Magnitude = as.numeric(tet_temp %>% filter(Year == 2040) %>% st_area()) /
(
tet_current %>% st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>% st_area()
),
Distance = as.numeric(
st_distance(
tet_current %>% st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>% st_centroid(),
tet_temp %>% filter(Year == 2040) %>% st_centroid(),
by_element = T
) / 1000
),
Year = 2040
)
tet_directions_2060 <-
tibble(
Binomial = (tet_temp %>% filter(Year == 2060))$Binomial,
Direction = st_azimuth2(
tet_temp %>% filter(Year == 2040) %>% st_centroid(),
tet_temp %>% filter(Year == 2060) %>% st_centroid()
),
Magnitude = as.numeric(tet_temp %>% filter(Year == 2060) %>% st_area()) /
(tet_temp %>% filter(Year == 2040) %>% st_area()),
Distance = as.numeric(
st_distance(
tet_temp %>% filter(Year == 2040) %>% st_centroid(),
tet_temp %>% filter(Year == 2060) %>% st_centroid(),
by_element = T
) / 1000
),
Year = 2060
)
tet_directions_2080 <-
tibble(
Binomial = (tet_temp %>% filter(Year == 2080))$Binomial,
Direction = st_azimuth2(
tet_temp %>% filter(Year == 2060) %>% st_centroid(),
tet_temp %>% filter(Year == 2080) %>% st_centroid()
),
Magnitude = as.numeric(tet_temp %>% filter(Year == 2080) %>% st_area()) /
(tet_temp %>% filter(Year == 2060) %>% st_area()),
Distance = as.numeric(
st_distance(
tet_temp %>% filter(Year == 2060) %>% st_centroid(),
tet_temp %>% filter(Year == 2080) %>% st_centroid(),
by_element = T
) / 1000
),
Year = 2080
)
tet_directions_2100 <-
tibble(
Binomial = (tet_temp %>% filter(Year == 2100))$Binomial,
Direction = st_azimuth2(
tet_temp %>% filter(Year == 2080) %>% st_centroid(),
tet_temp %>% filter(Year == 2100) %>% st_centroid()
),
Magnitude = as.numeric(tet_temp %>% filter(Year == 2100) %>% st_area()) /
(tet_temp %>% filter(Year == 2080) %>% st_area()),
Distance = as.numeric(
st_distance(
tet_temp %>% filter(Year == 2080) %>% st_centroid(),
tet_temp %>% filter(Year == 2100) %>% st_centroid(),
by_element = T
) / 1000
),
Year = 2100
)
tet_directions[[i]] <- bind_rows(
tet_directions_2040,
tet_directions_2060,
tet_directions_2080,
tet_directions_2100
) %>%
mutate(Model = str_split(models2[[i]], "_")[[1]][[1]],
Dispersal = str_split(models2[[i]], "_")[[1]][[2]])
}
tet_directions <-
bind_rows(tet_directions) %>% mutate(Magnitude = as.numeric(Magnitude))
tet_directions %>% write_csv("output/AA/tet_directions.csv")
tet_future2 <- bind_rows(
bind_rows(
tet_current %>% mutate(Year = "2040", Model = "SSP1.26", Dispersal = "Min"),
tet_current %>% mutate(Year = "2040", Model = "SSP1.26", Dispersal = "Mean"),
tet_current %>% mutate(Year = "2040", Model = "SSP1.26", Dispersal = "Max"),
tet_current %>% mutate(Year = "2040", Model = "SSP5.85", Dispersal = "Min"),
tet_current %>% mutate(Year = "2040", Model = "SSP5.85", Dispersal = "Mean"),
tet_current %>% mutate(Year = "2040", Model = "SSP5.85", Dispersal = "Max")
),
tet_future %>% filter(Year == "2040") %>% mutate(Year = "2060"),
tet_future %>% filter(Year == "2060") %>% mutate(Year = "2080"),
tet_future %>% filter(Year == "2080") %>% mutate(Year = "2100")
)
direct_future <- lapply(c(1:24), function(x)
oz_grid %>%
st_join(
tet_future2 %>%
left_join(tet_directions %>% mutate(Year = as.character(Year))) %>%
filter(Year == years[[x]],
Model == models[[x]],
Dispersal == disps[[x]]) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
)
) %>%
filter(!is.na(Binomial)) %>%
group_by(cellid) %>%
summarize(
Direction = mean(Direction, na.rm = T),
Magnitude = mean(Magnitude, na.rm = T),
Distance = mean(Distance, na.rm = T)
) %>%
mutate(Year = years[[x]],
Model = models[[x]],
Dispersal = disps[[x]]))
direct_future <- bind_rows(direct_future)
direct_future %>% write_sf("output/AA/direct_future.shp")Next we describe how many threatened species there are in each time-step under each scenario of climate change, and what predicts changes in threatened category. This chunk generates the plots making up Figures 1 and S5.
library(MASS) #install.packages("MASS")
library(gtsummary) #install.packages("gtsummary")
library(flextable) #install.packages("flextable")
library(ggh4x) #install.packages("ggh4x")
library(ggsankey) #devtools::install_github("davidsjoberg/ggsankey")
# generate barplot of distribution of threatened species in the future
tet_iucn <- read_csv("output/AA/tet_iucn.csv")
tet_iucn %>%
group_by(Year, Model, Dispersal, Class) %>%
count(category) %>%
filter(!Year == 2010) %>%
mutate(
category = case_when(is.na(category) ~ "EX",
TRUE ~ category),
category = factor(category, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
Dispersal = factor(Dispersal, levels = c("Min", "Mean", "Max"))
) %>%
filter(!category %in% c("LC", "NT")) %>%
ggplot(aes(x = interaction(Dispersal, as.character(Year)), y = n, fill = category)) +
geom_bar(stat = "identity") +
facet_grid(cols = vars(Model),
rows = vars(Class)) +
scale_fill_manual(values = c("brown1", "darkred", "black"),
labels = c("Vulnerable", "Endangered/Critically Endangered", "Extinct"),
name = "IUCN Category") +
scale_x_discrete(guide = "axis_nested") +
labs(x = "Dispersal / Year", y = "Number of Species") +
theme_bw() +
theme(legend.position = "bottom")
ggsave("output/tet_thrt.pdf", width = 16, height = 8, scale = .8)
# what predicts change in threat status?
tet_directions <-
read_csv("output/AA/tet_directions.csv")
# create new dataframe including all variables required for a polynomial regression
disps <- rep(c("Min", "Mean", "Max"), 2)
ssps <- rep(c("SSP1.26", "SSP5.85"), each = 3)
tet_directions_sum <-
tet_iucn %>%
mutate(
category = case_when(
category == "LC" ~ 1,
category == "NT" ~ 2,
category == "VU" ~ 3,
category == "EN/CR" ~ 4,
category == "EX" ~ 5
)
) %>%
select(-BIOME) %>%
drop_na(Dispersal) %>%
left_join(tet_directions) %>%
left_join(
tet_iucn %>%
mutate(
category = case_when(
category == "LC" ~ 1,
category == "NT" ~ 2,
category == "VU" ~ 3,
category == "EN/CR" ~ 4,
category == "EX" ~ 5
)
) %>%
filter(Year == 2010) %>%
dplyr::select(Binomial, category, Class) %>%
rename(current = category)
) %>%
group_by(Binomial, Model) %>%
arrange(Year) %>%
mutate(cat_lag = case_when(Year == 2040 ~ current,
T ~ lag(category)),
Model = paste(Model, Dispersal),
Direction = case_when(category == 5 ~ NA, T ~ Direction),
Magnitude = case_when(category == 5 ~ NA, T ~ Magnitude),
Distance = case_when(category == 5 ~ NA, T ~ Distance)) %>%
filter(cat_lag < 5) %>%
mutate(category = factor(category),
cat_lag = factor(cat_lag))
# model category ~ lag(category) + SSP * Model + Direction + Magnitude + Distance + Class
m <-
polr(
category ~ cat_lag + Model + Class,
data = tet_directions_sum,
Hess = TRUE
)
summary(m)
ctable <- coef(summary(m))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p)
ci <- confint(m)
exp(cbind(OR = coef(m), ci))
summary(update(m, method = "probit", Hess = TRUE), digits = 3)
summary(update(m, method = "logistic", Hess = TRUE), digits = 3)
polr_tab <- gtsummary::tbl_regression(m)
polr_tab %>%
gtsummary::as_flex_table() %>%
flextable::save_as_docx(path = "output/polr.docx")
tet_directions_sum %>%
ungroup() %>%
mutate(
LC = predict(m, tet_directions_sum, type = "p")[, 1],
NT = predict(m, tet_directions_sum, type = "p")[, 2],
VU = predict(m, tet_directions_sum, type = "p")[, 3],
EN = predict(m, tet_directions_sum, type = "p")[, 4],
EX = predict(m, tet_directions_sum, type = "p")[, 5]
) %>%
gather(cat, prob, c(LC, NT, VU, EN, EX)) %>%
select(Binomial, cat_lag, prob, cat, Model, Class) %>%
drop_na() %>%
ggplot(aes(
x = as.numeric(as.character(cat_lag)),
y = prob,
color = factor(cat, levels = c("LC", "NT", "VU", "EN", "EX")),
fill = factor(cat, levels = c("LC", "NT", "VU", "EN", "EX"))
)) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 3)) +
theme_bw() +
scale_x_continuous(name = "IUCN category in previous timestep",
labels = c("LC", "NT", "VU", "EN/CR")) +
scale_fill_manual(
values = c("lightgrey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
labels = c("LC", "NT", "VU", "EN/CR", "EX"),
name = "IUCN Category"
) +
scale_colour_manual(
values = c("lightgrey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
labels = c("LC", "NT", "VU", "EN/CR", "EX"),
name = "IUCN Category"
) +
ylab("Predicted Probability") +
facet_wrap(~Class,
nrow = 2)
ggsave("output/predicted_prob.pdf")
# prepare data for sankey plots
dat_ssp1 <- bind_rows(
tet_iucn %>%
filter(is.na(Model)) %>%
rename(`2020` = category) %>%
left_join(
tet_iucn %>%
filter(Model == "SSP1.26", Dispersal == "Min") %>%
select(-BIOME) %>%
pivot_wider(names_from = "Year", values_from = "category") %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`)
) %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`, Class) %>%
mutate_at(c(2:6), ~ replace_na(., "EX")) %>%
mutate(Dispersal = "Min"),
tet_iucn %>%
filter(is.na(Model)) %>%
rename(`2020` = category) %>%
left_join(
tet_iucn %>%
filter(Model == "SSP1.26", Dispersal == "Mean") %>%
select(-BIOME) %>%
pivot_wider(names_from = "Year", values_from = "category") %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`)
) %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`, Class) %>%
mutate_at(c(2:6), ~ replace_na(., "EX")) %>%
mutate(Dispersal = "Mean"),
tet_iucn %>%
filter(is.na(Model)) %>%
rename(`2020` = category) %>%
left_join(
tet_iucn %>%
filter(Model == "SSP1.26", Dispersal == "Max") %>%
select(-BIOME) %>%
pivot_wider(names_from = "Year", values_from = "category") %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`)
) %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`, Class) %>%
mutate_at(c(2:6), ~ replace_na(., "EX")) %>%
mutate(Dispersal = "Max")
)
dat_ssp5 <- bind_rows(
tet_iucn %>%
filter(is.na(Model)) %>%
rename(`2020` = category) %>%
left_join(
tet_iucn %>%
filter(Model == "SSP5.85", Dispersal == "Min") %>%
select(-BIOME) %>%
pivot_wider(names_from = "Year", values_from = "category") %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`)
) %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`, Class) %>%
mutate_at(c(2:6), ~ replace_na(., "EX")) %>%
mutate(Dispersal = "Min"),
tet_iucn %>%
filter(is.na(Model)) %>%
rename(`2020` = category) %>%
left_join(
tet_iucn %>%
filter(Model == "SSP5.85", Dispersal == "Mean") %>%
select(-BIOME) %>%
pivot_wider(names_from = "Year", values_from = "category") %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`)
) %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`, Class) %>%
mutate_at(c(2:6), ~ replace_na(., "EX")) %>%
mutate(Dispersal = "Mean"),
tet_iucn %>%
filter(is.na(Model)) %>%
rename(`2020` = category) %>%
left_join(
tet_iucn %>%
filter(Model == "SSP5.85", Dispersal == "Max") %>%
select(-BIOME) %>%
pivot_wider(names_from = "Year", values_from = "category") %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`)
) %>%
select(Binomial, `2040`, `2060`, `2080`, `2100`, Class) %>%
mutate_at(c(2:6), ~ replace_na(., "EX")) %>%
mutate(Dispersal = "Max")
)
# generate sankey plots
dat_ssp1 %>%
filter(Dispersal == "Min") %>%
make_long(`2040`, `2060`, `2080`, `2100`) %>%
mutate(Class = rep(dat_ssp1[1:dim(tet_iucn %>% filter(Year == 2010))[1], ]$Class, each = 4)) %>%
ggplot(
aes(
x = x,
next_x = next_x,
node = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
next_node = factor(next_node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
fill = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
label = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX"))
)
) +
geom_sankey(flow.alpha = 0.8) +
geom_sankey_label(size = 2, aes(color = factor(
node, levels = c("LC", "NT", "VU", "EN/CR", "EX")
))) +
theme_sankey(base_size = 16) +
scale_fill_manual(
values = c("light grey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
name = "IUCN Category"
) +
scale_colour_manual(values = c("black", "black", "black", "black", "white")) +
labs(x = "Year",
title = "Sustainability (SSP1.26), Minimal dispersal") +
theme(legend.position = "none") +
facet_wrap( ~ Class, scales = "free",
nrow = 2)
ggsave("output/ssp1_min.pdf", limitsize = F)
dat_ssp1 %>%
filter(Dispersal == "Mean") %>%
make_long(`2040`, `2060`, `2080`, `2100`) %>%
mutate(Class = rep(dat_ssp1[1:dim(tet_iucn %>% filter(Year == 2010))[1], ]$Class, each = 4)) %>%
ggplot(
aes(
x = x,
next_x = next_x,
node = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
next_node = factor(next_node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
fill = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
label = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX"))
)
) +
geom_sankey(flow.alpha = 0.8) +
geom_sankey_label(size = 2, aes(color = factor(
node, levels = c("LC", "NT", "VU", "EN/CR", "EX")
))) +
theme_sankey(base_size = 16) +
scale_fill_manual(
values = c("light grey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
name = "IUCN Category"
) +
scale_colour_manual(values = c("black", "black", "black", "black", "white")) +
labs(x = "Year",
title = "Sustainability (SSP1.26), Mean dispersal") +
theme(legend.position = "none") +
facet_wrap( ~ Class, scales = "free",
nrow = 2)
ggsave("output/ssp1_mean.pdf", limitsize = F)
dat_ssp1 %>%
filter(Dispersal == "Max") %>%
make_long(`2040`, `2060`, `2080`, `2100`) %>%
mutate(Class = rep(dat_ssp1[1:dim(tet_iucn %>% filter(Year == 2010))[1], ]$Class, each = 4)) %>%
ggplot(
aes(
x = x,
next_x = next_x,
node = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
next_node = factor(next_node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
fill = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
label = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX"))
)
) +
geom_sankey(flow.alpha = 0.8) +
geom_sankey_label(size = 2, aes(color = factor(
node, levels = c("LC", "NT", "VU", "EN/CR", "EX")
))) +
theme_sankey(base_size = 16) +
scale_fill_manual(
values = c("light grey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
name = "IUCN Category"
) +
scale_colour_manual(values = c("black", "black", "black", "black", "white")) +
labs(x = "Year",
title = "Sustainability (SSP1.26), Maximal dispersal") +
theme(legend.position = "none") +
facet_wrap( ~ Class, scales = "free",
nrow = 2)
ggsave("output/ssp1_max.pdf", limitsize = F)
dat_ssp5 %>%
filter(Dispersal == "Min") %>%
make_long(`2040`, `2060`, `2080`, `2100`) %>%
mutate(Class = rep(dat_ssp5[1:dim(tet_iucn %>% filter(Year == 2010))[1], ]$Class, each = 4)) %>%
ggplot(
aes(
x = x,
next_x = next_x,
node = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
next_node = factor(next_node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
fill = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
label = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX"))
)
) +
geom_sankey(flow.alpha = 0.8) +
geom_sankey_label(size = 2, aes(color = factor(
node, levels = c("LC", "NT", "VU", "EN/CR", "EX")
))) +
theme_sankey(base_size = 16) +
scale_fill_manual(
values = c("light grey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
name = "IUCN Category"
) +
scale_colour_manual(values = c("black", "black", "black", "black", "white")) +
labs(x = "Year",
title = "Fossil-fuelled development (SSP5.85), Minimal dispersal") +
theme(legend.position = "none") +
facet_wrap( ~ Class, scales = "free",
nrow = 2)
ggsave("output/ssp5_min.pdf", limitsize = F)
dat_ssp5 %>%
filter(Dispersal == "Mean") %>%
make_long(`2040`, `2060`, `2080`, `2100`) %>%
mutate(Class = rep(dat_ssp5[1:dim(tet_iucn %>% filter(Year == 2010))[1], ]$Class, each = 4)) %>%
ggplot(
aes(
x = x,
next_x = next_x,
node = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
next_node = factor(next_node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
fill = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
label = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX"))
)
) +
geom_sankey(flow.alpha = 0.8) +
geom_sankey_label(size = 2, aes(color = factor(
node, levels = c("LC", "NT", "VU", "EN/CR", "EX")
))) +
theme_sankey(base_size = 16) +
scale_fill_manual(
values = c("light grey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
name = "IUCN Category"
) +
scale_colour_manual(values = c("black", "black", "black", "black", "white")) +
labs(x = "Year",
title = "Fossil-fuelled development (SSP5.85), Mean dispersal") +
theme(legend.position = "none") +
facet_wrap( ~ Class, scales = "free",
nrow = 2)
ggsave("output/ssp5_mean.pdf", limitsize = F)
dat_ssp5 %>%
filter(Dispersal == "Max") %>%
make_long(`2040`, `2060`, `2080`, `2100`) %>%
mutate(Class = rep(dat_ssp5[1:dim(tet_iucn %>% filter(Year == 2010))[1], ]$Class, each = 4)) %>%
ggplot(
aes(
x = x,
next_x = next_x,
node = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
next_node = factor(next_node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
fill = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX")),
label = factor(node, levels = c("LC", "NT", "VU", "EN/CR", "EX"))
)
) +
geom_sankey(flow.alpha = 0.8) +
geom_sankey_label(size = 2, aes(color = factor(
node, levels = c("LC", "NT", "VU", "EN/CR", "EX")
))) +
theme_sankey(base_size = 16) +
scale_fill_manual(
values = c("light grey",
"darkgoldenrod1",
"brown1",
"darkred",
"black"),
name = "IUCN Category"
) +
scale_colour_manual(values = c("black", "black", "black", "black", "white")) +
labs(x = "Year",
title = "Fossil-fuelled development (SSP5.85), Maximal dispersal") +
theme(legend.position = "none") +
facet_wrap( ~ Class, scales = "free",
nrow = 2)
ggsave("output/ssp5_max.pdf", limitsize = F)Next we describe which species are predicted to go extinct and where. This chunk generates the plots making up Figures 2 & S6.
# import shapefile of terresterial ecoregions
biomes <- read_sf("data/RAW/wwf_terr_ecos.shp") %>%
group_by(BIOME) %>%
summarise() %>%
st_intersection(oz_map) %>%
dplyr::select(BIOME) %>%
mutate(BIOME = as.factor(BIOME))
# generate map of biomes
biomes %>%
filter(!BIOME == 11) %>%
ggplot() +
geom_sf(aes(fill = factor(
BIOME,
labels = c(
"Tropical and subtropical moist broadleaf forests",
"Temperate broadleaf and mixed forests",
"Tropical and subtropical grasslands, savannas, and shrublands",
"Temperate grasslands, savannas, and shrublands",
"Montane grasslands and shrublands",
"Mediterranean forests, woodlands, and scrub",
"Deserts and xeric shrublands"
)
))) +
scale_fill_manual(
values = c(
"darkgreen",
"lightgreen",
"darkolivegreen1",
"darkolivegreen3",
"grey95",
"darkolivegreen",
"burlywood1"
),
name = "Biome"
) +
theme_void() +
theme(legend.position = "none") +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggplot2::ggsave("output/biomes_map.pdf")
# characteristics of extinct species
tet_extinct <-
tet_iucn %>% filter(category == "EX") %>% dplyr::select(Binomial, Class) %>% unique()
# model extinction as function of range size, biome, and mass
glm_extinct <-
glm(
extinct ~ log10(area) * factor(BIOME) * log10(Mass),
data = tet_data %>%
filter(is.na(Model)) %>%
mutate(extinct = case_when(Binomial %in% tet_extinct$Binomial ~ 1,
T ~ 0)),
family = binomial("logit")
)
anova(glm_extinct, test = "Chi")
summary(glm_extinct)
gtsummary::tbl_regression(anova(glm_extinct, test = "Chi"))
# plot bar of distribution of range sizes of species predicted to go extinct
tet_data %>%
filter(is.na(Model)) %>%
mutate(extinct = case_when(Binomial %in% tet_extinct$Binomial ~ "Yes",
T ~ "No")) %>%
ggplot(aes(x = extinct, y = area, fill = extinct)) +
geom_boxplot() +
scale_y_continuous(trans = "log10") +
scale_fill_manual(values = c("light yellow", "black")) +
theme_bw() +
labs(x = "Species extinct by end of century?",
y = "Range Size (km2)") +
theme(legend.position = "none")
ggplot2::ggsave("output/extinct_range.pdf")
# plot barplot of distribution of species predicted to go extinct per biome
tet_data %>%
filter(is.na(Model)) %>%
mutate(extinct = case_when(Binomial %in% tet_extinct$Binomial ~ "Yes",
T ~ "No")) %>%
ggplot(aes(x = extinct, fill = factor(
BIOME,
labels = c(
"Tropical and subtropical moist broadleaf forests",
"Temperate broadleaf and mixed forests",
"Tropical and subtropical grasslands, savannas, and shrublands",
"Temperate grasslands, savannas, and shrublands",
"Montane grasslands and shrublands",
"Mediterranean forests, woodlands, and scrub",
"Deserts and xeric shrublands"
)
))) +
geom_bar(position = "fill") +
scale_fill_manual(
values = c(
"darkgreen",
"lightgreen",
"darkolivegreen1",
"darkolivegreen3",
"grey95",
"darkolivegreen",
"burlywood1"
),
name = "Biome (current)"
) +
theme_bw() +
labs(x = "Species extinct by end of century?",
y = "Proportion of species")
ggplot2::ggsave("output/extinct_biomes.pdf")
# create map extent (exclude Norfolk Island and Heard & McDonald)
map_ext <-
(oz_map %>%
st_cast("POLYGON"))[
-c(which.min((oz_map %>% st_cast("POLYGON") %>% st_centroid() %>% st_coordinates())[,2]),
which.max((oz_map %>% st_cast("POLYGON") %>% st_centroid() %>% st_coordinates())[,1])),
] %>% st_bbox
# plot map of species predicted to go extinct
future_extinct <- tet_current %>%
left_join(tet_iucn %>%
filter(category == "EX")) %>%
filter(!is.na(category)) %>%
mutate(Model = paste0(Model, ", ", Dispersal, " Dispersal")) %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
ggplot() +
geom_sf(data = oz_map,
fill = "white",
colour = "grey") +
geom_sf(fill = "black", colour = NA) +
theme_void() +
facet_grid(rows = vars(Year), cols = vars(Model)) +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggplot2::ggsave(
"output/future_extinct_tet.pdf",
future_extinct,
scale = 1.5,
dpi = 300,
width = 13,
height = 10
)The next chunk analyses shifts in threat hotspots and their overlap with protected areas, including analyses of congruence between classes. This is the main set of analyses, and will likely take a bit of time to run. This chunk generates the plots making up Figures 3 and 4.
library(cowplot) #install.packages("cowplot")
library(SpatialPack) #install.packages("SpatialPack")
library(feathers) #devtools::install_github(repo = "shandiya/feathers", ref = "main")
# generate maps of threatened species hotposts
thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(Class) %>%
mutate(num_thrt = num_thrt / max(num_thrt)) %>%
ggplot() +
geom_sf(data = oz_map, fill = "white") +
geom_sf(aes(fill = num_thrt), colour = NA) +
geom_sf(
fill = NA,
colour = "black",
size = 1.5,
data = thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(Class) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise()
) +
geom_sf(data = oz_map, fill = NA) +
facet_wrap( ~ Class, ncol = 4) +
theme_bw() +
scale_fill_gradient(low = "grey90",
high = "darkred",
name = "Proportion of threatened species") +
theme(legend.position = "bottom") +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggsave("output/threat_hotspots_2020.pdf")
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Year, Model, Class) %>%
filter(Year == "2100",
Dispersal == "Min") %>%
ungroup() %>%
group_by(Model, Class) %>%
mutate(num_thrt = num_thrt / max(num_thrt)) %>%
ggplot() +
geom_sf(data = oz_map, fill = "white") +
geom_sf(aes(fill = num_thrt), colour = NA) +
geom_sf(
fill = NA,
colour = "black",
size = 1.5,
data = thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Year, Model, Class) %>%
filter(
num_thrt > quantile(num_thrt, 0.9),
Year == "2100",
Dispersal == "Min"
) %>%
summarise()
) +
geom_sf(data = oz_map, fill = NA) +
facet_wrap(Model ~ Class, ncol = 4) +
theme_bw() +
scale_fill_gradient(low = "grey90",
high = "darkred",
name = "Proportion of threatened species") +
theme(legend.position = "none") +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggsave("output/threat_hotspots_2100_mindisp.pdf")
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Year, Model, Class) %>%
filter(Year == "2100",
Dispersal == "Mean") %>%
ungroup() %>%
group_by(Model, Class) %>%
mutate(num_thrt = num_thrt / max(num_thrt)) %>%
ggplot() +
geom_sf(data = oz_map, fill = "white") +
geom_sf(aes(fill = num_thrt), colour = NA) +
geom_sf(
fill = NA,
colour = "black",
size = 1.5,
data = thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Year, Model, Class) %>%
filter(
num_thrt > quantile(num_thrt, 0.9),
Year == "2100",
Dispersal == "Mean"
) %>%
summarise()
) +
geom_sf(data = oz_map, fill = NA) +
facet_wrap(Model ~ Class, ncol = 4) +
theme_bw() +
scale_fill_gradient(low = "grey90",
high = "darkred",
name = "Proportion of threatened species") +
theme(legend.position = "none") +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggsave("output/threat_hotspots_2100_meandisp.pdf")
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Year, Model, Class) %>%
filter(Year == "2100",
Dispersal == "Max") %>%
ungroup() %>%
group_by(Model, Class) %>%
mutate(num_thrt = num_thrt / max(num_thrt)) %>%
ggplot() +
geom_sf(data = oz_map, fill = "white") +
geom_sf(aes(fill = num_thrt), colour = NA) +
geom_sf(
fill = NA,
colour = "black",
size = 1.5,
data = thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Year, Model, Class) %>%
filter(
num_thrt > quantile(num_thrt, 0.9),
Year == "2100",
Dispersal == "Max"
) %>%
summarise()
) +
geom_sf(data = oz_map, fill = NA) +
facet_wrap(Model ~ Class, ncol = 4) +
theme_bw() +
scale_fill_gradient(low = "grey90",
high = "darkred",
name = "Proportion of threatened species") +
theme(legend.position = "none") +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggsave("output/threat_hotspots_2100_maxdisp.pdf")
plot_grid(
thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(Class) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
st_intersection() %>%
mutate(n.overlaps = factor(n.overlaps)) %>%
ggplot() +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = NA,
fill = "white"
) +
geom_sf(
aes(fill = n.overlaps),
color = NA,
size = .1,
alpha = .8
) +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = "black",
fill = NA,
size = .3
) +
theme_bw() +
theme(legend.position = "none") +
scale_fill_manual(
values = c("yellow", "red", "purple", "black"),
name = "No. of classes"
),
bind_rows(
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
ungroup() %>%
filter(
Year == "2100",
Model == "SSP1.26",
Dispersal == "Min"
) %>%
st_intersection(),
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
ungroup() %>%
filter(
Year == "2100",
Model == "SSP5.85",
Dispersal == "Min"
) %>%
st_intersection(),
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
ungroup() %>%
filter(
Year == "2100",
Model == "SSP1.26",
Dispersal == "Mean"
) %>%
st_intersection(),
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
ungroup() %>%
filter(
Year == "2100",
Model == "SSP5.85",
Dispersal == "Mean"
) %>%
st_intersection(),
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
ungroup() %>%
filter(
Year == "2100",
Model == "SSP1.26",
Dispersal == "Max"
) %>%
st_intersection(),
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
summarise() %>%
ungroup() %>%
filter(
Year == "2100",
Model == "SSP5.85",
Dispersal == "Max"
) %>%
st_intersection()
) %>%
mutate(n.overlaps = factor(n.overlaps)) %>%
ggplot() +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = NA,
fill = "white"
) +
geom_sf(
aes(fill = n.overlaps),
color = NA,
size = .1,
alpha = .8
) +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = "black",
fill = NA,
size = .3
) +
theme_void() +
facet_grid(cols = vars(Model), rows = vars(Dispersal)) +
theme(legend.position = "right") +
scale_fill_manual(
values = c("yellow", "red", "purple", "black"),
name = "No. of classes"
) +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"])),
nrow = 1
)
ggsave("output/thrt_hotspots_topten.pdf", width = 12, height = 6, scale = .8)
## Congruence in threatened species hotspots
# Congruence between classes:
### CURRENT
# calculate spatially corrected Pearson's correlation
all_cells <- thrt_rich_current %>%
filter(!is.na(Class)) %>%
spread(Class, num_thrt) %>%
mutate_at(vars(Amphibian, Bird, Mammal, Reptile), function(x)
replace_na(x, 0))
all_cells_cor <-
all_cells %>%
filter(sum(Amphibian, Bird) > 0)
cor_amph_bird <-
cor.spatial(
x = all_cells_cor$Amphibian,
y = all_cells_cor$Bird,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Amphibian, Reptile) > 0)
cor_amph_rep <-
cor.spatial(
x = all_cells_cor$Amphibian,
y = all_cells_cor$Reptile,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Amphibian, Mammal) > 0)
cor_amph_mam <-
cor.spatial(
x = all_cells_cor$Amphibian,
y = all_cells_cor$Mammal,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Bird, Mammal) > 0)
cor_bird_mam <-
cor.spatial(
x = all_cells_cor$Bird,
y = all_cells_cor$Mammal,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Bird, Reptile) > 0)
cor_bird_rep <-
cor.spatial(
x = all_cells_cor$Bird,
y = all_cells_cor$Reptile,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Mammal, Reptile) > 0)
cor_mam_rep <-
cor.spatial(
x = all_cells_cor$Mammal,
y = all_cells_cor$Reptile,
coords = st_coordinates(st_centroid(all_cells_cor))
)
cong_current <-
tibble(
comp = c(
"amph_bird",
"amph_rep",
"amph_mam",
"bird_mam",
"bird_rep",
"mam_rep"
),
year = rep("2020", 6),
model = rep(NA, 6),
dispersal = rep(NA, 6),
r = c(
as.numeric(cor_amph_bird),
as.numeric(cor_amph_rep),
as.numeric(cor_amph_mam),
as.numeric(cor_bird_mam),
as.numeric(cor_bird_rep),
as.numeric(cor_mam_rep)
),
var = c(
attributes(cor_amph_bird)[[1]],
attributes(cor_amph_rep)[[1]],
attributes(cor_amph_mam)[[1]],
attributes(cor_bird_mam)[[1]],
attributes(cor_bird_rep)[[1]],
attributes(cor_mam_rep)[[1]]
)
)
# generate vectors with names of models and years of projection
models <-
rep(c(
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP5.85",
"SSP5.85",
"SSP5.85",
"SSP5.85"
), 3)
years <-
rep(c(
"2040",
"2060",
"2080",
"2100",
"2040",
"2060",
"2080",
"2100"
), 3)
disps <-
rep(c("Min", "Mean", "Max"), each = 8)
### FUTURE
cong_future <- list()
for (i in 1:24) {
all_cells <- thrt_rich_future %>%
filter(!is.na(Class),
Year == years[[i]],
Model == models[[i]],
Dispersal == disps[[i]]) %>%
spread(Class, num_thrt) %>%
mutate_at(vars(Amphibian, Bird, Mammal, Reptile), function(x)
replace_na(x, 0))
# calculate spatially corrected Pearson's correlation
all_cells_cor <-
all_cells %>%
filter(sum(Amphibian, Bird) > 0)
cor_amph_bird <-
cor.spatial(
x = all_cells_cor$Amphibian,
y = all_cells_cor$Bird,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Amphibian, Reptile) > 0)
cor_amph_rep <-
cor.spatial(
x = all_cells_cor$Amphibian,
y = all_cells_cor$Reptile,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Amphibian, Mammal) > 0)
cor_amph_mam <-
cor.spatial(
x = all_cells_cor$Amphibian,
y = all_cells_cor$Mammal,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Bird, Mammal) > 0)
cor_bird_mam <-
cor.spatial(
x = all_cells_cor$Bird,
y = all_cells_cor$Mammal,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Bird, Reptile) > 0)
cor_bird_rep <-
cor.spatial(
x = all_cells_cor$Bird,
y = all_cells_cor$Reptile,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <-
all_cells %>%
filter(sum(Mammal, Reptile) > 0)
cor_mam_rep <-
cor.spatial(
x = all_cells_cor$Mammal,
y = all_cells_cor$Reptile,
coords = st_coordinates(st_centroid(all_cells_cor))
)
cong_future[[i]] <-
tibble(
comp = c(
"amph_bird",
"amph_rep",
"amph_mam",
"bird_mam",
"bird_rep",
"mam_rep"
),
year = rep(years[[i]], 6),
model = rep(models[[i]], 6),
dispersal = rep(disps[[i]], 6),
r = c(
as.numeric(cor_amph_bird),
as.numeric(cor_amph_rep),
as.numeric(cor_amph_mam),
as.numeric(cor_bird_mam),
as.numeric(cor_bird_rep),
as.numeric(cor_mam_rep)
),
var = c(
attributes(cor_amph_bird)[[1]],
attributes(cor_amph_rep)[[1]],
attributes(cor_amph_mam)[[1]],
attributes(cor_bird_mam)[[1]],
attributes(cor_bird_rep)[[1]],
attributes(cor_mam_rep)[[1]]
)
)
}
cong_future <- bind_rows(cong_future)
# combine to single dataframe
cong <- bind_rows(bind_rows(lapply(1:6, function(x)
cong_current %>%
mutate(
model = replace_na(
as.character(model),
c(
"SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP5.85",
"SSP5.85",
"SSP5.85"
)[x]
),
dispersal = replace_na(
as.character(dispersal),
c(
"Min",
"Mean",
"Max",
"Min",
"Mean",
"Max"
)[x]
)
))),
cong_future)
# generate lineplot of congruence through time
cong %>%
mutate(comp2 = case_when(str_detect(comp, "mam") ~ comp,
T ~ "other")) %>%
#filter(dispersal == "Mean") %>%
ggplot(aes(x = as.numeric(year),
y = r,
colour = comp)) +
geom_boxplot(aes(x = as.numeric(year), group = year)) +
geom_point(shape = 21, fill = "white") +
geom_smooth(method = "lm", alpha = .1) +
theme_bw() +
scale_color_manual(values = get_pal("bee_eater"),
name = "Comparison",
labels = c("Amphibians-Birds",
"Amphibians-Mammals",
"Amphibians-Reptiles",
"Birds-Mammals",
"Birds-Reptiles",
"Mammals-Reptiles")) +
labs(x = "Year",
y = "Spatially corrected Pearson's correlation coefficient") +
scale_x_continuous(breaks = c(2020, 2040, 2060, 2080, 2100)) +
theme(legend.position = "bottom")
ggsave("output/congruence_class.pdf")
as_tibble(TukeyHSD(aov(lm(scale(r) ~ scale(as.numeric(year)) * comp, data = cong)), "comp")$comp, rownames = "comparison") %>%
mutate_if(is.numeric, round, 3) %>%
as.data.frame(row.names = rownames(TukeyHSD(aov(lm(scale(r) ~ scale(as.numeric(year)) * comp, data = cong)), "comp")$comp)) %>%
gt::gt() %>%
gt::gtsave(filename = "cong_tukey.docx", path = "output")
aov(lm(scale(r) ~ scale(as.numeric(year)) * comp, data = cong)) %>%
gtsummary::tbl_regression() %>%
gtsummary::as_flex_table() %>%
flextable::save_as_docx(path = "output/cong_aov.docx")
summary(lm(scale(r) ~ scale(as.numeric(year)) * comp, data = cong)) %>%
gtsummary::tbl_regression(estimate_fun = purrr::partial(gtsummary::style_ratio, digits = 3),
pvalue_fun = purrr::partial(gtsummary::style_sigfig, digits = 3)) %>%
gtsummary::as_flex_table() %>%
flextable::save_as_docx(path = "output/cong_summary.docx")
emmeans::test(emmeans::emtrends(lm(scale(r) ~ scale(year) * comp, data = cong %>% mutate(year = as.numeric(year))), specs = "comp", var = "year")) %>%
as_tibble() %>%
mutate_if(is.numeric, round, 3) %>%
gt::gt() %>%
gt::gtsave(filename = "cong_slopes.docx", path = "output")
# calculate total congruence (all tetrapods) between years:
thrt_rich_all_current <- thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(cellid) %>%
summarise(num_thrt = sum(num_thrt))
thrt_rich_all_future <- thrt_rich_future %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
filter(!is.na(Class)) %>%
group_by(cellid, Year, Model, Dispersal) %>%
summarise(num_thrt = sum(num_thrt))
thrt_rich_all <- thrt_rich_all_current %>%
rename(`2020` = num_thrt) %>%
full_join(thrt_rich_all_future %>%
spread(Year, num_thrt)) %>%
filter(!is.na(Model))
models <-
c("SSP1.26",
"SSP1.26",
"SSP1.26",
"SSP5.85",
"SSP5.85",
"SSP5.85")
disps <-
c("Min",
"Mean",
"Max",
"Min",
"Mean",
"Max")
cong_years <- list()
for (i in 1:6) {
# calculate spatially corrected Pearson's correlation
all_cells_cor <- thrt_rich_all %>%
filter(Model == models[[i]],
Dispersal == disps[[i]]) %>%
rowwise() %>%
mutate(sum = sum(`2020`, `2040`)) %>%
filter(sum > 0)
filt_cells <- unique(c(which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,1])), which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,2]))))
all_cells_cor <- all_cells_cor[filt_cells,]
cor_2040 <-
cor.spatial(
x = all_cells_cor$`2020`,
y = all_cells_cor$`2040`,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <- thrt_rich_all %>%
filter(Model == models[[i]],
Dispersal == disps[[i]]) %>%
rowwise() %>%
mutate(sum = sum(`2040`, `2060`)) %>%
filter(sum > 0)
filt_cells <- unique(c(which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,1])), which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,2]))))
all_cells_cor <- all_cells_cor[filt_cells,]
cor_2060 <-
cor.spatial(
x = all_cells_cor$`2040`,
y = all_cells_cor$`2060`,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <- thrt_rich_all %>%
filter(Model == models[[i]],
Dispersal == disps[[i]]) %>%
rowwise() %>%
mutate(sum = sum(`2060`, `2080`)) %>%
filter(sum > 0)
filt_cells <- unique(c(which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,1])), which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,2]))))
all_cells_cor <- all_cells_cor[filt_cells,]
cor_2080 <-
cor.spatial(
x = all_cells_cor$`2060`,
y = all_cells_cor$`2080`,
coords = st_coordinates(st_centroid(all_cells_cor))
)
all_cells_cor <- thrt_rich_all %>%
filter(Model == models[[i]],
Dispersal == disps[[i]]) %>%
rowwise() %>%
mutate(sum = sum(`2080`, `2100`)) %>%
filter(sum > 0)
filt_cells <- unique(c(which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,1])), which(!is.na(st_coordinates(st_centroid(all_cells_cor))[,2]))))
all_cells_cor <- all_cells_cor[filt_cells,]
cor_2100 <-
cor.spatial(
x = all_cells_cor$`2080`,
y = all_cells_cor$`2100`,
coords = st_coordinates(st_centroid(all_cells_cor))
)
cong_years[[i]] <-
tibble(
comp = c("2040", "2060", "2080", "2100"),
model = models[[i]],
dispersal = disps[[i]],
r = c(
as.numeric(cor_2040),
as.numeric(cor_2060),
as.numeric(cor_2080),
as.numeric(cor_2100)
),
var = c(
attributes(cor_2040)[[1]],
attributes(cor_2060)[[1]],
attributes(cor_2080)[[1]],
attributes(cor_2100)[[1]]
)
)
}
cong_years <- bind_rows(cong_years)
# generate lineplot of congruence through time
cong_years %>%
ggplot(aes(x = as.numeric(comp), y = r)) +
geom_hline(aes(yintercept = 0),
colour = "red",
linetype = "dashed") +
geom_line() +
geom_point(shape = 21, fill = "white") +
facet_grid(rows = vars(dispersal),
cols = vars(model)) +
theme_bw() +
labs(x = "Year",
y = "Spatially corrected Pearson's correlation coefficient")
ggsave("output/congruence_years.pdf")
## Protected areas
aus_PA <- read_sf("data/RAW/CAPAD2020_terrestrial.shp") %>%
st_transform(st_crs(oz_map)) %>%
st_crop(st_bbox(oz_map)) %>%
mutate(IUCN = case_when(IUCN %in% c("NA", "NAS", "NR") ~ "Not Reported",
T ~ IUCN)) %>%
mutate(IUCN = factor(
IUCN,
levels = c("Ia", "Ib", "II", "III", "IV", "V", "VI", "Not Reported")
))
# # generate map of protected areas and threat hotspots
# plot_grid(
# aus_PA %>%
# st_transform(
# "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
# ) %>%
# ggplot() +
# geom_sf(
# data = oz_map %>%
# st_transform(
# "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
# ),
# colour = NA,
# fill = "white"
# ) +
# geom_sf(
# data = thrt_rich_current %>%
# filter(!is.na(Class)) %>%
# group_by(Class) %>%
# filter(num_thrt > quantile(num_thrt, 0.9)) %>%
# ungroup() %>%
# summarise(),
# fill = "dark red",
# colour = "black"
# ) +
# geom_sf(fill = "darkgrey", colour = NA, alpha = .5) +
# geom_sf(
# data = oz_map %>%
# st_transform(
# "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
# ),
# colour = "black",
# fill = NA,
# size = .3
# ) +
# theme_bw() +
# # scale_fill_manual(values = get_pal("oriole")) +
# theme(legend.position = "none"),
# aus_PA %>%
# st_transform(
# "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
# ) %>%
# ggplot() +
# geom_sf(
# data = oz_map %>%
# st_transform(
# "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
# ),
# colour = NA,
# fill = "white"
# ) +
# geom_sf(
# data = thrt_rich_future %>%
# filter(!is.na(Class)) %>%
# group_by(Class, Year, Model, Dispersal) %>%
# filter(num_thrt > quantile(num_thrt, 0.9)) %>%
# group_by(Year, Model, Dispersal) %>%
# summarise() %>%
# filter(
# Year == "2100"
# ),
# fill = "dark red",
# colour = "black"
# ) +
# geom_sf(fill = "darkgrey", colour = NA, alpha = .5) +
# geom_sf(
# data = oz_map %>%
# st_transform(
# "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
# ),
# colour = "black",
# fill = NA,
# size = .3
# ) +
# theme_bw() +
# # scale_fill_manual(values = get_pal("oriole")) +
# facet_grid(rows = vars(Dispersal), cols = vars(Model)) +
# theme(legend.position = "bottom"),
# ncol = 1
# )
#
# ggsave("output/PA_thrt.pdf")
# analyses of overlap of protected areas with threat hotspots
aus_PA_sum <- aus_PA %>%
group_by(IUCN) %>%
summarise()
# current
PA_overlap_current <- thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(Class) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
st_centroid() %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
st_intersects(
aus_PA_sum %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
sparse = F
)
thrt_PA_current <- bind_rows(
lapply(1:8, function(x)
thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(Class) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
ungroup() %>%
st_centroid() %>%
filter(PA_overlap_current[, x]) %>%
count(Class) %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
mutate(IUCN = aus_PA_sum$IUCN[x]))
) %>%
left_join(
thrt_rich_current %>%
filter(!is.na(Class)) %>%
group_by(Class) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
count(Class) %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
rename(total = n)
) %>%
rowwise() %>%
mutate(overlap = n / total)
# future
PA_overlap_future <- thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
st_centroid() %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
) %>%
st_intersects(
aus_PA_sum %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
sparse = F
)
thrt_PA_future <- bind_rows(
lapply(1:8, function(x)
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
ungroup() %>%
st_centroid() %>%
filter(PA_overlap_future[, x]) %>%
group_by(Year, Model, Dispersal) %>%
count(Class) %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
mutate(IUCN = aus_PA_sum$IUCN[x]))
) %>%
left_join(
thrt_rich_future %>%
filter(!is.na(Class)) %>%
group_by(Class, Year, Model, Dispersal) %>%
filter(num_thrt > quantile(num_thrt, 0.9)) %>%
group_by(Year, Model, Dispersal) %>%
count(Class) %>%
as_tibble() %>%
dplyr::select(-geometry) %>%
rename(total = n)
) %>%
rowwise() %>%
mutate(overlap = n / total)
thrt_PA_future %>%
group_by(Class, Year, Model, Dispersal) %>%
summarise(overlap = sum(overlap)) %>%
ggplot() +
geom_boxplot(aes(x = Year, y = overlap, fill = Class)) +
geom_hline(
data = thrt_PA_current %>%
group_by(Class) %>%
summarise(overlap = sum(overlap)),
aes(yintercept = overlap),
colour = "blue",
linetype = "dashed",
size = 1
) +
facet_grid(cols = vars(Class)) +
scale_fill_manual(values = get_pal("eastern_rosella")) +
theme_bw()
ggsave("output/PA_overlap.pdf", width = 12, height = 6, scale = .8)Finally, we have a code chunk that produces all of the plots making up Figures S1-S4 and S7-S12. These include several sensitivity analyses, but mostly just more visualisations of data and trends.
# plot distribution of AUC values
tet_summary <- read_csv("output/SDM/tet_summary.csv") %>% left_join(tet_iucn %>% select(Binomial, Class) %>% unique()) %>% drop_na(Class)
tet_summary %>%
ggplot(aes(x = AUC)) +
geom_density(aes(fill = Class), alpha = .8) +
geom_vline(
aes(xintercept = mean_AUC, colour = Class),
data = tet_summary %>%
group_by(Class) %>%
summarise(mean_AUC = mean(AUC)),
linetype = "dashed",
size = 1
) +
theme_bw() +
facet_wrap( ~ Class, ncol = 1) +
scale_fill_manual(values = get_pal("eastern_rosella"), name = "Class") +
scale_colour_manual(values = get_pal("eastern_rosella"), name = "Class")
ggsave("output/auc_summary.pdf")
# generate map of shifts in range centroids
p_future_direct <- direct_future %>%
ggplot() +
geom_sf(aes(fill = Direction), colour = NA) +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = "black",
fill = NA
) +
scale_fill_gradientn(
colors = c(
"#00FEFF",
"#00FF6D",
"#0AFF00",
"#7BFF00",
"#EDFF00",
"#FF8400",
"#FF0100",
"#FF0092",
"#F500FF",
"#8400FF",
"#1200FF",
"#007BFF",
"#00FEFF"
),
name = "Direction Shifted",
limits = c(0, 360)
) +
theme_void() +
theme(legend.position = "bottom") +
labs(x = "Longitude",
y = "Latitude") +
facet_grid(rows = vars(Year),
cols = vars(Model)) +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
ggsave(
"output/AA/direction_iucn_tet_maps.pdf",
p_future_direct,
scale = 1.5,
dpi = 300,
width = 13,
height = 10
)
# generate map of range expansion/contraction
p_future_magnitude <- direct_future %>%
ggplot() +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = "black",
fill = "white"
) +
geom_sf(aes(fill = log(Magnitude)), colour = NA) +
geom_sf(
data = oz_map %>%
st_transform(
"+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=61.55,-10.87,-40.19,-39.4924,-32.7221,-32.8979,-9.99400000001316 +units=m +no_defs +type=crs"
),
colour = "black",
fill = NA
) +
theme_void() +
theme(legend.position = "bottom",
legend.text = element_text(angle = 45, vjust = 0)) +
scale_fill_gradient2(
low = "red",
mid = "white",
high = "blue",
name = "Change in Range Size"
) +
labs(x = "Longitude",
y = "Latitude") +
facet_grid(rows = vars(Year),
cols = vars(Model)) +
coord_sf(xlim = c(map_ext["xmin"], map_ext["xmax"]),
ylim = c(map_ext["ymin"], map_ext["ymax"]))
# plot distributions of range contractions/expansions
ggsave(
"output/AA/magnitude_iucn_tet_maps.pdf",
p_future_magnitude,
scale = 1.5,
dpi = 300,
width = 13,
height = 10
)
tet_directions_sum %>%
drop_na(Class) %>%
mutate(Magnitude = log(Magnitude)) %>%
ggplot(aes(x = Magnitude, fill = Class)) +
geom_density() +
geom_vline(
aes(xintercept = mag_mean, colour = type),
linetype = "dashed",
data = tet_directions_sum %>%
drop_na(Class) %>%
mutate(Magnitude = log(Magnitude)) %>%
group_by(Class, Year) %>%
summarise(mag_mean = mean(Magnitude, na.rm = T)) %>%
mutate(type = case_when(
mag_mean < 0 ~ "Contraction",
mag_mean > 0 ~ "Expansion"
))
) +
facet_wrap(Class ~ Year, scales = "free_y") +
coord_cartesian(xlim = c(-2.5, 2.5)) +
theme_bw() +
scale_fill_manual(values = get_pal("eastern_rosella")) +
scale_colour_manual(values = c("red", "blue"), name = "Mean range change") +
labs(x = "Range contraction/expansion")
ggplot2::ggsave("output/ranges_magnitude.pdf",
width = 13)
tet_sum <- read_csv("output/tet_lu_hp.csv")
# plot distributions of overlap with human population density and modified land
pdf("output/lu_prop.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = lu.prop)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(limits = c(0, 0.25), position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = lu.prop, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 0.25)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat >= 50% human-modified"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/lu_mean.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = lu.mean)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(limits = c(0, 1), position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = lu.mean, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 1)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Average proportion of human-modified suitable habitat"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/hp_prop.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = hp.prop)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(limits = c(0, 0.25), position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = hp.prop, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 0.25)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat with human population >= 100/km2"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/hp_mean.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = hp.mean)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(limits = c(0, 200), position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = hp.mean, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 200)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Average human population /km2"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/overlap_camel.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = Camelus_dromedarius)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = Camelus_dromedarius, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 1)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat overlapping with camels"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/overlap_cat.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = Felis_catus)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = Felis_catus, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 1)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat overlapping with cats"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/overlap_rabbit.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = Oryctolagus_cuniculus)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = Oryctolagus_cuniculus, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 1)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat overlapping with European rabbits"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/overlap_toad.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = Rhinella_marina)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = Rhinella_marina, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 1)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat overlapping with cane toads"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
pdf("output/overlap_fox.pdf", useDingbats = F)
cowplot::plot_grid(
tet_sum %>%
filter(is.na(Model)) %>%
ggplot(aes(x = Vulpes_vulpes)) +
geom_density(alpha = .75, fill = "black") +
theme_bw() +
scale_x_continuous(position = "top") +
labs(x = "Background"),
tet_sum %>%
filter(!is.na(Model)) %>%
ggplot(aes(x = Vulpes_vulpes, fill = Model)) +
geom_density(alpha = .75) +
facet_grid(rows = vars(Year), cols = vars(Dispersal)) +
theme_bw() +
scale_x_continuous(limits = c(0, 1)) +
scale_fill_manual(
values = c("yellow", "red")
) +
theme(legend.position = "bottom") +
labs(x = "Proportion of suitable habitat overlapping with red foxes"),
ncol = 1,
rel_heights = c(1, 4)
)
dev.off()
# generate barplot and table of feature importance
features <-
bind_rows(lapply(dir("output/AA", "importance", full.names = T), function(x)
read_csv(x) %>% mutate(Model = x))) %>%
drop_na() %>%
# filter(!Feature %in% c("X", "Y")) %>%
mutate(
Type = case_when(
Feature %in% c("area", "Mass") ~ "Species traits",
Feature %in% c("hp.mean", "hp.prop", "lu.mean", "lu.prop") ~ "Habitat",
Feature %in% c(
"Rhinella_marina",
"Felis_catus",
"Oryctolagus_cuniculus",
"Vulpes_vulpes",
"Camelus_dromedarius"
) ~ "Introduced species",
Feature %in% c("biomeDRY", "biomeTEM", "biomeTRO") ~ "Geography",
Feature %in% c("classReptile", "classAmphibian", "classMammal", "classBird") ~ "Taxonomy"
),
Feature = case_when(
Feature == "area" ~ "Range Size",
Feature == "hp.mean" ~ "Mean HPD",
Feature == "hp.prop" ~ "% HPD",
Feature == "lu.mean" ~ "Mean HML",
Feature == "lu.prop" ~ "% HML",
Feature == "Rhinella_marina" ~ "Cane toads",
Feature == "Felis_catus" ~ "Cats",
Feature == "Oryctolagus_cuniculus" ~ "European rabbits",
Feature == "Vulpes_vulpes" ~ "Red foxes",
Feature == "Camelus_dromedarius" ~ "Dromedary camels",
Feature == "biomeDRY" ~ "Dry biomes",
Feature == "biomeTEM" ~ "Temperate biomes",
Feature == "biomeTRO" ~ "Tropical biomes",
Feature == "classReptile" ~ "Reptiles",
Feature == "classMammal" ~ "Mammals",
Feature == "classAmphibian" ~ "Amphibians",
Feature == "classBird" ~ "Birds",
T ~ Feature
),
Model = case_when(
Model == "output/AA/42en_cr_importance.csv" ~ "EN/CR vs VU",
Model == "output/AA/42nt_importance.csv" ~ "NT vs LC",
Model == "output/AA/42threatened_importance.csv" ~ "Threatened vs Non-threatened"
)
)
plot_grid(
features %>%
ggplot(aes(
x = fct_reorder(Feature, Gain, .fun = median, .desc = T),
y = Gain,
fill = Type
)) +
geom_boxplot() +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
) +
labs(x = "Feature") +
scale_fill_manual(values = get_pal("superb_fairy_wren")),
features %>%
ggplot(aes(
x = fct_reorder(Type, Gain, .fun = median, .desc = T),
y = Gain,
fill = Type,
)) +
geom_boxplot() +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
) +
labs(x = "Feature type") +
scale_fill_manual(values = get_pal("superb_fairy_wren")),
ncol = 1,
labels = "AUTO"
)
ggplot2::ggsave("output/feature_importance.pdf", limitsize = F)
features %>%
dplyr::select(Feature, Type, Model, Gain) %>%
group_by(Feature) %>%
mutate(mean_gain = mean(Gain)) %>%
ungroup() %>%
arrange(desc(mean_gain), desc(Gain)) %>%
mutate(Gain = round(Gain, 3)) %>%
dplyr::select(-mean_gain) %>%
write_csv("output/feature_importance.csv")
as_tibble(TukeyHSD(aov(Gain ~ Type, features))$Type, rownames = "comparison") %>%
mutate_if(is.numeric, round, 3) %>%
as.data.frame(row.names = rownames(TTukeyHSD(aov(Gain ~ Type, features))$Type)) %>%
gt::gt() %>%
gt::gtsave(filename = "features_tukey.docx", path = "output")
# run sensitivity analyses on omitted species
drop_sp <-
read_csv("output/AA/tet_data_AA.csv") %>%
filter(is.na(Model) & category %in% c("DD", "NE")) %>%
pull(Binomial)
tet_omitted <- tet_dist %>%
left_join(read_csv("output/SDM/tet_summary.csv")) %>%
filter(Binomial %in%
drop_sp |
is.na(threshold)) %>%
as_tibble() %>%
dplyr::select(Binomial, category) %>%
mutate(
class = case_when(
Binomial %in% readRDS("data/mam_dist.RDS")$Binomial ~ "Mammal",
Binomial %in% readRDS("data/bird_dist.RDS")$Binomial ~ "Bird",
Binomial %in% readRDS("data/rep_dist.RDS")$Binomial ~ "Reptile",
Binomial %in% readRDS("data/amph_dist.RDS")$Binomial ~ "Amphibian"
)
)
bind_rows(
tet_omitted %>%
count(class, category) %>%
left_join(tet_omitted %>% count(class) %>% rename(total = n)) %>%
mutate(n = n / total) %>%
mutate(
category = case_when(is.na(category) ~ "NE",
T ~ category),
category = factor(category,
levels = c("LC", "NT", "VU", "EN", "CR", "DD", "NE")),
dat = "Not analysed"
),
read_csv("output/SDM/tet_summary.csv") %>%
select(Binomial) %>%
mutate(
class = case_when(
Binomial %in% readRDS("data/mam_dist.RDS")$Binomial ~ "Mammal",
Binomial %in% readRDS("data/bird_dist.RDS")$Binomial ~ "Bird",
Binomial %in% readRDS("data/rep_dist.RDS")$Binomial ~ "Reptile",
Binomial %in% readRDS("data/amph_dist.RDS")$Binomial ~ "Amphibian"
)
) %>%
left_join(tet_dist %>%
as_tibble() %>%
select(Binomial, category)) %>%
count(class, category) %>%
left_join(
read_csv("output/SDM/tet_summary.csv") %>%
select(Binomial) %>%
mutate(
class = case_when(
Binomial %in% readRDS("data/mam_dist.RDS")$Binomial ~ "Mammal",
Binomial %in% readRDS("data/bird_dist.RDS")$Binomial ~ "Bird",
Binomial %in% readRDS("data/rep_dist.RDS")$Binomial ~ "Reptile",
Binomial %in% readRDS("data/amph_dist.RDS")$Binomial ~ "Amphibian"
)
) %>%
left_join(tet_dist %>%
as_tibble() %>%
select(Binomial, category)) %>% count(class) %>% rename(total = n)
) %>%
mutate(n = n / total) %>%
mutate(
category = case_when(is.na(category) ~ "NE",
T ~ category),
category = factor(category,
levels = c("LC", "NT", "VU", "EN", "CR", "DD", "NE")),
dat = "Analysed"
)
) %>%
ggplot(aes(x = class, y = n, fill = category)) +
geom_bar(stat = "identity", position = "stack") +
scale_fill_manual(
values = c(
"lightgrey",
"darkgoldenrod1",
"brown1",
"brown3",
"darkred",
"grey15",
"grey5"
),
name = "IUCN Category"
) +
theme_bw() +
labs(x = "Class", y = "Proportion of species") +
theme(legend.position = "bottom") +
facet_wrap( ~ dat)
ggplot2::ggsave("output/sensitivity_iucn.pdf")
tet_dist %>%
mutate(area = units::set_units(st_area(geometry), "km^2")) %>%
left_join(read_csv("output/SDM/tet_summary.csv")) %>%
mutate(analysed = case_when(is.na(threshold) ~ F,
T ~ T)) %>%
ggplot(aes(x = analysed, y = as.numeric(area), fill = analysed)) +
geom_boxplot() +
scale_y_continuous(trans = "log10") +
theme_bw() +
scale_fill_manual(values = c("grey", "light blue")) +
labs(x = "Included in Analysis",
y = "Range Size (km2)") +
theme(legend.position = "none")
ggplot2::ggsave("output/sensitivity_range.pdf")