## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", # fig.path = "man/figures/README-", # out.width = "100%", fig.width = 10, fig.height = 6, # dpi = 200, echo = TRUE, warning = FALSE, eval = TRUE, message = FALSE, warning = FALSE ) ## ----dependencies, include = FALSE-------------------------------------------- # devtools::install_github('sjevelazco/flexsdm') library(knitr) ## ----packages, eval = FALSE--------------------------------------------------- # # Install and load required packages # if (!require(flexsdm)) { # remotes::install_github("flexsdm") # } # if (!require(terra)) { # install.packages("terra") # } # if (!require(geodata)) { # install.packages("geodata") # } # if (!require(tidyterra)) { # install.packages("tidyterra") # } # if (!require(ggspatial)) { # install.packages("ggspatial") # } # if (!require(ggplot2)) { # install.packages("ggplot2") # } # if (!require(dplyr)) { # install.packages("dplyr") # } # if (!require(tidyr)) { # install.packages("tidyr") # } # if (!require(readr)) { # install.packages("readr") # } # if (!require(sf)) { # install.packages("sf") # } # if (!require(reshape2)) { # install.packages("reshape2") # } # if (!require(stringr)) { # install.packages("stringr") # } # if (!require(purrr)) { # install.packages("purrr") # } # # sapply( # c( # "flexsdm", # "terra", # "geodata", # "tidyterra", # "ggplot2", # "dplyr", # "tidyr", # "readr", # "sf", # "ggspatial", # "reshape2", # "stringr", # "purrr" # ), # require, # character.only = TRUE # ) # # # Set seed for reproducibility # set.seed(123) ## ----directory_setup, eval = FALSE-------------------------------------------- # # Future scenarios (we will use the function geodata::cmip6_world() to download future climate data but this will be highly dependent on the source of your future data # # # Define scenario components # ssps <- c("370", "585") # two Shared Socioeconomic Pathway scenarios # time <- c("2041-2060") # just one time frame but you can have multiple # models <- c("CNRM-CM6-1", "MIROC6") # two climate models # # # Expand into all combinations # scenarios <- expand.grid( # model = models, # ssp = ssps, # time = time, # stringsAsFactors = FALSE # ) # # # Create identifiers for each scenario combination # scenarios$filename <- apply(scenarios, 1, function(x) { # paste0( # tolower(gsub("-", "", x["model"])), # "_", # x["ssp"], # "_", # gsub("-", "_", x["time"]) # ) # }) # # # # Let's create a directory to create a directory system to store models inputs and outputs # temp_dir <- tempdir() # This directory is temporary. We advise using a non-temporary directory to keep outputs saved when you test or adapt these codes. # proj_dir <- file.path(temp_dir, "1_SDM") # dir.create(proj_dir) # # proj_dir <- flexsdm::sdm_directory( # main_dir = proj_dir, # projections = scenarios$filename, # algorithm = c( # "gam", # "gau", # "glm", # "gbm", # "max", # "net", # "raf", # "svm", # "mean", # "median" # ), # threshold = FALSE, # return_vector = TRUE # ) # # # You can now use the proj_dir object as a shortcut for saving outputs # # # It will be helpful to create a folder for saving figures # # Figure directory # dir_fig <- file.path(proj_dir[1], "2_Outputs", "3_Figures") # dir.create(dir_fig) # # # Path to save data based on points (i.e. presences, pseudo-absences, and background points) # points_dir <- proj_dir[grep( # paste0("1_Inputs/1_Occurrences"), # proj_dir # )] # points_dir ## ----save_occurrences, eval = FALSE------------------------------------------- # # Load palms occurrence data # data(palms) # palms # # # Save the data to a file ./1_SDM/1_Inputs/1_Occurrence # points_dir # readr::write_tsv(palms, file.path(points_dir, "0_south_br_palms.txt")) ## ----load_data, eval = FALSE-------------------------------------------------- # # Load species data # species_data <- readr::read_tsv(file.path(points_dir, "0_south_br_palms.txt")) # # # Check the number of occurrences available for each species # species_counts <- table(species_data$species) %>% sort() # species_counts # # # Let's remove species with fewer than 15 occurrences # species_data <- species_data %>% # dplyr::group_by(species) %>% # dplyr::filter(n() >= 15) %>% # dplyr::ungroup() # # # Plot occurrence points to visualize the data # ggplot(data = species_data, aes(x, y)) + # geom_point(aes(color = species)) + # theme_minimal() + # labs(title = "Palm species occurrences in Southern Brazil") + # facet_grid(. ~ species) ## ----get_environmental_data, eval = FALSE------------------------------------- # # Define area of interest (bounding box for South Brazil) # bbox <- c(-70, -30, -36, -5) # xmin, xmax, ymin, ymax <- define this based on your study extent of all your species # # # If you have a polygon for your study extent, you can define like bbox <- st_bbox(study_extent) # # # Download WorldClim bioclimatic variables at 10 arc-minutes resolution (you can adjust the spatial resolution with the res argument) # env_data <- geodata::worldclim_global(var = "bio", res = 10, path = temp_dir) # tempdir() # # # Crop to our area of interest -- retain all variables for now -- we will correct for multicollinearity later # env_data <- terra::crop(env_data, terra::ext(bbox)) # # # Clean variables names # names(env_data) <- names(env_data) %>% # gsub("wc2.1_10m_", "", .) # # # Plot a raster layers and points # plot(env_data[[1]]) # points(species_data %>% dplyr::select(x, y)) # # # Save raster for current conditions # env_path <- proj_dir[grep( # paste0("1_Inputs/2_Predictors/1_Current"), # proj_dir # )] # terra::writeRaster( # env_data, # file.path(env_path, "bio_var.tif"), # overwrite = TRUE # ) ## ----plot_variables, eval = FALSE--------------------------------------------- # # Plot the environmental variables # # ---- Temperature Variables (1–11) ---- # plot( # env_data[[1:11]], # main = c( # "Annual Mean Temperature", # "Mean Diurnal Range", # "Isothermality", # "Temperature Seasonality", # "Max Temperature of Warmest Month", # "Min Temperature of Coldest Month", # "Temperature Annual Range", # "Mean Temperature of Wettest Quarter", # "Mean Temperature of Driest Quarter", # "Mean Temperature of Warmest Quarter", # "Mean Temperature of Coldest Quarter" # ) # ) # # # # ---- Precipitation Variables (12–19) ---- # plot( # env_data[[12:19]], # main = c( # "Annual Precipitation", # "Precipitation of Wettest Month", # "Precipitation of Driest Month", # "Precipitation Seasonality", # "Precipitation of Wettest Quarter", # "Precipitation of Driest Quarter", # "Precipitation of Warmest Quarter", # "Precipitation of Coldest Quarter" # ) # ) ## ----download_fut_data, eval = FALSE------------------------------------------ # # Download future climate change projecitons of the bioclim variables # fut_env <- list() # # for (i in 1:nrow(scenarios)) { # print(scenarios$filename[i]) # # fut <- geodata::cmip6_world( # model = scenarios$model[i], # ssp = scenarios$ssp[i], # time = scenarios$time[i], # var = "bioc", # res = 10, # path = tempdir() # ) # # fut_env[[i]] <- terra::crop(fut, terra::ext(bbox)) # if you want your future data cropped to the extent we defined earlier # # names(fut_env[[i]]) <- paste0("bio_", 1:nlyr(fut_env[[i]])) # # path <- file.path( # temp_dir, # "1_SDM/1_Inputs/2_Predictors/2_Projection", # scenarios$filename[i] # ) # # # Save raster for future scenario # terra::writeRaster( # fut_env[[i]], # file.path(path, paste0(scenarios$filename[i], "_bio_var.tif")), # overwrite = TRUE # ) # } # # names(fut_env) <- scenarios$filename # # # Compare bio1 across the four scenarios # # Get global min/max for bio_1 across all scenarios # bio1_range <- range(unlist(lapply(fut_env, function(r) minmax(r[["bio_1"]])))) # # # Set up plotting window (2x2 for 4 scenarios) # par(mfrow = c(2, 2)) # # # Plot with consistent zlim # for (i in seq_along(fut_env)) { # plot(fut_env[[i]][["bio_1"]], main = scenarios$filename[i], zlim = bio1_range) # } ## ----calib_area, eval = FALSE------------------------------------------------- # # A vector with species names that will be modeled # sp_names <- unique(species_data$species) %>% sort() # # # Directory to save training areas to shapefiles # calib_dir <- proj_dir[grep( # paste0("1_Inputs/3_Calibration_area"), # proj_dir # )] # calib_dir # # # # Create training areas for each species # for (i in 1:length(sp_names)) { # message("Creating training area for species: ", sp_names[i]) # # # Define training area for the i-th species # ca <- calib_area( # data = species_data %>% dplyr::filter(species == sp_names[i]), # x = "x", # y = "y", # method = c("buffer", width = 300000), # 300 km around presences # crs = crs(env_data) # ) # # # Save shapefile with tranning area # terra::writeVector( # ca, # filename = file.path(calib_dir, paste0(sp_names[i], ".gpkg")), # overwrite = TRUE # ) # } # # # Check the training areas in the directory # list.files(calib_dir, full.names = TRUE) ## ----multicollinearity, eval = FALSE------------------------------------------ # # Extract environmental values at species occurrence points # species_env <- flexsdm::sdm_extract( # data = species_data, # x = "x", # y = "y", # env_layer = env_data, # filter_na = TRUE # ) # # sp_names # Vector with species to be modeled # # # List of species data frames # species_data_list <- species_env %>% # dplyr::mutate(pr_ab = 1) %>% # dplyr::group_by(species) %>% # dplyr::group_split() # # names(species_data_list) <- sp_names # # head(species_data_list) ## ----multicollinearity_2, eval = FALSE---------------------------------------- # # Vector with variable names # var_names <- names(env_data) # # # First let's visualize predictor collinearity at species' occurrences for each species # for (i in 1:length(sp_names)) { # cor_data <- species_env %>% # dplyr::filter(species == sp_names[i]) %>% # dplyr::select(all_of(var_names)) %>% # cor(use = "pairwise.complete.obs", method = "pearson") %>% # reshape2::melt() # # cor_data <- cor_data[as.numeric(cor_data$Var1) > as.numeric(cor_data$Var2), ] # # ggplot(cor_data, aes(x = Var1, y = Var2, fill = value)) + # geom_tile() + # scale_fill_gradient2( # low = "blue", # mid = "white", # high = "red", # midpoint = 0 # ) + # geom_text(aes(label = sprintf("%.2f", value)), size = 4) + # theme_bw() + # labs( # title = paste("Correlation: ", sp_names[i]), # fill = "Correlation" # ) + # theme( # axis.text.x = element_text(angle = 45, hjust = 1), # axis.title = element_blank() # ) # # # save each heatmap (adjust file path as needed) # ggsave( # filename = file.path( # dir_fig, # paste0(sp_names[i], "_correlation_heatmap.png") # ), # width = 8, # height = 6 # ) # } ## ----correct collinearity-pca one species, eval = FALSE----------------------- # # (Optional) Access the help page for the function to learn more about its parameters: # # ?correct_colinvar # # # Generate a list of file paths to each species' training area. # # Each training area is stored as a .gpkg (GeoPackage) file in calib_dir. # calib_a <- list.files(calib_dir, full.names = TRUE, pattern = ".gpkg") # # # Assign species names as the names of the elements in calib_a for easy reference. # # This extracts the file name without the ".gpkg" extension and sets it as the list name. # names(calib_a) <- basename(calib_a) %>% gsub(".gpkg", "", .) # # # Print the named vector to verify the paths and their associated species names. # calib_a # # # --- Demonstrate the process for all species ---- # # # ----------------------------------------------------------- # # Example: Running PCA-based collinearity correction # # for multiple species across multiple scenarios # # ----------------------------------------------------------- # # # List training area polygons (.gpkg files) # calib_a <- list.files(calib_dir, full.names = TRUE, pattern = ".gpkg") # # # Assign species names to each element for reference # names(calib_a) <- basename(calib_a) %>% gsub(".gpkg", "", .) # # # Inspect what we have # calib_a # # # --- Demonstrate the process for a single species (the first one) --- # # # Load the vector (polygon) corresponding to the first species’ training area. # # This will be used to restrict PCA to the relevant region. # v <- terra::vect(calib_a[sp_names[1]]) # # # Plot the training area polygon to visually confirm it's correct. # plot(v) # # # Run the correct_colinvar function to reduce multicollinearity among predictors for this species. # # - env_layer: the stack of environmental predictor rasters. # # - restric_to_region: restricts the analysis to the training area polygon. # # - restric_pca_proj: ensures PCA is projected only within the training area. # # - method: use "pca" (Principal Component Analysis) for dimensionality reduction. # # - based_on_points: FALSE means PCA is based on the entire environmental layer within the region, not just the occurrence points. # pca_out <- flexsdm::correct_colinvar( # env_layer = env_data, # restric_to_region = v, # Restrict to the training area polygon # restric_pca_proj = TRUE, # Apply PCA projection only to the training area # method = c("pca"), # Use Principal Component Analysis # based_on_points = FALSE # Base PCA on environmental data, not just points # ) # # # Print the results of the PCA correction to examine the structure of the output object. # pca_out # # # Plot the resulting raster of principal components. # # This raster will be used as the new, uncorrelated set of predictors for SDM for this species. # plot(pca_out$env_layer) # # # Display the cumulative variance explained by each principal component. # # This helps determine how many components are needed to capture most of the environmental variation. # pca_out$cumulative_variance # # # Show the PCA loadings (coefficients) for each original variable. # # These indicate the contribution of each environmental variable to each principal component. # pca_out$coefficients ## ----correct collinearity -pca all species, eval = FALSE---------------------- # # Lets create a new folder to store variables of each species # save_dir <- file.path(proj_dir[2], "4_Predictors_by_sp") # dir.create(save_dir) # save_dir # # # path with future scenarios # scenario_path <- proj_dir[grep( # paste0("1_Inputs/2_Predictors/2_Projection"), # proj_dir # )][1] # # # Separate directories for current and future # folders <- c("1_Current", "2_Projection") # # # Create subfolders inside the base folder # folder_paths <- sapply( # folders, # function(f) { # dir_path <- file.path(save_dir, f) # if (!dir.exists(dir_path)) { # dir.create(dir_path, recursive = TRUE) # } # dir_path # return the created path # }, # USE.NAMES = TRUE # ) # # env_new <- list() # # # Perform PCA for each species and save the outputs # for (i in 1:length(sp_names)) { # message("Processing species: ", sp_names[i]) # # v <- calib_a[i] %>% terra::vect() # # # species-specific predictor projection folders # # fut_folder <- file.path(folder_paths[2], sp_names[i]) # dir.create(fut_folder) # # # Run PCA collinearity correction # env_new[[i]] <- flexsdm::correct_colinvar( # env_layer = env_data, # proj = scenario_path, # save_proj = fut_folder, # restric_to_region = v, # restric_pca_proj = TRUE, # method = "pca", # based_on_points = FALSE # ) # # # Save tables # # env_new[[i]]$cumulative_variance %>% # readr::write_tsv(file.path( # folder_paths[1], # paste0(sp_names[i], "-cum_var.txt") # )) # # env_new[[i]]$coefficients %>% # readr::write_tsv(file.path( # folder_paths[1], # paste0(sp_names[i], "-coefficients.txt") # )) # # # Saver raster for current conditions # terra::writeRaster( # env_new[[i]]$env_layer, # file.path(folder_paths[1], paste0(sp_names[i], ".tif")), # overwrite = TRUE # ) # } ## ----sample_bias_correction, eval = FALSE------------------------------------- # # 'species_data' holds all occurrence records for all species to be modeled. # species_data # Data frame: species names, coordinates, etc. # # # 'sp_names' contains the names of the species to be included in the analysis. # sp_names # Character vector: species to process # # # 'save_dir' is the directory where each species' PCA-reduced environmental predictor files are saved. # save_dir # Path to directory with species-specific predictors # # # 'folders' are the directories where the PCA-reduced environmental predictors are split between Current and Projections # file.path(save_dir, folders) # # # List all .tif files in folders (one per species, each containing PCA environmental layers for the baseline time). # sp_pc <- list.files( # file.path(save_dir, folders[1]), # full.names = TRUE, # pattern = ".tif" # ) # # # Assign species names to each file in 'sp_pc' by removing the '.tif' extension from the filename. # names(sp_pc) <- basename(sp_pc) %>% gsub(".tif", "", .) # Now names(sp_pc) matches sp_names # # # Print the named vector to verify that the PCA files are matched to the correct species. # sp_pc # Named vector: path to each species' PCA raster # # # Initialize empty lists to store filtered occurrences and bin selection results for each species. # filtered_occurrences <- list() # bins_results <- list() # # # Loop over each species to perform occurrence thinning in environmental space. # for (s in 1:length(sp_names)) { # message("Filtering occurrences for species: ", s) # # # Extract occurrence data for the current species. # occ_data <- species_data %>% # dplyr::filter(species == sp_names[s]) # # # Add a unique identifier to each occurrence record (required for filtering). # occ_data$idd <- seq_len(nrow(occ_data)) # # # Load the PCA-reduced environmental raster for this species. # env_data <- sp_pc[sp_names[s]] %>% terra::rast() # # # Filter occurrences in environmental space using a range of bin numbers. # # This step thins points so that only one point is retained per environmental bin, reducing sampling bias. # filtered_dif_bins <- flexsdm::occfilt_env( # data = occ_data, # x = "x", # y = "y", # id = "idd", # env_layer = env_data, # nbins = c(8, 10, 12, 14, 18, 20, 24, 26, 30) # Try different numbers of bins # ) # # # Select the optimal bin size (trade-off between spatial independence and number of retained records). # # The function returns the filtered occurrences and summary statistics for each bin size. # occ_selected <- flexsdm::occfilt_select( # occ_list = filtered_dif_bins, # x = "x", # y = "y", # env_layer = env_data, # filter_prop = TRUE # Optimize for spatial independence and number of points retained # ) # # # Save the filtered occurrences and bin selection results for this species. # filtered_occurrences[[s]] <- occ_selected$occ # bins_results[[s]] <- occ_selected$filter_prop # # # Clean up memory (especially useful in large loops). # gc() # } # # # Assign species names as list names for easier downstream handling. # names(filtered_occurrences) <- sp_names # names(bins_results) <- sp_names # # # Combine filtered occurrences from all species into a single data frame with a 'species' column. # filtered_occurrences <- dplyr::bind_rows(filtered_occurrences, .id = "species") # filtered_occurrences # # # Combine bin selection summaries for all species into a single data frame. # bins_results <- dplyr::bind_rows(bins_results, .id = "species") # bins_results # # # Save bin selection results and filtered occurrences to tab-separated files for reproducibility and downstream use. # readr::write_tsv( # bins_results, # file.path(points_dir, "0_filt_prop_selected.txt") # ) # readr::write_tsv( # filtered_occurrences, # file.path(points_dir, "0_occurrences_filtered.txt") # ) # # # Compare the number of original (unfiltered) and filtered occurrence records per species. # # This helps assess how much data was removed by the filtering process. # left_join( # species_data %>% group_by(species) %>% summarise(unfiltered = n()), # filtered_occurrences %>% group_by(species) %>% summarise(filtered = n()), # by = "species" # ) ## ----pres_pabsences, eval = FALSE--------------------------------------------- # # Load the filtered occurrence data from file. # # This table contains the bias-corrected presence records for all species. # filtered_occurrences <- readr::read_tsv(file.path( # proj_dir[2], # "1_Occurrences/0_occurrences_filtered.txt" # )) # # # Add a column indicating presences (pr_ab = 1), as all records here are presences. # filtered_occurrences$pr_ab <- 1 # # # Prepare to load environmental PCA predictors for each species. # # Find all .tif files in the '4_Predictors_by_sp' directory (one per species). # sp_pc <- file.path(proj_dir[2], "4_Predictors_by_sp/1_Current/") %>% # list.files(., full.names = TRUE, pattern = ".tif") # # # Set the names of the vector to the species names (by stripping ".tif" extension from the file names). # names(sp_pc) <- basename(sp_pc) %>% gsub(".tif$", "", .) # # # Print out the named vector to confirm the mapping of species to file paths. # sp_pc # # # Initialize empty lists to store results for all species: # # - pres_pabsences: for data with presences and pseudo-absences # # - background: for background points (for Maxent models) # pres_pabsences <- list() # background <- list() # # # Loop over each species to generate pseudo-absences and background points. # for (i in 1:length(sp_names)) { # message("species ", i) # # # Extract filtered occurrence records for the current species. # coord <- dplyr::filter(filtered_occurrences, species == sp_names[i]) # # # Load the PCA-reduced environmental raster for the current species. # pred <- terra::rast(sp_pc[sp_names[i]]) # # set.seed(15) # Set random seed for reproducibility. # # # Generate pseudo-absence points for the current species. # # n = twice the number of presences (a common SDM practice). # # Points are sampled randomly within the valid region of the predictor raster. # pres_pabsences[[i]] <- # flexsdm::sample_pseudoabs( # data = coord, # x = "x", # y = "y", # n = nrow(coord) * 2, # 2x number of presences # method = "random", # Random selection # rlayer = pred # Within the species-specific training area # ) %>% # dplyr::mutate( # species = sp_names[i] # Label with species name # ) %>% # dplyr::bind_rows(coord, .) # Combine presences and pseudo-absences # # set.seed(15) # Reset seed for reproducibility. # # # Generate 10,000 random background points for the current species (for Maxent). # background[[i]] <- flexsdm::sample_background( # data = coord, # x = "x", # y = "y", # n = 10000, # method = "random", # rlayer = pred # ) %>% # dplyr::mutate( # species = sp_names[i] # ) # } # # Assign species names to each list element for easy reference later. # names(background) <- sp_names # # # Combine the lists of data frames into single data frames for all species. # pres_pabsences <- dplyr::bind_rows(pres_pabsences) # background <- dplyr::bind_rows(background) # # # Save the combined presence + pseudo-absence data and background data to disk as .txt files. # readr::write_tsv( # pres_pabsences, # file.path(points_dir, "0_occurrences_psabsences.txt") # ) # readr::write_tsv(background, file.path(points_dir, "0_background.txt")) ## ----plot pseudo-absences, eval = FALSE--------------------------------------- # # This section generates maps visualizing the spatial distribution of: # # 1) species presence locations, # # 2) pseudo-absence points, and # # 3) background points for each species. # # The resulting maps will help to visually assess the sampling design and spatial coverage of your data. # # dir_fig # Directory where figure files will be saved # # # List all training area shapefiles (.gpkg) for each species. # # These polygons define the modeling extent (calibration area) for each species. # calib_a <- list.files(calib_dir, full.names = TRUE, pattern = ".gpkg") # # # Assign species names to each element of the calib_a vector, # # by extracting the filename and removing the ".gpkg" extension. # names(calib_a) <- basename(calib_a) %>% gsub(".gpkg", "", .) # # # Loop through species and generate maps # for (i in 1:length(sp_names)) { # # Isolate species data # pa_data <- pres_pabsences %>% dplyr::filter(species == sp_names[i]) # bg_data <- background %>% dplyr::filter(species == sp_names[i]) # ca <- terra::vect(calib_a[sp_names[i]]) # # # Plot # p <- ggplot() + # geom_spatvector( # data = ca, # fill = NA, # color = "black", # linetype = "dashed" # ) + # geom_point( # data = bg_data, # aes(x = x, y = y), # color = "gray60", # size = 0.7, # alpha = 0.5 # ) + # geom_point(data = pa_data, aes(x = x, y = y, color = as.factor(pr_ab))) + # scale_color_manual( # values = c("1" = "blue", "0" = "red"), # ) + # coord_sf() + # theme_bw() + # labs( # title = paste("Species:", sp_names), # subtitle = "Blue = Presence | Red = Pseudo-Absence | Gray = Background", # x = "Longitude", # y = "Latitude" # ) + # annotation_scale(location = "br", width_hint = 0.3) + # annotation_north_arrow( # location = "bl", # style = north_arrow_fancy_orienteering, # height = unit(1.5, "cm"), # width = unit(1.5, "cm") # ) + # labs(color = element_blank()) # # # Print the plot to the console # p # # # Save # ggsave( # filename = file.path(dir_fig, paste0(sp_names[i], "_pr-ab-bg_map.png")), # plot = p, # width = 8, # height = 6, # dpi = 300 # ) # } ## ----ddata-partition, eval = FALSE-------------------------------------------- # # Initialize lists to store partitioned data for each species. # # - partitioned_psa: will hold partitioned presence and pseudo-absence records. # # - partitioned_bg: will hold partitioned background points. # partitioned_psa <- list() # partitioned_bg <- list() # # # Loop through each species to perform spatial partitioning for cross-validation. # for (i in seq_along(sp_names)) { # message("Processing spatial partitioning for: ", i) # # # Select presence + pseudo-absence data for the current species. # psa_select <- pres_pabsences %>% dplyr::filter(species == sp_names[i]) # # # Select background points for the current species. # bg_select <- background %>% dplyr::filter(species == sp_names[i]) # # # Load PCA-reduced environmental raster for the species (used for spatial partitioning). # env_subset <- terra::rast(sp_pc[sp_names[i]]) # # # Attempt spatial partitioning using spatial blocks (grids). # # This method divides the study area into blocks, assigning data in each block to the same fold. # part <- # flexsdm::part_sblock( # env_layer = env_subset, # data = psa_select, # x = "x", # y = "y", # pr_ab = "pr_ab", # n_part = 4, # Number of partitions/folds # min_res_mult = 4, # Controls minimum block size (relative to raster resolution) # max_res_mult = 300, # Controls maximum block size # num_grids = 60, # Number of candidate grids to test # min_occ = 5, # Minimum number of occurrences per fold # prop = 0.9 # Proportion of data to be included in the partition (for exclusion of sparse blocks) # ) # # # If block partitioning fails (e.g., too few records per block), try latitudinal bands. # if (length(part) < 3) { # part <- # flexsdm::part_sband( # env_layer = env_subset, # data = psa_select, # x = "x", # y = "y", # pr_ab = "pr_ab", # type = "lat", # Split into horizontal (latitudinal) bands # n_part = 4, # min_bands = 2, # max_bands = 60, # min_occ = 5, # prop = 0.9 # ) # } # # If latitudinal bands also fail, try longitudinal bands. # if (length(part) < 3) { # part <- # flexsdm::part_sband( # env_layer = env_subset, # data = psa_select, # x = "x", # y = "y", # pr_ab = "pr_ab", # type = "lon", # Split into vertical (longitudinal) bands # n_part = 4, # min_bands = 2, # max_bands = 60, # min_occ = 5, # prop = 0.9 # ) # } # # If all spatial partitioning methods fail, revert to standard random k-fold partitioning. # if (length(part) < 3) { # # Randomly assign presences/pseudo-absences to 5 folds. # partitioned_psa[[i]] <- flexsdm::part_random( # data = psa_select, # pr_ab = "pr_ab", # method = c(method = "kfold", folds = 5) # ) # # # Randomly assign background points to 5 folds. # partitioned_bg[[i]] <- flexsdm::part_random( # data = bg_select, # pr_ab = "pr_ab", # method = c(method = "kfold", folds = 5) # ) # } # # # If spatial partitioning was successful (contains "best_part_info"), store results and save partition details. # if (any("best_part_info" == names(part))) { # # Store partitioned presence/pseudo-absence data, adding the species name. # partitioned_psa[[i]] <- part$part %>% dplyr::mutate(species = sp_names[i]) # # # Save the spatial partition grid (as a raster) for reference or plotting later. # terra::writeRaster( # part$grid, # file.path(points_dir, paste0(sp_names[i], "_best_part.tif")), # overwrite = TRUE # ) # # # Save partitioning summary (e.g., number of records per fold, SD, etc.). # readr::write_tsv( # part$best_part_info, # file.path(points_dir, paste0(sp_names[i], "_best_part.txt")) # ) # # # Assign each background point to a partition/fold according to the spatial partitioning. # partitioned_bg[[i]] <- # flexsdm::sdm_extract( # data = bg_select, # x = "x", # y = "y", # env_layer = part$grid # ) # } # # # Run garbage collection to free memory before next iteration (useful in large species datasets). # gc() # } # # # Combine all species' partitioned presence/pseudo-absence data into a single data frame. # partitioned_psa <- dplyr::bind_rows(partitioned_psa) # # # Combine all species' partitioned background data into a single data frame. # partitioned_bg <- dplyr::bind_rows(partitioned_bg) # # # Save the final partitioned datasets to disk for model training and evaluation. # readr::write_tsv( # partitioned_psa, # file.path(points_dir, "0_partitioned_occurrences_psabsences.txt") # ) # readr::write_tsv( # partitioned_bg, # file.path(points_dir, "0_partitioned_background.txt") # ) ## ----model_evaluation - output directories, eval = FALSE---------------------- # # Define the directory where model evaluation metrics and performance summaries will be saved. # dir_perf <- file.path(temp_dir, "1_SDM/2_Outputs/0_Model_performance") # # # Define the directory where model predictions for current environmental conditions will be saved. # # This folder will store the prediction rasters for each algorithm. # dir_current <- file.path(temp_dir, "1_SDM/2_Outputs/1_Current/Algorithm") # # # 'dir_fig' is assumed to be previously defined, and will be used for saving figures such as plots and maps. # dir_fig # # # List all environmental raster files (principal component layers) for each species. # # These are the PCA-reduced predictor rasters stored in the specified directory. # sp_pc <- list.files( # file.path(temp_dir, "1_SDM/1_Inputs/4_Predictors_by_sp/1_Current"), # full.names = TRUE, # pattern = ".tif" # ) # # # List all future environmental raster files (principal component layers) for each species. # # These are the PCA-reduced predictor rasters stored in the specified directory. # dir_futur_pca <- file.path(proj_dir[2], "4_Predictors_by_sp/2_Projection/") # # # Get all species folders (exclude the base folder itself) # species_dirs <- list.dirs(dir_futur_pca, recursive = FALSE) # # # Initialize nested list # fut_sp_pc_list <- list() # # for (sp_dir in species_dirs) { # sp_name <- basename(sp_dir) # # Get all scenario folders for this species # scenario_dirs <- list.dirs(sp_dir, recursive = FALSE) # # Initialize list for species # fut_sp_pc_list[[sp_name]] <- list() # for (sc_dir in scenario_dirs) { # sc_name <- basename(sc_dir) # # List all .tif rasters in this scenario folder # tif_files <- list.files(sc_dir, pattern = "\\.tif$", full.names = TRUE) # # Load all rasters as a raster stack # if (length(tif_files) > 0) { # fut_sp_pc_list[[sp_name]][[sc_name]] <- rast(tif_files) # } else { # fut_sp_pc_list[[sp_name]][[sc_name]] <- NULL # } # } # } # # # Check structure # str(fut_sp_pc_list, max.level = 3) # # # Assign species names to each raster file in 'sp_pc' by stripping the '.tif' extension from the file name. # # This makes it easy to match species names to their respective PCA raster paths. # names(sp_pc) <- basename(sp_pc) %>% gsub(".tif", "", .) # # # Load the partitioned presence and pseudo-absence dataset from file. # # Each row corresponds to an occurrence or pseudo-absence, along with partition/fold assignments for cross-validation. # psa <- readr::read_tsv(file.path( # proj_dir[1], # "1_Inputs/1_Occurrences/0_partitioned_occurrences_psabsences.txt" # )) # # # Load the partitioned background points dataset from file. # # Background points are used in algorithms like Maxent and should also have partition assignments. # bkgrnd <- readr::read_tsv(file.path( # proj_dir[1], # "1_Inputs/1_Occurrences/0_partitioned_background.txt" # )) ## ----model_evaluation - species loop, eval = FALSE---------------------------- # # Get a unique list of species names to iterate through and process each one. # sp_names <- unique(psa$species) # # # Loop through each species name in the 'sp_names' vector. # for (i in 1:length(sp_names)) { # # Get the current species name for this iteration. # s_name <- sp_names[i] # # Print a message to the console to track progress. # message(paste0("\n===== Processing species: ", s_name, " =====")) # # ## %######################################################%## # #### 1. Prepare Data for the Current Species #### # ## %######################################################%## # # # Load the environmental raster data specific to the species. # env_r <- terra::rast(sp_pc[s_name]) # # Get the names of the environmental variables (predictors). # var_names <- names(env_r) # # # Filter the presence-absence data for the current species. # sp_pa <- psa %>% dplyr::filter(species == s_name) # # Filter the background data for the current species. # sp_bg <- bkgrnd %>% dplyr::filter(species == s_name) # # # Extract environmental data for presence-absence points, removing any points with NA values. # sp_pa <- sdm_extract( # sp_pa, # x = "x", # y = "y", # env_layer = env_r, # filter_na = TRUE # ) # # Extract environmental data for background points. # sp_bg <- sdm_extract( # sp_bg, # x = "x", # y = "y", # env_layer = env_r, # filter_na = TRUE # ) # # ## %######################################################%## # #### 2. Fit Models #### # ## %######################################################%## # message(" > Fitting models...") # # # ---- Random Forest (tune_raf) ---- # # Tune hyperparameters for a Random Forest model. # m_raf <- tune_raf( # data = sp_pa, # Presence-absence data # response = "pr_ab", # The column with presence/absence (1/0) # predictors = var_names, # Names of environmental variables # partition = ".part", # Column for data partitioning (training/testing) # grid = expand.grid( # # Define the grid of hyperparameters to test # mtry = seq(1, length(var_names), 1), # Number of variables to sample at each split # ntree = c(400, 600, 800, 1000) # Number of trees to grow # ), # thr = "max_sorensen", # Threshold to convert continuous prediction to binary # metric = "SORENSEN", # Metric to evaluate the model performance # n_cores = 4 # Number of CPU cores for parallel processing # ) # # # ---- Maxent (tune_max) ---- # # Tune hyperparameters for a Maxent model. # m_max <- tune_max( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # background = sp_bg, # Background data is required for Maxent # partition = ".part", # grid = expand.grid( # # Define the grid of hyperparameters to test # regmult = seq(0.1, 2, 0.5), # Regularization multiplier to prevent overfitting # classes = c("lq", "lqh", "lqhp", "lqhpt") # Feature classes (l=linear, q=quadratic, etc.) # ), # thr = "max_sorensen", # metric = "SORENSEN", # n_cores = 4 # ) # # # ---- Boosted Regression Tree (tune_gbm) ---- # # Tune hyperparameters for a Boosted Regression Tree model. # m_gbm <- tune_gbm( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # partition = ".part", # grid = expand.grid( # # Define the grid of hyperparameters to test # n.trees = seq(150, 200, 10), # Number of trees # shrinkage = seq(1, 1.5, 0.2), # Learning rate # n.minobsinnode = seq(1, 15, 2) # Minimum number of observations in a terminal node # ), # thr = "max_sorensen", # metric = "SORENSEN", # n_cores = 4 # ) # # # ---- Generalized Additive Model (fit_gam) ---- # # Determine the number of training samples. # n_t <- flexsdm:::n_training(data = sp_pa, partition = ".part") # # Set a starting number for basis functions (k). # candidate_k <- 20 # # Reduce 'k' until it's smaller than the number of training samples to avoid errors. # while ( # any( # n_t < # flexsdm:::n_coefficients( # data = sp_pa, # predictors = var_names, # k = candidate_k # ) # ) # ) { # candidate_k <- candidate_k - 3 # } # # Fit a GAM model with the adjusted 'k'. # m_gam <- fit_gam( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # partition = ".part", # thr = "max_sorensen", # k = candidate_k # Basis functions parameter # ) # # # ---- Generalized Linear Model (fit_glm) ---- # # Fit a GLM with a second-degree polynomial term. # m_glm <- fit_glm( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # partition = ".part", # thr = "max_sorensen", # poly = 2 # Degree of polynomial functions # ) # # # ---- Gaussian Process (fit_gau) ---- # # Fit a Gaussian Process model. # m_gau <- fit_gau( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # partition = ".part", # thr = "max_sorensen" # ) # # # ---- Neural Network (tune_net) ---- # # Tune hyperparameters for a Neural Network model. # m_net <- tune_net( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # partition = ".part", # grid = expand.grid( # # Define the grid of hyperparameters to test # size = 2:length(var_names), # Number of units in the hidden layer # decay = c(seq(0.01, 1, 0.05), 1, 3, 4, 5, 6) # Weight decay for regularization # ), # thr = "max_sorensen", # metric = "SORENSEN", # n_cores = 4 # ) # # # ---- Support Vector Machine (tune_svm) ---- # # Tune hyperparameters for a Support Vector Machine model. # m_svm <- tune_svm( # data = sp_pa, # response = "pr_ab", # predictors = var_names, # partition = ".part", # grid = expand.grid( # C = seq(2, 60, 5), # Cost parameter for regularization # sigma = c(seq(0.001, 0.3, 0.01)) # Kernel parameter # ), # thr = "max_sorensen", # metric = "SORENSEN", # n_cores = 4 # ) # # # Gather all fitted model objects (those starting with "m_") into a single named list. # model_obj_list <- mget(grep("^m_", ls(), value = TRUE)) # # Clean up the names in the list (e.g., "m_raf" becomes "raf"). # names(model_obj_list) <- gsub("m_", "", names(model_obj_list)) # # ## %######################################################%## # #### 3. Partial Dependence & Variable Importance #### # ## %######################################################%## # message(" > Calculating PDP and Variable Importance...") # # # Initialize an empty list to store variable importance results for each model. # varimp <- list() # # # Loop through each fitted model to generate plots and calculate importance. # for (i in 1:length(model_obj_list)) { # mod_name <- names(model_obj_list)[i] # model_obj <- model_obj_list[[i]] # # # Generate Partial Dependence Plots (PDP) to visualize variable effects. # # The 'p_pdp' function extracts the core model (e.g., randomForest) from the flexsdm object. # p <- p_pdp(model = model_obj$model, training_data = sp_pa) + # ggtitle(paste(s_name, "-", mod_name)) # Add a title to the plot # # # If the plot was created successfully, save it as a PNG file. # if (!is.null(p)) { # ggsave( # plot = p, # filename = file.path( # dir_fig, # paste0(s_name, "_", mod_name, "_pdp.png") # ), # width = 8, # height = 6, # dpi = 300 # ) # } # # # Calculate variable importance using a permutation-based method. # varimp[[i]] <- flexsdm::sdm_varimp( # model = model_obj, # data = sp_pa, # response = "pr_ab", # predictors = var_names, # n_sim = 30, # Number of permutations to run # n_cores = 2, # thr = "max_sorensen", # clamp = TRUE, # pred_type = "cloglog" # Use cloglog output for predictions # ) # } # # # Combine the variable importance results from all models into a single data frame. # varimp <- dplyr::bind_rows(varimp) # # Save the combined variable importance table to a text file. # readr::write_tsv(varimp, file.path(dir_perf, paste0(s_name, "_varimp.txt"))) # # # Create a bar plot to visualize variable importance across models. # varimp2 <- varimp %>% # tidyr::pivot_longer( # cols = TPR:IMAE, # names_to = "metric", # values_to = "value" # ) %>% # filter(metric %in% c("AUC", "TSS", "SORENSEN")) # Filter for specific metrics # # p_vi <- ggplot(varimp2, aes(x = predictors, y = value, fill = model)) + # scale_fill_brewer(palette = "Set1") + # geom_bar(stat = "identity", position = position_dodge()) + # facet_grid(model ~ metric, scales = "free") + # coord_flip() + # labs( # x = "Predictor variables", # y = "Variable importance" # ) + # theme_bw() # # print(p_vi) # ggsave( # filename = file.path(dir_fig, paste0(s_name, "_varImp.png")), # plot = p_vi, # width = 12, # height = 14, # dpi = 300, # unit = "cm", # scale = 1.5 # ) # # # Summarize performance metrics (e.g., SORENSEN, TSS) for all individual models. # model_perf <- flexsdm::sdm_summarize(model_obj_list) # # Save the performance summary to a text file. # readr::write_tsv( # model_perf, # file.path(dir_perf, paste0(s_name, "_models_performance.txt")) # ) # # ## %######################################################%## # #### 4. Filter Models for Ensemble #### # ## %######################################################%## # message(" > Filtering models to perform ensemble...") # # # Filter models based on performance criteria to select the "best" ones for the ensemble. # # Here, we select models with a mean Sorensen score >= 0.7. # # We also exclude models with thresholds of 0 or 1, as they might indicate overfitting. # good_models_names <- model_perf$model[ # model_perf$SORENSEN_mean >= 0.7 & # model_perf$thr_value > 0 & # model_perf$thr_value < 1 # ] # # # Create a new list containing only the selected "good" models. # good_models <- model_obj_list[good_models_names] # # ## %######################################################%## # #### 5. Ensemble Modeling #### # ## %######################################################%## # # Initialize an empty list to store ensemble models. # ensemble_list <- list() # # # Only proceed if there is more than one "good" model to create an ensemble. # if (length(good_models) > 1) { # message(" > Fitting model ensemble...") # # Create an ensemble by averaging the predictions of the good models. # m_mean <- flexsdm::fit_ensemble( # good_models, # ens_method = "mean", # Ensemble method: mean, weighted_mean, median, etc. # thr_model = "max_sorensen", # Thresholding method for individual models # thr = "max_sorensen", # Thresholding method for the final ensemble # metric = "SORENSEN" # ) # # # Create an ensemble using the median of the predictions. # m_median <- # flexsdm::fit_ensemble( # good_models, # ens_method = "median", # thr_model = "max_sorensen", # thr = "max_sorensen", # metric = "SORENSEN" # ) # # # Combine the ensemble models into a list. # ensemble_list <- list(mean = m_mean, median = m_median) # # Summarize the performance of the ensemble models. # ens_perf <- flexsdm::sdm_summarize(ensemble_list) # # # Combine the performance tables of individual models and ensemble models. # model_perf_with_ens <- dplyr::bind_rows(model_perf, ens_perf) # # Overwrite the previous performance file to include the ensemble results. # readr::write_tsv( # model_perf_with_ens, # file.path(dir_perf, paste0(s_name, "_models_performance.txt")) # ) # } # # ## %######################################################%## # #### 6. Spatial Predictions #### # ## %######################################################%## # message(" > Predicting models ...") # # names(good_models) <- NULL # # # Predict habitat suitability for each of the "good" individual models. # if (length(good_models) > 0) { # prd_list <- flexsdm::sdm_predict( # models = good_models, # List of good models # pred = env_r, # Environmental raster layers to predict onto # thr = "max_sorensen", # clamp = TRUE, # Restrict predictions to the range of training data # con_thr = TRUE, # Produce continuous and binary (thresholded) predictions # pred_type = "cloglog" # ) # # # Loop through the list of predictions and save each one as a GeoTIFF raster file. # for (alg in names(prd_list)) { # terra::writeRaster( # prd_list[[alg]], # filename = file.path(dir_current, alg, paste0(s_name, ".tif")), # overwrite = TRUE, # filetype = "GTiff" # ) # } # } # # # Predict habitat suitability for the ensemble models. # # if (length(ensemble_list) > 0) { # for (i in seq_along(ensemble_list)) { # ens_prd <- sdm_predict( # models = ensemble_list[[i]], # pred = env_r, # thr = "max_sorensen", # con_thr = TRUE, # clamp = TRUE, # pred_type = "cloglog" # ) # # # Save prediction as a GeoTIFF. # writeRaster( # ens_prd[[1]], # filename = file.path( # dir_current, # names(ensemble_list[i]), # paste0(s_name, ".tif") # ), # overwrite = TRUE, # filetype = "GTiff" # ) # } # } # # ## %######################################################%## # #### 7. Future projections #### # ## %######################################################%## # # message(" > Predicting models for future projections ...") # # fut_env <- fut_sp_pc_list[[s_name]] # # # Loop over each scenario raster stack (fut_sp_pc_list already has species->scenario->raster) # for (sc_name in names(fut_env)) { # env_r <- fut_env[[sc_name]] # raster stack for this scenario # # # Create output folder if it doesn't exist # dir_future <- proj_dir[grep(paste0("2_Outputs/2_Projection"), proj_dir)][1] # # # --- Predict individual models --- # if (length(good_models) > 0) { # prd_list <- flexsdm::sdm_predict( # models = good_models, # pred = env_r, # thr = "max_sorensen", # clamp = TRUE, # con_thr = TRUE, # pred_type = "cloglog" # ) # # # Save individual model predictions # for (alg in names(prd_list)) { # terra::writeRaster( # prd_list[[alg]], # file.path( # dir_future, # sc_name, # paste0("Algorithm"), # alg, # paste0(s_name, ".tif") # ), # overwrite = TRUE, # filetype = "GTiff" # ) # } # } # # # --- Predict ensemble models --- # if (length(ensemble_list) > 0) { # for (i in seq_along(ensemble_list)) { # ens_prd <- flexsdm::sdm_predict( # models = ensemble_list[[i]], # pred = env_r, # thr = "max_sorensen", # clamp = TRUE, # con_thr = TRUE, # pred_type = "cloglog" # ) # # # Save ensemble prediction # writeRaster( # ens_prd[[1]], # filename = file.path( # dir_future, # sc_name, # paste0("Algorithm"), # paste0(names(ensemble_list[i])), # paste0(s_name, ".tif") # ), # overwrite = TRUE, # filetype = "GTiff" # ) # } # } # } # # ## %######################################################%## # #### 8. Clean Up for Next Iteration #### # ## %######################################################%## # # Remove all model objects (starting with "m_") to free up memory. # rm(list = grep("^m_", ls(), value = TRUE)) # # Remove another temporary object if it exists. # if (exists("alg_prd_cont")) { # rm(alg_prd_cont) # } # # Trigger garbage collection to release memory explicitly. # gc() # } # End of the main species loop. ## ----visualization_spatial_pred, eval = FALSE--------------------------------- # for (s_name in sp_names) { # message(paste(" - Processing maps for:", s_name)) # # # 1. Search recursively for all prediction raster files (.tif) for the current species. # # Each file represents a prediction from a different model or ensemble. # # The pattern ensures only files ending with the species name and .tif extension are matched. # species_files <- list.files( # path = dir_current, # pattern = paste0(s_name, ".tif$"), # '$' anchors the match to the end of the filename # recursive = TRUE, # full.names = TRUE # ) # # # If no prediction files are found for this species, output a message and skip to the next species. # if (length(species_files) == 0) { # message(paste( # " - No prediction files were found for", # s_name, # ". Skipping." # )) # next # } # # # 2. Load all prediction raster files and convert them to a format suitable for plotting with ggplot2. # r_cont <- list() # list to store continuous models # r_semibin <- list() # list to store semi-binary models # # for (i in 1:length(species_files)) { # r_cont[[i]] <- terra::rast(species_files[i])[[1]] # r_semibin[[i]] <- terra::rast(species_files[i])[[2]] # } # r_cont <- do.call("c", r_cont) # r_semibin <- do.call("c", r_semibin) # names(r_semibin) <- names(r_cont) # # # 3. Plot the continuous habitat suitability maps for this species, with one facet per model. # p_map_cont <- ggplot() + # tidyterra::geom_spatraster(data = r_cont) + # facet_wrap(~lyr) + # scale_fill_princess_c(palette = "america") + # theme_bw() # # # 4. Plot the semibinary habitat suitability maps for this species, with one facet per model. # p_map_semibin <- ggplot() + # tidyterra::geom_spatraster(data = r_semibin) + # facet_wrap(~lyr) + # scale_fill_princess_c(palette = "america") + # theme_bw() # # # 5. Save the generated map as a PNG file in the figures directory, one file per species. # ggsave( # filename = file.path(dir_fig, paste0(s_name, "_cont_map", ".png")), # plot = p_map_cont, # width = 12, # height = 10, # dpi = 300 # ) # ggsave( # filename = file.path(dir_fig, paste0(s_name, "_semibin_map", ".png")), # plot = p_map_semibin, # width = 12, # height = 10, # dpi = 300 # ) # message(paste(" - Map saved for:", s_name)) # } # # message("\nCheck figures output in 'dir_fig' path") # message("\nScript successfully completed!") ## ----vis_future_proj, eval = FALSE-------------------------------------------- # # Base folder with future projections # dir_future <- proj_dir[grep(paste0("2_Outputs/2_Projection"), proj_dir)][1] # # # List of scenarios (folders) # scenario_folders <- list.dirs(dir_future, recursive = FALSE, full.names = TRUE) # # for (s_name in sp_names) { # message("Processing species: ", s_name) # # for (sc_folder in scenario_folders) { # sc_name <- basename(sc_folder) # # # List algorithm folders # alg_folders <- list.dirs( # file.path(sc_folder, "Algorithm"), # recursive = FALSE, # full.names = TRUE # ) # # r_cont_list <- list() # r_semibin_list <- list() # alg_names <- character() # # for (alg_folder in alg_folders) { # alg_name <- basename(alg_folder) # r_file <- list.files( # alg_folder, # pattern = paste0(s_name, ".tif$"), # full.names = TRUE # ) # if (length(r_file) == 0) { # next # } # # r_stack <- terra::rast(r_file) # r_cont_list[[alg_name]] <- r_stack[[1]] # r_semibin_list[[alg_name]] <- r_stack[[2]] # alg_names <- c(alg_names, alg_name) # } # # # Skip if no raster files found # if (length(r_cont_list) == 0) { # message( # " - No raster files found for species ", # s_name, # " in scenario ", # sc_name, # ". Skipping." # ) # next # } # # # Combine rasters into a single SpatRaster for ggplot faceting # r_cont <- r_cont_list[[1]] # for (i in 2:length(r_cont_list)) { # r_cont <- c(r_cont, r_cont_list[[i]]) # } # # r_semibin <- r_semibin_list[[1]] # for (i in 2:length(r_semibin_list)) { # r_semibin <- c(r_semibin, r_semibin_list[[i]]) # } # # names(r_semibin) <- names(r_cont) # # # Continuous map # p_cont <- ggplot() + # tidyterra::geom_spatraster(data = r_cont) + # facet_wrap(~lyr) + # scale_fill_princess_c(palette = "america") + # theme_bw() + # ggtitle(paste(s_name, sc_name, "Continuous")) # # ggsave( # file.path(dir_fig, paste0(s_name, "_", sc_name, "_cont.png")), # plot = p_cont, # width = 12, # height = 10, # dpi = 300 # ) # # # Semibinary map # p_semibin <- ggplot() + # tidyterra::geom_spatraster(data = r_semibin) + # facet_wrap(~lyr) + # scale_fill_princess_c(palette = "america") + # theme_bw() + # ggtitle(paste(s_name, sc_name, "Semibinary")) # # ggsave( # file.path(dir_fig, paste0(s_name, "_", sc_name, "_semibin.png")), # plot = p_semibin, # width = 12, # height = 10, # dpi = 300 # ) # # message(" - Faceted maps saved for scenario ", sc_name) # } # }