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

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

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

# 1. Load reference niche (nicheR_ellipsoid object)
data("ref_ellipse", package = "nicheR")

# 2. Load pre-calculated prediction surfaces (from previous vignettes)
# These SpatRasters contain "suitability", "Mahalanobis", "suitability_trunc", etc.
pred <- terra::rast(system.file("extdata", "predictions_rast.tif", package = "nicheR"))
pred_3d <- terra::rast(system.file("extdata", "predictions_3d_rast.tif", package = "nicheR"))

## ----results='hide', message=FALSE, warning=FALSE-----------------------------
# Load a sample bias layer containing 'sp_richness' and 'nighttime'
biases_file <- system.file("extdata", "ma_biases.tif", package = "nicheR")
raw_bias <- terra::rast(biases_file)

# --- Plotting the Output ---
par(mfrow = c(1, 2), mar = c(4, 4, 3, 2))

# Plot the raw inputs
terra::plot(raw_bias[["sp_richness"]], main = "Species Richness")
terra::plot(raw_bias[["nighttime"]], main = "Nighttime Lights")

## ----results='hide', message=FALSE, warning=FALSE-----------------------------
# Prepare a composite bias surface mapping unique directions to each layer
prep_composite <- prepare_bias(bias_surface = raw_bias, 
                               effect_direction = c("direct", "inverse"), 
                               verbose = FALSE)
# Plot the resulting unified bias probability surface
terra::plot(prep_composite$composite_surface, main = "Composite Bias Surface")

## -----------------------------------------------------------------------------
# Apply the composite bias to our suitability layer
applied_bias <- apply_bias(prepared_bias = prep_composite, 
                           prediction = pred, 
                           prediction_layer = "suitability",
                           effect_direction = "direct")

# --- Plotting the Output ---
par(mfrow = c(1, 2), mar = c(4, 4, 3, 2))

# Original Biological Suitability
terra::plot(pred[["suitability"]], main = "Habitat Suitability")

# Suitability mathematically restricted by our composite sampling bias
terra::plot(applied_bias[[1]], main = "Suitability + Composite Bias")

## -----------------------------------------------------------------------------
# Apply the composite bias to our 3D suitability layer
applied_bias_3d <- apply_bias(prepared_bias = prep_composite, 
                              prediction = pred_3d, 
                              prediction_layer = "suitability",
                              effect_direction = "direct")

# --- Plotting the Output ---
par(mfrow = c(1, 2), mar = c(4, 4, 3, 2))

terra::plot(pred_3d[["suitability"]], main = "3D Habitat Suitability")
terra::plot(applied_bias_3d[[1]], main = "3D Suitability + Composite Bias")

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

## -----------------------------------------------------------------------------
# Save the final biased prediction layers to a temporary directory
temp_rast <- file.path(tempdir(), "applied_bias_rast.tif")
temp_rast_3d <- file.path(tempdir(), "applied_bias_3d_rast.tif")

terra::writeRaster(applied_bias[[1]], filename = temp_rast, overwrite = TRUE)
terra::writeRaster(applied_bias_3d[[1]], filename = temp_rast_3d, overwrite = TRUE)