# R scripts for short course "Species distribution models with Bayesian statistics in R"
# A. Marcia Barbosa (https://modtools.wordpress.com)


# SETUP ####

# load required packages:

library(terra)


# set working directory:

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


# import required data:

# (you can skip this section if you just ran the entire script 1a for the same species and region)

# example terrestrial species:
my_species <- "Loxodonta africana"  # if you choose another species, change it here, not all over the script!

# example marine species:
# my_species <- "Phocoena sinus"

output_data_folder <- paste0("../outputs/", my_species, "/occ_data")

download_window <- read.csv(paste0(output_data_folder, "/download_window.csv"))[,1]

presences <- read.csv(paste0(output_data_folder, "/presences_all.csv"))

presences <- terra::vect(presences, geom = c("decimalLongitude", "decimalLatitude"), keepgeom = TRUE, crs = "EPSG:4326")

countries <- geodata::world(path = "../outputs")

layers <- geodata::worldclim_global(var = "bio", res = 5, path = "../outputs")

# or for MARINE layers:
# marine_layer_names <- c("Chlorophyll", "Current.Velocity", "Ice.thickness", "Salinity", "Temperature")
# layers <- terra::rast(lapply(marine_layer_names, function(v) geodata::bio_oracle(path = "../outputs", var = v, stat = "Mean", benthic = FALSE)))

# simplify variable names by replacing the prefix with nothing:
names(layers) <- gsub("wc2.1_5m_|Present.", "", names(layers))  # for the terrestrial example data


# DELIMIT MODELLING REGION ####

plot(presences, background = "lightblue", ext = terra::ext(presences) + 2)

# add countries map:
plot(countries, col = "sienna4", alpha = 0.3, add = TRUE)

# select a modelling region:
pres_buff <- fuzzySim::getRegion(presences, type = "width")  # see help file and EXPLORE OTHER OPTIONS! this choice will strongly influence the modelling results

plot(pres_buff, col = "yellow", border = NA, background = "lightblue")
plot(presences, cex = 0.3, add = TRUE)
plot(countries, col = "sienna4", alpha = 0.3, lwd = 0.1, add = TRUE)

mod_region <- pres_buff

# IF YOU USED A LIMITED WINDOW OF COORDINATES to download the occurrence data, or if your occurrence data come from a limited surveyed region, you must limit the modelling region also to that spatial area from which data were obtained, e.g.:
mod_region <- terra::crop(mod_region, download_window)

plot(mod_region, border = "darkred", lwd = 2, add = TRUE)

# (you can also use 'crop' to limit the modelling region to a smaller window of your own preference)
# this is actually preferable during the course, otherwise your models and (especially) subsequent analyses may take too long to compute here!
my_window <- terra::ext(10, 45, -38, 0)  # edit for your own region!
# for the marine example:
# my_window <- terra::ext(-118, -112, 28, 34)
mod_region <- terra::crop(mod_region, my_window)
plot(mod_region, border = "darkblue", lwd = 3, add = TRUE)


# you can also limit the modelling region to user-defined polygons, e.g. (IF YOUR SPECIES IS NOT MARINE!) to apparently well-surveyed countries:
countries <- terra::disagg(countries)  # disaggregate into separate polygons for mainlands and islands
countries_with_points <- countries[presences, ]
countries_with_points$NAME_0
plot(countries_with_points, col = "tan", background = "lightblue")
plot(countries, col = "tan", alpha = 0.5, add = TRUE)
plot(presences, col = "darkblue", cex = 0.3, add = TRUE)
text(countries_with_points, countries_with_points$NAME_0, cex = 0.6, halo = TRUE)
# judge if some countries are visibly better surveyed (in this dataset) for this species; compare with occurrence data from other sources, e.g. https://www.iucnredlist.org, national atlases, data papers
# so, assuming we can't complement the dataset with occurrences from other sources, we can select the apparently better surveyed countries for modelling:
include <- c("Namibia", "Botswana", "South Africa", "Zimbabwe", "Zambia")
countries_include <- subset(countries_with_points, countries_with_points$NAME_0 %in% include)
mod_region <- terra::crop(mod_region, countries_include)

# you can (instead or additionally) select polygons to exclude due to apparently insufficient survey:
exclude <- c("Angola", "Mozambique")
# mod_region <- terra::erase(mod_region, subset(countries_with_points, countries_with_points$NAME_0 %in% exclude))
mod_region <- crop(subset(countries, !(countries$NAME_0 %in% exclude)), mod_region)

# plot the modelling region:
plot(pres_buff, col = "yellow", border = NA, background = "lightblue")
# plot(terra::ext(download_window), border = "darkred", lwd = 2, add = TRUE)  # if applicable
plot(my_window, border = "darkred", add = TRUE)
plot(presences, cex = 0.3, add = TRUE)
plot(countries, col = "sienna4", lwd = 0.1, alpha = 0.3, add = TRUE)
plot(mod_region, col = "blue", alpha = 0.3, add = TRUE)


# now import and cut (crop + mask) the variable maps to the extent of the modelling region:
layers_cut <- terra::crop(layers, mod_region, snap = "out", mask = TRUE)
plot(layers_cut[[1:4]])

plot(layers_cut[[1]], main = my_species)
plot(countries, add = TRUE)
plot(mod_region, border = "blue", lwd = 2, add = TRUE)
plot(presences, pch = 20, cex = 0.3, col = "magenta", add = TRUE)
# check that everything overlaps correctly!


#### ASSIGNMENT 2: if you've tried different data than the example, post the map that you got with the above commands


# SET SPATIAL RESOLUTION ####

terra::res(layers_cut)  # this is the spatial resolution in map units, which are geographical degrees in this CRS
# but this can mean very different pixel sizes across latitudes!
# see e.g. https://modtools.wordpress.com/2024/01/23/actual-pixel-sizes-of-unprojected-raster-maps

# check mean pixel size in this study region:
source("https://raw.githubusercontent.com/AMBarbosa/ecotrends/refs/heads/master/R/pixelArea.R")
pixelArea(layers_cut, unit = "km")

# if necessary (e.g. to make computation times more reasonable, or to model more coarsely gridded occurrence data), you can aggregate the raster layers, to e.g. a 3-times coarser resolution:
layers_aggr <- terra::aggregate(layers_cut, fact = 3, fun = mean, na.rm = TRUE)
terra::res(layers_aggr)
dim(layers_aggr)
plot(layers_aggr[[1]])
plot(presences, cex = 0.3, add = TRUE)

# run the command below ONLY IF YOU DO NEED / DECIDE to aggregate the layers (NOT FOR THE MARINE EXAMPLE DATA). ASK ME if in doubt!
layers_cut <- layers_aggr

pixelArea(layers_cut, unit = "km")


# save the cut layers to a folder on disk (we'll need them later!):

layers_folder <- paste0("../outputs/", my_species, "/layers_cut")
dir.create(layers_folder, showWarnings = FALSE)
terra::writeRaster(layers_cut, filename = paste0(layers_folder, "/layers_cut.tif"), gdal = c("COMPRESS=DEFLATE"), overwrite = TRUE)


# GET PRESENCE/(PSEUDO)ABSENCE DATASET FOR MODELLING ####

# now make a dataframe of the species occurrence data gridded to the resolution of the raster variables
# i.e., one row per pixel with the values of the variables and the presence/absence of species records:

head(presences)
gridded_data <- fuzzySim::gridRecords(rst = layers_cut, pres.coords = presences, plot = TRUE)
# if you have (enough) absences, you can provide them as 'abs.coords' to this function, instead of using all of the pixels without presence

head(gridded_data)

nrow(gridded_data)  # (if no NAs) should be the same as:
sum(!is.na(values(layers_cut[[1]])))

# NOTE: to keep computation times reasonable DURING THE COURSE, this number should not be (much) larger than 8000; consider making your region extent smaller, or your spatial resolution coarser!

# you can also select a subset of absence pixels in a way that addresses survey bias:
bias_layer <- geodata::footprint(path = "../outputs")  # TERRESTRIAL ONLY! for marine species, see e.g. https://oceanhealthindex.org/resources/data/cumulative-human-impacts
bias_layer <- terra::project(bias_layer, layers_cut)
plot(bias_layer, main = "Human footprint")
plot(presences, cex = 0.3, col = "magenta", add = TRUE)

# you may define regions to exclude from pseudo-absence generation, e.g. (FOR THE TERRESTRIAL SPECIES EXAMPLE):
unsurveyed <- subset(countries, countries$NAME_0 %in% c("Lesotho", "Eswatini"))  # you can (additionally or instead) use your own polygon(s), imported with terra::vect()
bias_layer[unsurveyed] <- NA

# select a sample of the pixels with absence of records:
gridded_data$selected <- fuzzySim::selectAbsences(gridded_data, sp.cols = "presence", coord.cols = c("x", "y"), mult.p = 10, seed = 1234, df = FALSE, bias = bias_layer)
table(gridded_data$selected)

# map the gridded data and variables to check everything looks OK:
plot(gridded_data, 5, cex = 0.4, type = "continuous")  # coloured according to the 5th column
plot(gridded_data, "presence", cex = 0.4)  # coloured according to the "presence" column
plot(subset(gridded_data, gridded_data$selected == 1), "presence", cex = 0.4, main = my_species)
add_legend("topleft", legend = paste(c("N.abs =", "N.pres ="), table(data.frame(gridded_data)[gridded_data$selected, "presence"])), cex = 0.6, bg = "white", y.intersp = 0.5)


#### ASSIGNMENT 3: if you're using different data than the example, post the map that you got with the above commands


# save gridded data to disk:

write.csv(data.frame(gridded_data), paste0(output_data_folder, "/presences_gridded.csv"), row.names = FALSE)