## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 7,
  dpi = 100,
  out.width = "95%"
)

## ----get_ready, results='hide', message=FALSE, warning=FALSE------------------
# Load packages
library(nicheR)
library(terra)

# Current directory
getwd()

# Define new directory
#setwd("YOUR/DIRECTORY")  # modify if setting a new directory

# Saving original plotting parameters
original_par <- par(no.readonly = TRUE)

## ----data---------------------------------------------------------------------
# Reference niche
data("ref_ellipse", package = "nicheR")

# Background data
data("back_data", package = "nicheR")

# Raster layers for predictions
ma_bios <- terra::rast(system.file("extdata", "ma_bios.tif",
                                   package = "nicheR"))

## ----data_check---------------------------------------------------------------
# Check reference niche
print(ref_ellipse)

# Check background data
head(back_data)

# Check the raster layers
ma_bios

## ----data_plot----------------------------------------------------------------
# Pick the variables for the background data
vars <- c("bio_1", "bio_12")

# Plotting the background data to visualize the environmental space
mars <- c(4, 4, 2, 1)
par(mar = mars)  # adjust margins for better visualization

plot_ellipsoid(ref_ellipse, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797", lwd = 2,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Reference Niche and Background Environmental Space")

## ----random-------------------------------------------------------------------
# Simulating the community
rand_comm <- random_ellipses(object = ref_ellipse,
                             background = back_data[, vars],
                             n = 20)  # number of species in the community

# check the a few details from the generated community
names(rand_comm)  # elements in the community object

print(rand_comm)  # a summary of the elements in the community object

## ----random_plot--------------------------------------------------------------
# Plotting the community
par(mar = mars)  # adjust margins for better visualization

plot_community(rand_comm, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Community of Random Ellipses")

## ----random_den---------------------------------------------------------------
# Simulating the community with the full background
rand_comm_full <- random_ellipses(object = ref_ellipse,
                                  background = back_data[, vars], n = 25)

# Simulating the community with a thinned background
rand_comm_thin <- random_ellipses(object = ref_ellipse,
                                  background = back_data[, vars], n = 25,
                                  thin_background = TRUE, resolution = 20)

## ----random_den_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"------
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with the full background
plot_community(rand_comm_full, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community from Full Background")

## Plotting the community with the thinned background
plot_community(rand_comm_thin, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community from Thinned Background")

## ----random_prop--------------------------------------------------------------
# Community with both arguments small
rand_comm1 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.1,
                              largest_proportion = 0.5)

# Community with both arguments large
rand_comm2 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.7,
                              largest_proportion = 1.5)

# Community with smallest_proportion large and largest_proportion small
rand_comm3 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.5,
                              largest_proportion = 0.7)

# Community with smallest_proportion small and largest_proportion large
rand_comm4 <- random_ellipses(object = ref_ellipse,
                              background = back_data[, vars], n = 25,
                              thin_background = TRUE, resolution = 20,
                              smallest_proportion = 0.1,
                              largest_proportion = 1.5)

## ----random_prop_plot---------------------------------------------------------
# Plotting the communities
par(mfrow = c(2, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Community with small proportions
plot_community(rand_comm1, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Small Proportions")

## Community with large proportions
plot_community(rand_comm2, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Large Proportions")

## Community with large smallest_proportion and small largest_proportion
plot_community(rand_comm3, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Large and Small")

## Community with small smallest_proportion and large largest_proportion
plot_community(rand_comm4, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Random Community with Small and Large")

## ----nested-------------------------------------------------------------------
# Simulating the community
nest_comm <- nested_ellipses(object = ref_ellipse, n = 20)

# check the a few details from the generated community
print(nest_comm)  # a summary of the elements in the community object

## ----nested_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"----------
# Plotting the community
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community of nested ellipses
plot_community(nest_comm,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Community of Nested Ellipses")

## Plotting the community of nested ellipses with the background
plot_community(nest_comm, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested Community with Background")

## ----nested_prop--------------------------------------------------------------
# Simulating the community with a small smallest_proportion
nest_comm_small <- nested_ellipses(object = ref_ellipse, n = 20,
                                   smallest_proportion = 0.1)

# Simulating the community with a large smallest_proportion
nest_comm_large <- nested_ellipses(object = ref_ellipse, n = 20,
                                   smallest_proportion = 0.6)

# Lets check the volume stats for the communities for comparison
## Communities with small smallest_proportion
mean(sapply(nest_comm_small$ellipse_community, function(x) x$volume))

## Communities with large smallest_proportion
mean(sapply(nest_comm_large$ellipse_community, function(x) x$volume))

## ----nested_prop_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"-----
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with a small smallest_proportion
plot_community(nest_comm_small,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Small Proportion")

## Plotting the community with a large smallest_proportion
plot_community(nest_comm_large,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Large Proportion")

## ----nested_bias--------------------------------------------------------------
# Simulating the community with a bias towards the border
nest_comm_small_bias <- nested_ellipses(object = ref_ellipse, n = 20,
                                        smallest_proportion = 0.1, bias = 0.2)

# Simulating the community with a bias towards the centroid
nest_comm_large_bias <- nested_ellipses(object = ref_ellipse, n = 20,
                                        smallest_proportion = 0.1, bias = 2)

# Lets check the volume stats for the communities for comparison
## Community with bias towards the border
mean(sapply(nest_comm_small_bias$ellipse_community, function(x) x$volume))

## Community with bias towards the centroid
mean(sapply(nest_comm_large_bias$ellipse_community, function(x) x$volume))

## ----nested_bias_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"-----
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with a bias towards the border
plot_community(nest_comm_small_bias,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Bias Towards Border")

## Plotting the community with a bias towards the centroid
plot_community(nest_comm_large_bias,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Nested with Bias Towards Centroid")

## ----conserved----------------------------------------------------------------
# Simulating the community
cons_comm <- conserved_ellipses(object = ref_ellipse,
                                background = back_data[, vars],
                                n = 20)

# check the a few details from the generated community
print(cons_comm)  # a summary of the elements in the community object

## ----conserved_plot-----------------------------------------------------------
# Plotting the community withe the background
par(mar = mars)  # adjust margins for better visualization

plot_community(cons_comm, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Community of Conserved Ellipses")

## ----cons_den-----------------------------------------------------------------
# Simulating the community with the full background
cons_comm_full <- conserved_ellipses(object = ref_ellipse,
                                     background = back_data[, vars], n = 20)

# Simulating the community with a thinned background
cons_comm_thin <- conserved_ellipses(object = ref_ellipse,
                                     background = back_data[, vars], n = 20,
                                     thin_background = TRUE, resolution = 10)

## ----cons_den_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"--------
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Plotting the community with the full background
plot_community(cons_comm_full, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community from Full Background")

## Plotting the community with the thinned background
plot_community(cons_comm_thin, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community from Thinned Background")

## ----cons_prop----------------------------------------------------------------
# Community with both arguments small
cons_comm1 <- conserved_ellipses(object = ref_ellipse,
                                 background = back_data[, vars], n = 20,
                                 thin_background = TRUE, resolution = 10,
                                 smallest_proportion = 0.1,
                                 largest_proportion = 0.5)

# Community with smallest_proportion small and largest_proportion large
cons_comm2 <- conserved_ellipses(object = ref_ellipse,
                                 background = back_data[, vars], n = 20,
                                 thin_background = TRUE, resolution = 10,
                                 smallest_proportion = 0.1,
                                 largest_proportion = 1.5)

## ----cons_prop_plot, fig.width = 7, fig.height = 3.5, out.width = "95%"-------
# Plotting the communities
par(mfrow = c(1, 2), cex = 0.6, mar = mars)  # set up the plotting area

## Community with small proportions
plot_community(cons_comm1, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community with Small Proportions")

## Community with small smallest_proportion and large largest_proportion
plot_community(cons_comm2, background = back_data[, vars],
               pch = ".", col_bg = "#9a9797",
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Conserved Community with Small and Large")

## ----predict------------------------------------------------------------------
# Predicting Mahalanobis distances
maha_cons_pred <- predict(cons_comm, newdata = back_data[, vars],
                          prediction = "Mahalanobis")

# Predicting suitability
suit_cons_pred <- predict(cons_comm, newdata = back_data[, vars],
                          prediction = "suitability")

# Check Mahalanobis predictions
maha_cons_pred[1:5, 1:5]

# Check suiatbility predictions
suit_cons_pred[1:5, 1:5]

## ----predict_plot-------------------------------------------------------------
# Plotting the results
## Plotting area parameters
par(mfrow = c(2, 2), cex = 0.6, mar = mars)  # adjust margins for visualization

## plot mahalanobis distance predictions for ellipse 1
plot_ellipsoid(object = cons_comm[[3]][[1]],  # the relevant ellipse
               prediction = maha_cons_pred,  # mahalanobis distance prediction
               col_layer = "ell_1",  # ellipse prediction to use for colors
               pal = hcl.colors(100, palette = "Oslo"),  # color palette
               rev_pal = TRUE,  # reversing the palette
               col_ell = "#e10000", lwd = 2,  # color and line for ellipse
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Mahalanobis Distance Ellipse 1")

# plot mahalanobis distance predictions for ellipse 2
plot_ellipsoid(object = cons_comm[[3]][[2]],
               prediction = maha_cons_pred,
               col_layer = "ell_2",
               pal = hcl.colors(100, palette = "Oslo"),
               rev_pal = TRUE,
               col_ell = "#e10000", lwd = 2,  # color and line for ellipse
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Mahalanobis Distance Ellipse 2")

# plot suitability predictions for ellipse 1
plot_ellipsoid(object = cons_comm[[3]][[1]],
               prediction = suit_cons_pred,  # suitability prediction
               col_layer = "ell_1",
               pal = hcl.colors(100, palette = "Viridis"),
               col_ell = "#e10000", lwd = 2,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Suitability Ellipse 1")

# plot suitability predictions for ellipse 2
plot_ellipsoid(object = cons_comm[[3]][[2]],
               prediction = suit_cons_pred,  # suitability prediction
               col_layer = "ell_2",
               pal = hcl.colors(100, palette = "Viridis"),
               col_ell = "#e10000", lwd = 2,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Suitability Ellipse 2")

## ----predictr-----------------------------------------------------------------
# Predicting Mahalanobis distances
maha_cons_predr <- predict(cons_comm, newdata = ma_bios,
                           prediction = "Mahalanobis")

# Predicting suitability
suit_cons_predr <- predict(cons_comm, newdata = ma_bios,
                           prediction = "suitability")

# Check Mahalanobis predictions
maha_cons_predr

# Check suiatbility predictions
suit_cons_predr

## ----predictr_plot, fig.width = 7, fig.height = 4.5, out.width = "95%"--------
# Plotting area parameters
par(mfrow = c(2, 2), cex = 0.6)  # adjust margins for visualization
marsr <- c(0.5, 0.5, 2, 4)

# Plots
terra::plot(maha_cons_predr$ell_1,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Mahalanobis Distance Ellipse 1")

terra::plot(maha_cons_predr$ell_2,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Mahalanobis Distance Ellipse 2")

terra::plot(suit_cons_predr$ell_1,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Ellipse 1")

terra::plot(suit_cons_predr$ell_2,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Ellipse 2")

## ----trunc--------------------------------------------------------------------
# Predicting suitability truncated using a data frame
suit_cons_predt <- predict(cons_comm, newdata = back_data[, vars],
                           prediction = "suitability_trunc")

# Predicting suitability truncated using raster data
suit_cons_predrt <- predict(cons_comm, newdata = ma_bios,
                            prediction = "suitability_trunc")

# Check predictions in data.frame
suit_cons_predt[1:5, 1:5]

# Check predictions in raster
suit_cons_predrt

## ----trunc_plot, fig.width = 7, fig.height = 3, out.width = "95%"-------------
# Plotting area parameters
par(mfrow = c(1, 2), cex = 0.6, mar = mars)

## the truncated results for the background
plot_ellipsoid(object = cons_comm[[3]][[1]],
               prediction = suit_cons_predt,  # suitability prediction
               col_layer = "ell_1",
               pal = hcl.colors(100, palette = "Viridis"),
               col_ell = "NA", lwd = 0,
               xlab = "Bio1 (Mean Annual Temperature)",
               ylab = "Bio12 (Annual Precipitation)",
               main = "Suitability Trunc. E Space")

## custom color palette for raster predictions
tr_cols <- c("gray", hcl.colors(100, palette = "Viridis"))
tr_breacks <- c(0, seq(1e-6, 1, length.out = 101))

## the truncated raster
terra::plot(suit_cons_predrt$ell_1,
            col = tr_cols, breaks = tr_breacks,
            type = "continuous",
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Trunc. G Space")

## ----trunc_plot2, fig.width = 7, fig.height = 3, out.width = "95%"------------
# Obtaining values of zero and one
## results in data.frame
bin_suit_cons_predt <- suit_cons_predt
bin_suit_cons_predt[, -(1:2)] <- (bin_suit_cons_predt[, -(1:2)] > 0) * 1

## results in raster
bin_suit_cons_predrt <- suit_cons_predrt
bin_suit_cons_predrt <- (bin_suit_cons_predrt > 0) * 1


# Plotting
## Colors for suitability
bincol <- c("#c9c9c9", "#0004d5")

## Plotting area parameters
par(mfrow = c(1, 2), cex = 0.6, mar = mars)

## the binary for the background
plot(bin_suit_cons_predt[, vars],  # plot the variables as points
     col = bincol[as.factor(bin_suit_cons_predt$ell_1)],
     pch = 16, xlab = "Bio1 (Mean Annual Temperature)",
     ylab = "Bio12 (Annual Precipitation)",
     main = "Suitability Binary E Space")

## the binary raster
terra::plot(bin_suit_cons_predrt$ell_1, col = bincol,
            axes = FALSE, box = TRUE, mar = marsr,
            main = "Suitability Binary G Space")

## ----pam_plot-----------------------------------------------------------------
# Exclude sites with no species to make the plot easier
pam <- bin_suit_cons_predt[, -(1:2)]
pam <- as.matrix(pam[!rowSums(pam) == 0, ])

# Plot PAM using image
par(mar = c(1, 1, 5, 1), cex = 0.6)

image(1:ncol(pam), 1:nrow(pam), t(pam[nrow(pam):1, ]),
      col = bincol, axes = FALSE)
box()
text(x = 1:ncol(pam), y = nrow(pam) * 1.01,
     labels = colnames(pam),
     srt = 90, adj = 0, xpd = TRUE)

title(main = "Species Presence-Absence Matrix", line = 3.5)

## ----richness, fig.width = 7, fig.height = 4.5, out.width = "95%"-------------
# Compute richness
richness <- terra::app(bin_suit_cons_predrt, sum)

# Plot richness
terra::plot(richness, mar = marsr, col = rev(heat.colors(20)),
            main = "Species Richness Map")

## ----par_reset----------------------------------------------------------------
# Reset plotting parameters
par(original_par)

## ----save_import, eval=FALSE--------------------------------------------------
# # file name (in a temporary directory for demonstration purposes)
# temp_file <- file.path(tempdir(), "conserved_community.rds")
# 
# # Save the community object to a local directory
# save_nicheR(cons_comm, file = temp_file)
# 
# # Import the community object from a local directory
# read_com <- read_nicheR(temp_file)

## ----save_import2, eval=FALSE-------------------------------------------------
# # file names (in a temporary directory for demonstration purposes)
# temp_df_file <- file.path(tempdir(), "df_cons_com_predictions.csv")
# temp_raster <- file.path(tempdir(), "raster_cons_com_predictions.tif")
# 
# # Save predictions in data.frame objects
# write.csv(suit_cons_predt, file = temp_df_file, row.names = FALSE)
# 
# # Save predictions in raster objects
# terra::writeRaster(suit_cons_predt, filename = temp_raster)
# 
# # Import predictions as data.frame objects
# read_com_pred_df <- read.csv(temp_df_file)
# 
# # Import predictions as SpatRaster objects
# read_com_pred_ras <- terra::rast(temp_raster)