From 18f41007dfb7334e896c0d55ebc7c38a63c77b16 Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Mon, 11 Nov 2024 22:51:41 -0800
Subject: [PATCH 1/8] Refactor: Extract reusable functions into /R

---
 R/exposure.R             | 34 +++++++++++++++++++++
 R/exposure_times.R       | 65 ++++++++++++++++++++++++++++++++++++++++
 R/extract_climate_data.R | 16 ++++++++++
 R/get_niche_limits.R     | 63 ++++++++++++++++++++++++++++++++++++++
 R/prepare_range.R        | 37 +++++++++++++++++++++++
 5 files changed, 215 insertions(+)
 create mode 100644 R/exposure.R
 create mode 100644 R/exposure_times.R
 create mode 100644 R/extract_climate_data.R
 create mode 100644 R/get_niche_limits.R
 create mode 100644 R/prepare_range.R

diff --git a/R/exposure.R b/R/exposure.R
new file mode 100644
index 0000000..a362cea
--- /dev/null
+++ b/R/exposure.R
@@ -0,0 +1,34 @@
+#' Calculate exposure for each species and grid cell
+#' @param data Index of the species
+#' @param species_range List of grid cell IDs for each species
+#' @param climate_data Climate data tibble
+#' @param niche Niche limits tibble
+#' @return A data frame with exposure results
+exposure <- function(data, species_range, climate_data, niche) {
+  # Get data for the current species
+  spp_data <- species_range[[data]]
+  spp_name <- names(species_range)[[data]]
+
+  # Filter climate data for the species' grid cells
+  spp_matrix <- climate_data %>%
+    filter(world_id %in% spp_data) %>%
+    na.omit()
+
+  # Extract niche limits for the species
+  spp_niche <- niche %>%
+    filter(species %in% spp_name)
+
+  # Compute exposure (1 if suitable, 0 if unsuitable)
+  spp_matrix <- spp_matrix %>%
+    mutate(across(2:ncol(spp_matrix), ~ case_when(
+      . <= spp_niche$niche_max ~ 1,
+      . > spp_niche$niche_max ~ 0
+    )))
+
+  # Add species column and rearrange
+  spp_matrix$species <- spp_name
+  spp_matrix <- spp_matrix %>%
+    relocate(species)
+
+  return(spp_matrix)
+}
diff --git a/R/exposure_times.R b/R/exposure_times.R
new file mode 100644
index 0000000..d07899a
--- /dev/null
+++ b/R/exposure_times.R
@@ -0,0 +1,65 @@
+#' Calculate exposure and de-exposure times for species
+#' @param data A row of exposure data (species, world_id, and years)
+#' @param original.state Initial exposure state (0 or 1)
+#' @param consecutive.elements Number of consecutive years for state change
+#' @return A tibble with exposure times and durations for the species
+exposure_times <- function(data, original.state, consecutive.elements) {
+    # Extract species and world_id
+    species <- data[1]
+    world_id <- data[2]
+
+    # Extract year data as numeric vector
+    n <- as.numeric(data[-c(1, 2)])
+
+    # Calculate shift sequences
+    rle_x <- data.frame(unclass(rle(n)))
+
+    # Add year column to represent time steps
+    rle_x$year <- 2015 + cumsum(rle_x$lengths) - rle_x$lengths
+
+    # Filter sequences with sufficient consecutive elements
+    rle_x <- rle_x[rle_x$lengths >= consecutive.elements,]
+
+    # Add a line for the original state to ensure valid transitions
+    rle_x <- rbind(c(1, original.state, 2000), rle_x)
+
+    # Remove unnecessary state repetitions
+    rle_x <- rle_x[c(TRUE, diff(rle_x$values) != 0),]
+
+    # Remove the first line (original state or duplicate)
+    rle_x <- rle_x[-1,]
+
+    # Handle cases with no valid exposure sequences
+    if (nrow(rle_x) == 0) {
+        return(tibble(species, world_id, exposure = NA, deexposure = NA, duration = NA))
+    }
+
+    # Handle cases where all values are 0 (exposure with no de-exposure)
+    if (length(unique(rle_x$values)) == 1 && unique(rle_x$values) == 0) {
+        exposure <- rle_x$year[1]
+        deexposure <- 2101  # Indicates de-exposure did not occur
+        duration <- deexposure - exposure
+        return(tibble(species, world_id, exposure, deexposure, duration))
+    }
+
+    # Handle cases with both exposure (0) and de-exposure (1)
+    if (length(unique(rle_x$values)) == 2) {
+        exposure <- rle_x %>%
+            filter(values == 0) %>%
+            pull(year)
+
+        deexposure <- rle_x %>%
+            filter(values == 1) %>%
+            pull(year)
+
+        # If there are more exposures than deexposures, add a placeholder for deexposure
+        if (length(exposure) > length(deexposure)) {
+            deexposure[length(exposure)] <- 2101
+        }
+
+        # Calculate duration of exposure
+        duration <- deexposure - exposure
+
+        return(tibble(species, world_id, exposure, deexposure, duration))
+    }
+}
diff --git a/R/extract_climate_data.R b/R/extract_climate_data.R
new file mode 100644
index 0000000..1a58bff
--- /dev/null
+++ b/R/extract_climate_data.R
@@ -0,0 +1,16 @@
+#' Extract climate data for each grid cell
+#' @param climate_data Raster data for climate variables
+#' @param grid Grid data (100 km x 100 km)
+#' @return A tibble with climate data for each grid cell
+extract_climate_data <- function(climate_data, grid) {
+    climate <- project(climate_data, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
+    climate <- rotate(climate)
+    climate <- climate - 273.15
+
+    df <- exact_extract(climate, grid, fun = "mean", weights = "area")
+    df <- as_tibble(df) %>%
+        mutate(world_id = grid$world_id) %>%
+        relocate(world_id)
+
+    return(df)
+}
diff --git a/R/get_niche_limits.R b/R/get_niche_limits.R
new file mode 100644
index 0000000..cd8e409
--- /dev/null
+++ b/R/get_niche_limits.R
@@ -0,0 +1,63 @@
+#' Compute the thermal niche limits for each species
+#' @param species_ranges List of grid cell IDs for each species
+#' @param climate_df Tibble of climate data
+#' @return A tibble with niche_max and niche_min for each species
+get_niche_limits <- function(species_ranges, climate_df) {
+  # Filter climate data for the species ranges
+  data <- climate_df %>%
+    filter(world_id %in% species_ranges) %>%
+    select(-world_id) %>%
+    na.omit()
+
+  # Return NA when no data is available
+  if (nrow(data) == 0) {
+    return(tibble(niche_max = NA, niche_min = NA))
+  }
+
+  # Calculate mean and standard deviation
+  means <- apply(data, 1, mean)
+  sds <- apply(data, 1, sd) * 3
+
+  # Define upper and lower limits
+  upper_limit <- means + sds
+  lower_limit <- means - sds
+
+  # Remove outliers
+  upper_outliers <- sweep(data, 1, upper_limit)
+  lower_outliers <- sweep(data, 1, lower_limit)
+  data[upper_outliers > 0] <- NA
+  data[lower_outliers < 0] <- NA
+
+  # Compute max and min for each row
+  row_max <- apply(data, 1, max, na.rm = TRUE)
+  row_min <- apply(data, 1, min, na.rm = TRUE)
+
+  # Calculate overall niche limits
+  row_max_mean <- mean(row_max)
+  row_max_sd <- sd(row_max) * 3
+
+  row_min_mean <- mean(row_min)
+  row_min_sd <- sd(row_min) * 3
+
+  if (!is.na(row_max_sd)) {
+    # Handle outlier removal for max and min
+    row_max_upper <- row_max_mean + row_max_sd
+    row_max_lower <- row_max_mean - row_max_sd
+
+    row_min_upper <- row_min_mean + row_min_sd
+    row_min_lower <- row_min_mean - row_min_sd
+
+    pre_max <- row_max[which(row_max <= row_max_upper & row_max >= row_max_lower)]
+    pre_min <- row_min[which(row_min <= row_min_upper & row_min >= row_min_lower)]
+
+    niche_max <- max(pre_max)
+    niche_min <- min(pre_min)
+  } else {
+    # Fallback calculation
+    niche_max <- apply(data, 1, max, na.rm = TRUE)
+    niche_min <- apply(data, 1, min, na.rm = TRUE)
+  }
+
+  # Return results as a tibble
+  return(tibble(niche_max, niche_min))
+}
diff --git a/R/prepare_range.R b/R/prepare_range.R
new file mode 100644
index 0000000..c5b7d2c
--- /dev/null
+++ b/R/prepare_range.R
@@ -0,0 +1,37 @@
+#' Prepare range data to match the grid
+#' @param range_data Polygon data of species distribution
+#' @param grid Grid data (100 km x 100 km)
+#' @return A list of grid cell IDs for each species
+prepare_range <- function(range_data, grid) {
+    # Filter presence (extant), origin (native and reintroduced), and seasonal (resident and breeding)
+    range_filtered <- range_data %>%
+        dplyr::filter(
+            presence == 1,
+            origin %in% c(1, 2),
+            seasonal %in% c(1, 2)
+        )
+
+    # Enable parallel processing
+    plan("multisession", workers = availableCores() - 1)
+
+    res <- future_map(
+        st_geometry(range_filtered),
+        possibly(function(x) {
+            y <- st_intersects(x, grid)
+            y <- unlist(y)
+            y <- grid %>%
+                slice(y) %>%
+                pull(world_id)
+            y
+        }, quiet = TRUE),
+        .progress = TRUE
+    )
+
+    names(res) <- range_filtered$sci_name
+    res <- discard(res, is.null)
+
+    # Combine elements with the same name
+    res_final <- tapply(unlist(res, use.names = FALSE), rep(names(res), lengths(res)), FUN = c)
+
+    return(res_final)
+}

From ea4bd8ff2bb23ef3bd1ad5673ad34e93a68c5ddf Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Mon, 11 Nov 2024 22:51:57 -0800
Subject: [PATCH 2/8] Refactor: Update driver script to use reusable functions

---
 scripts/VISS_Sample_Data.R | 362 ++++---------------------------------
 1 file changed, 38 insertions(+), 324 deletions(-)

diff --git a/scripts/VISS_Sample_Data.R b/scripts/VISS_Sample_Data.R
index fff1c67..fbec661 100644
--- a/scripts/VISS_Sample_Data.R
+++ b/scripts/VISS_Sample_Data.R
@@ -1,357 +1,71 @@
 library(tidyverse)
 library(furrr)
 library(terra)
-library(exactextractr)
 library(pbapply)
-library(sf)
 library(parallel)
 
-# set the folder containing the files as the working directory
+# Source functions from the /R folder
+source("R/prepare_range.R")
+source("R/extract_climate_data.R")
+source("R/get_niche_limits.R")
+source("R/exposure.R")
+source("R/exposure_times.R")
+
+# Set the folder containing the files as the working directory
 path <- "data-raw/"
 
-# load data
+# Load input data
 historical_climate <- readRDS(paste0(path, "historical_climaate_data.rds"))
 future_climate <- readRDS(paste0(path, "future_climaate_data.rds"))
 grid <- readRDS(paste0(path, "grid.rds"))
 primates_shp <- readRDS(paste0(path, "primates_shapefiles.rds"))
 
-# historical_climate : spatraster containing monthly climate simulations of surface mean temperature from 1850 to 2014.
-# future_climate : spatraster containing monthly climate simulations  of surface mean temperature from 2015 to 2100.
-# primates_shp : distribution polygons of world primates
-# grid : 100 km x 100 km grid used to extract and match climate and distribution data
-
-############################################################
-# 1. Transform the distribution polygons to match the grid
-# The function prepare_range transform the polygons to the grid format.
-# It returns a list in which each elements contains a vector of integers,
-# which represents the IDs of the grid cells that overlap with the polygons
-
-# function
-prepare_range <- function(range_data, grid){
-
-
-  # filter presence (extant), origin (native and reintroduced), and seasonal (resident and breeding)
-  range_filtered <- range_data %>%
-    dplyr::filter(presence == 1,
-                  origin %in% c(1,2),
-                  seasonal %in% c(1,2))
-
-
-  plan("multisession", workers = availableCores() - 1)
-
-  res <- future_map(st_geometry(range_filtered), possibly(function(x){
-    y <- st_intersects(x, grid)
-    y <- unlist(y)
-    y <- grid %>%
-      slice(y) %>%
-      pull(world_id)
-    y
-
-  }, quiet = T), .progress = TRUE)
-
-
-  names(res) <- range_filtered$sci_name
-
-  res <- discard(res, is.null)
-
-  # combining elements with same name
-  res_final <- tapply(unlist(res, use.names = FALSE), rep(names(res), lengths(res)), FUN = c)
-
-  return(res_final)
-
-}
-
-# run
+# Prepare range data
 primates_range_data <- prepare_range(primates_shp, grid)
 
-############################################################
-# 2. Extract climate data using the grid
-# The function extract_climate_data extract the climate data associate with each cell in the grid.
-# It creates a data frame (tibble) in which each row
-# represents a grid cell, and each columns a time step
-
-# function
-extract_climate_data <- function(climate_data, grid){
-
-  climate <- project(climate_data, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
-  climate <- rotate(climate)
-  climate <- climate - 273.15
-
-  df <- exact_extract(climate, grid, fun = "mean", weights = "area")
-  df <- as_tibble(df) %>%
-    mutate(world_id = grid$world_id) %>%
-    relocate(world_id)
-
-  return(df)
-
-}
-
-# run
+# Extract climate data for historical and future datasets
 historical_climate_df <- extract_climate_data(historical_climate, grid)
 future_climate_df <- extract_climate_data(future_climate, grid)
 
-# rename columns
-colnames(historical_climate_df) <- c("world_id", 1850:2014)
-colnames(future_climate_df) <- c("world_id", 2015:2100)
-
-############################################################
-# 3. Compute the thermal limits for each species
-# The function get_niche_limits calculates upper and lower niche limits.
-# It creates a data frame with maximum and minimum niche estimates for each species
-
-# functions
-get_niche_limits <- function(species_ranges, climate_df){
-
-  data <- climate_df %>%
-    filter(world_id %in% species_ranges) %>%
-    select(-world_id) %>%
-    na.omit()
-
-  # when the range don't overlap with the climate data, return NA
-  if(nrow(data) == 0) return(tibble(niche_max = NA, niche_min = NA))
-
-  means <- apply(data, 1, mean)
-  sds <- apply(data, 1, sd) * 3
-
-  upper_limit <- means + sds
-  lower_limit <- means - sds
-
-  upper_outliers <- sweep(data ,1, upper_limit)
-  lower_outliers <- sweep(data ,1, lower_limit)
-
-  data[upper_outliers > 0] <- NA
-  data[lower_outliers < 0] <- NA
-
-  row_max <- apply(data, 1, max, na.rm = T)
-  row_min <- apply(data, 1, min, na.rm = T)
-
-  row_max_mean <- mean(row_max)
-  row_max_sd <- sd(row_max) * 3
-
-  row_min_mean <- mean(row_min)
-  row_min_sd <- sd(row_min) * 3
-
-  if(!is.na(row_max_sd)){
-
-    row_max_upper <- row_max_mean + row_max_sd
-    row_max_lower <- row_max_mean - row_max_sd
-
-    row_min_upper <- row_min_mean + row_min_sd
-    row_min_lower <- row_min_mean - row_min_sd
-
-    pre_max <- row_max[which(row_max <= row_max_upper & row_max >= row_max_lower)]
-    pre_min <- row_min[which(row_min <= row_min_upper & row_min >= row_min_lower)]
-
-    niche_max <- max(pre_max)
-    niche_min <- min(pre_min)
-
-    res <- data.frame(niche_max,niche_min)
-
-  } else {
-
-    niche_max <- apply(data, 1, max, na.rm = T)
-    niche_min <- apply(data, 1, min, na.rm = T)
-
-    res <- data.frame(niche_max,niche_min)
-
-  }
-
-  return(as_tibble(res))
-
-}
-
-# run
+# Compute thermal niche limits for each species
 plan("multisession", workers = availableCores() - 1)
-niche_limits <- future_map_dfr(primates_range_data, ~ get_niche_limits(.x, historical_climate_df), .id = "species", .progress = T)
-
-
-############################################################
-# 4. Calculate exposure
-# The function exposure calculates the years in which climate change exceeds the
-# niche limits of the species. The code produce data frames in which rows
-# represent the species occurring in a grid cell and columns are the years in the time series.
-# The cell is assigned with 1 if the climate is suitable for the species in a given year (i.e. below the maximum niche limit),
-# and 0 if the climate is unsuitable (above the maximum niche limit)
-
-# function
-exposure <- function(data, species_range, climate_data, niche){
-
-  spp_data <- species_range[[data]]
-  spp_name <- names(species_range)[[data]]
-
-  spp_matrix <- climate_data %>%
-    filter(world_id %in% spp_data) %>%
-    na.omit()
-
-  spp_niche <- niche %>%
-    filter(species %in% spp_name)
-
-
-  spp_matrix <- spp_matrix %>%
-    mutate(across(2:ncol(spp_matrix), ~ case_when(. <= spp_niche$niche_max ~ 1,
-                                                  . > spp_niche$niche_max ~ 0)))
-
-  spp_matrix$species <- spp_name
-  spp_matrix <- spp_matrix %>%
-    relocate(species)
-
-  return(spp_matrix)
-
-}
-
-# run
-plan("multisession", workers = availableCores() - 1)
-exposure_list <- future_map(1:length(primates_range_data), ~ exposure(.x, primates_range_data, future_climate_df, niche_limits), .progress = T)
-names(exposure_list) <- names(primates_range_data)
-
-
-############################################################
-# 4. Calculate exposure
-# The function exposure calculates the years in which climate change exceeds the
-# niche limits of the species. The code produce data frames in which rows
-# represent the species occurring in a grid cell and columns are the years in the time series.
-# The cell is assigned with 1 if the climate is suitable for the species in a given year (i.e. below the maximum niche limit),
-# and 0 if the climate is unsuitable (above the maximum niche limit)
-
-# function
-exposure <- function(data, species_range, climate_data, niche){
-
-  spp_data <- species_range[[data]]
-  spp_name <- names(species_range)[[data]]
-
-  spp_matrix <- climate_data %>%
-    filter(world_id %in% spp_data) %>%
-    na.omit()
-
-  spp_niche <- niche %>%
-    filter(species %in% spp_name)
-
-
-  spp_matrix <- spp_matrix %>%
-    mutate(across(2:ncol(spp_matrix), ~ case_when(. <= spp_niche$niche_max ~ 1,
-                                                  . > spp_niche$niche_max ~ 0)))
-
-  spp_matrix$species <- spp_name
-  spp_matrix <- spp_matrix %>%
-    relocate(species)
-
-  return(spp_matrix)
-
-}
-
-# run
-plan("multisession", workers = availableCores() - 1)
-exposure_list <- future_map(1:length(primates_range_data), ~ exposure(.x, primates_range_data, future_climate_df, niche_limits), .progress = T)
+niche_limits <- future_map_dfr(
+  primates_range_data,
+  ~ get_niche_limits(.x, historical_climate_df),
+  .id = "species",
+  .progress = TRUE
+)
+
+# Calculate exposure for each species and grid cell
+exposure_list <- future_map(
+  1:length(primates_range_data),
+  ~ exposure(.x, primates_range_data, future_climate_df, niche_limits),
+  .progress = TRUE
+)
 names(exposure_list) <- names(primates_range_data)
 
-
-############################################################
-# 5. Calculate exposure times
-# The function calculates in which year exposure occurs. A population (i.e. a species occurrence in a grid cell)
-# is classified as exposed when the temperature exceeds the upper niche limit for at least
-# five consecutive years. The function returns a data frame containing the species name,
-# the ID of the grid cell (world_id), the year of exposure, the year of de-exposure, and duration of exposure.
-# To become de-exposed, the same species must experience five consecutive years under
-# temperatures within the thermal limits.
-
-# function
-exposure_times <- function(data, original.state, consecutive.elements){
-
-  species <- data[1]
-  world_id <- data[2]
-
-  n <- as.numeric(data[-c(1,2)])
-
-  # Calculate shift sequences
-  rle_x <- data.frame(unclass(rle(n)))
-
-  # Add year
-  rle_x$year <- 2015 + cumsum(rle_x$lengths) - rle_x$lengths
-
-  # Select only shifts with n or more consecuitve elements
-  rle_x <- rle_x[rle_x$lengths >= consecutive.elements,]
-
-  # Add line with original state
-  rle_x <- rbind(c(1, original.state, 2000), rle_x)
-
-  # Remove lines with shifts that are not separated for n consecutive elements
-  rle_x <- rle_x[c(TRUE, diff(rle_x$values) != 0),]
-
-  # Remove first line because the first line is either the original state
-  # or the same value as the original state
-  rle_x <- rle_x[-1,]
-
-  # if there are no rows in rle_x, it means no exposure
-  if(nrow(rle_x) == 0) {
-
-    exposure <- NA
-    deexposure <- NA
-    duration <- NA
-
-    return(tibble(species, world_id, exposure, deexposure, duration))
-
-  }
-
-
-  # if the only value in x$values is 0, it means that there was a single exposure event
-  # with no de-exposure
-
-  if(length(unique(rle_x$values)) == 1){
-    if(unique(rle_x$values) == 0){
-
-      exposure <- rle_x$year[1]
-      deexposure <- 2101 # if deexposure = 2101 it means that deexposure did not occur
-      duration <- deexposure - exposure
-      return(tibble(species, world_id, exposure, deexposure, duration))
-    }
-  }
-
-  # the remaining data will always have 0 and 1 on rle_x$values
-  if(length(unique(rle_x$values)) == 2){
-
-    exposure <- rle_x %>%
-      filter(values == 0) %>%
-      pull(year)
-
-    deexposure <- rle_x %>%
-      filter(values == 1) %>%
-      pull(year)
-
-    # if(length(deexposure) == 0) deexposure <- 2201
-    if(length(exposure) > length(deexposure))  deexposure[length(exposure)] <- 2101
-
-    duration <- deexposure - exposure
-
-    return(tibble(species, world_id, exposure, deexposure, duration))
-
-  }
-}
-
-# run
-exposure_df <- exposure_list %>%
-  bind_rows() %>%
+# Combine exposure data into a single data frame
+exposure_df <- bind_rows(exposure_list) %>%
   mutate(sum = rowSums(select(., starts_with("2")))) %>%
-  filter(sum < 82) %>%  # Select only cells with < 82 suitable years (>= 82 years means no exposure).This is done to improve computational time by avoiding calcluting exposure for species that are not exposed.
+  filter(sum < 82) %>%  # Exclude species with no exposure
   select(-sum)
 
+# Calculate exposure times (using parallel processing)
 cl <- makeCluster(availableCores() - 1)
 clusterEvalQ(cl, library(dplyr))
 clusterExport(cl, "exposure_times")
 
+res_final <- pbapply(
+  X = exposure_df,
+  MARGIN = 1,
+  FUN = function(x) exposure_times(data = x, original.state = 1, consecutive.elements = 5),
+  cl = cl
+)
 
+stopCluster(cl)  # Stop the parallel cluster
 
-res_final <- pbapply(X = exposure_df,
-                     MARGIN = 1,
-                     FUN = function(x) exposure_times(data = x, original.state = 1, consecutive.elements = 5),
-                     cl = cl)
-
-
-res_final <- res_final %>%
-  bind_rows() %>%
+# Combine results and clean up
+res_final <- bind_rows(res_final) %>%
   na.omit()
 
-
-stopCluster(cl)
-
-# Final data frame with exposure times for each species at each grid cell
 res_final

From 4bfb2e556e00532fb4726b8ebb31ae48eb2a0d3a Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Mon, 11 Nov 2024 23:37:04 -0800
Subject: [PATCH 3/8] refactor: update imports, documentation, and package
 metadata

---
 DESCRIPTION                 | 10 +++++++++-
 NAMESPACE                   | 21 +++++++++++++++++++++
 R/exposure.R                | 12 +++++++-----
 R/exposure_times.R          | 13 ++++++++-----
 R/extract_climate_data.R    | 11 ++++++++---
 R/get_niche_limits.R        | 10 +++++++---
 R/globals.R                 |  3 +++
 R/prepare_range.R           | 12 +++++++++---
 man/exposure.Rd             | 23 +++++++++++++++++++++++
 man/exposure_times.Rd       | 21 +++++++++++++++++++++
 man/extract_climate_data.Rd | 19 +++++++++++++++++++
 man/get_niche_limits.Rd     | 19 +++++++++++++++++++
 man/prepare_range.Rd        | 19 +++++++++++++++++++
 13 files changed, 173 insertions(+), 20 deletions(-)
 create mode 100644 R/globals.R
 create mode 100644 man/exposure.Rd
 create mode 100644 man/exposure_times.Rd
 create mode 100644 man/extract_climate_data.Rd
 create mode 100644 man/get_niche_limits.Rd
 create mode 100644 man/prepare_range.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index 5b56788..66d8fb1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -9,7 +9,15 @@ Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 Imports:
-    magrittr
+    magrittr,
+    dplyr,
+    purrr,
+    sf,
+    tibble,
+    terra,
+    exactextractr,
+    future,
+    furrr
 Suggests:
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index ea39623..fbdcf74 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,4 +1,25 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
+importFrom(dplyr,across)
+importFrom(dplyr,case_when)
+importFrom(dplyr,filter)
+importFrom(dplyr,mutate)
+importFrom(dplyr,pull)
+importFrom(dplyr,relocate)
+importFrom(dplyr,select)
+importFrom(exactextractr,exact_extract)
+importFrom(furrr,future_map)
+importFrom(future,availableCores)
+importFrom(future,plan)
 importFrom(magrittr,"%>%")
+importFrom(purrr,discard)
+importFrom(purrr,possibly)
+importFrom(sf,st_geometry)
+importFrom(sf,st_intersects)
+importFrom(stats,na.omit)
+importFrom(stats,sd)
+importFrom(terra,project)
+importFrom(terra,rotate)
+importFrom(tibble,as_tibble)
+importFrom(tibble,tibble)
diff --git a/R/exposure.R b/R/exposure.R
index a362cea..393a22e 100644
--- a/R/exposure.R
+++ b/R/exposure.R
@@ -1,9 +1,11 @@
-#' Calculate exposure for each species and grid cell
-#' @param data Index of the species
+#' Calculate species exposure to climate changes
+#'
+#' @param data Data for a single species
 #' @param species_range List of grid cell IDs for each species
-#' @param climate_data Climate data tibble
-#' @param niche Niche limits tibble
-#' @return A data frame with exposure results
+#' @param climate_data Data frame of climate data by grid cell
+#' @param niche Niche limits for each species
+#' @return A data frame with exposure data
+#' @importFrom dplyr filter mutate across relocate case_when
 exposure <- function(data, species_range, climate_data, niche) {
   # Get data for the current species
   spp_data <- species_range[[data]]
diff --git a/R/exposure_times.R b/R/exposure_times.R
index d07899a..7239f31 100644
--- a/R/exposure_times.R
+++ b/R/exposure_times.R
@@ -1,8 +1,11 @@
-#' Calculate exposure and de-exposure times for species
-#' @param data A row of exposure data (species, world_id, and years)
-#' @param original.state Initial exposure state (0 or 1)
-#' @param consecutive.elements Number of consecutive years for state change
-#' @return A tibble with exposure times and durations for the species
+#' Calculate exposure times for each species
+#'
+#' @param data A row of exposure data
+#' @param original.state Initial exposure state
+#' @param consecutive.elements Minimum consecutive years for state change
+#' @return A tibble with exposure and de-exposure times
+#' @importFrom dplyr filter pull
+#' @importFrom tibble tibble
 exposure_times <- function(data, original.state, consecutive.elements) {
     # Extract species and world_id
     species <- data[1]
diff --git a/R/extract_climate_data.R b/R/extract_climate_data.R
index 1a58bff..15f10e8 100644
--- a/R/extract_climate_data.R
+++ b/R/extract_climate_data.R
@@ -1,7 +1,12 @@
 #' Extract climate data for each grid cell
-#' @param climate_data Raster data for climate variables
-#' @param grid Grid data (100 km x 100 km)
-#' @return A tibble with climate data for each grid cell
+#'
+#' @param climate_data Raster data of climate variables
+#' @param grid Grid data frame
+#' @return A tibble with climate data by grid cell
+#' @importFrom dplyr mutate relocate
+#' @importFrom tibble as_tibble
+#' @importFrom terra project rotate
+#' @importFrom exactextractr exact_extract
 extract_climate_data <- function(climate_data, grid) {
     climate <- project(climate_data, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
     climate <- rotate(climate)
diff --git a/R/get_niche_limits.R b/R/get_niche_limits.R
index cd8e409..74728de 100644
--- a/R/get_niche_limits.R
+++ b/R/get_niche_limits.R
@@ -1,7 +1,11 @@
-#' Compute the thermal niche limits for each species
+#' Calculate thermal niche limits for each species
+#'
 #' @param species_ranges List of grid cell IDs for each species
-#' @param climate_df Tibble of climate data
-#' @return A tibble with niche_max and niche_min for each species
+#' @param climate_df Data frame of climate data by grid cell
+#' @return A tibble with upper and lower niche limits
+#' @importFrom dplyr filter select
+#' @importFrom tibble tibble
+#' @importFrom stats na.omit sd
 get_niche_limits <- function(species_ranges, climate_df) {
   # Filter climate data for the species ranges
   data <- climate_df %>%
diff --git a/R/globals.R b/R/globals.R
new file mode 100644
index 0000000..272f9fc
--- /dev/null
+++ b/R/globals.R
@@ -0,0 +1,3 @@
+utils::globalVariables(c(
+  "world_id", "species", "values", "year", "presence", "origin", "seasonal"
+))
diff --git a/R/prepare_range.R b/R/prepare_range.R
index c5b7d2c..2a7cf36 100644
--- a/R/prepare_range.R
+++ b/R/prepare_range.R
@@ -1,7 +1,13 @@
 #' Prepare range data to match the grid
-#' @param range_data Polygon data of species distribution
-#' @param grid Grid data (100 km x 100 km)
-#' @return A list of grid cell IDs for each species
+#' @param range_data Data frame of species ranges
+#' @param grid Grid data frame for spatial matching
+#' @return A list of matched ranges
+#' @importFrom dplyr filter mutate select pull
+#' @importFrom sf st_geometry st_intersects
+#' @importFrom purrr discard possibly
+#' @importFrom future plan availableCores
+#' @importFrom furrr future_map
+
 prepare_range <- function(range_data, grid) {
     # Filter presence (extant), origin (native and reintroduced), and seasonal (resident and breeding)
     range_filtered <- range_data %>%
diff --git a/man/exposure.Rd b/man/exposure.Rd
new file mode 100644
index 0000000..9f1162f
--- /dev/null
+++ b/man/exposure.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/exposure.R
+\name{exposure}
+\alias{exposure}
+\title{Calculate species exposure to climate changes}
+\usage{
+exposure(data, species_range, climate_data, niche)
+}
+\arguments{
+\item{data}{Data for a single species}
+
+\item{species_range}{List of grid cell IDs for each species}
+
+\item{climate_data}{Data frame of climate data by grid cell}
+
+\item{niche}{Niche limits for each species}
+}
+\value{
+A data frame with exposure data
+}
+\description{
+Calculate species exposure to climate changes
+}
diff --git a/man/exposure_times.Rd b/man/exposure_times.Rd
new file mode 100644
index 0000000..9cbf679
--- /dev/null
+++ b/man/exposure_times.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/exposure_times.R
+\name{exposure_times}
+\alias{exposure_times}
+\title{Calculate exposure times for each species}
+\usage{
+exposure_times(data, original.state, consecutive.elements)
+}
+\arguments{
+\item{data}{A row of exposure data}
+
+\item{original.state}{Initial exposure state}
+
+\item{consecutive.elements}{Minimum consecutive years for state change}
+}
+\value{
+A tibble with exposure and de-exposure times
+}
+\description{
+Calculate exposure times for each species
+}
diff --git a/man/extract_climate_data.Rd b/man/extract_climate_data.Rd
new file mode 100644
index 0000000..f6cd081
--- /dev/null
+++ b/man/extract_climate_data.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/extract_climate_data.R
+\name{extract_climate_data}
+\alias{extract_climate_data}
+\title{Extract climate data for each grid cell}
+\usage{
+extract_climate_data(climate_data, grid)
+}
+\arguments{
+\item{climate_data}{Raster data of climate variables}
+
+\item{grid}{Grid data frame}
+}
+\value{
+A tibble with climate data by grid cell
+}
+\description{
+Extract climate data for each grid cell
+}
diff --git a/man/get_niche_limits.Rd b/man/get_niche_limits.Rd
new file mode 100644
index 0000000..277fcac
--- /dev/null
+++ b/man/get_niche_limits.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/get_niche_limits.R
+\name{get_niche_limits}
+\alias{get_niche_limits}
+\title{Calculate thermal niche limits for each species}
+\usage{
+get_niche_limits(species_ranges, climate_df)
+}
+\arguments{
+\item{species_ranges}{List of grid cell IDs for each species}
+
+\item{climate_df}{Data frame of climate data by grid cell}
+}
+\value{
+A tibble with upper and lower niche limits
+}
+\description{
+Calculate thermal niche limits for each species
+}
diff --git a/man/prepare_range.Rd b/man/prepare_range.Rd
new file mode 100644
index 0000000..772780d
--- /dev/null
+++ b/man/prepare_range.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prepare_range.R
+\name{prepare_range}
+\alias{prepare_range}
+\title{Prepare range data to match the grid}
+\usage{
+prepare_range(range_data, grid)
+}
+\arguments{
+\item{range_data}{Data frame of species ranges}
+
+\item{grid}{Grid data frame for spatial matching}
+}
+\value{
+A list of matched ranges
+}
+\description{
+Prepare range data to match the grid
+}

From 45a52c1a181d670c3bf0f50f936834be0591564e Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Mon, 18 Nov 2024 15:10:10 -0800
Subject: [PATCH 4/8] Export all main functions for user accessibility

---
 NAMESPACE                | 5 +++++
 R/exposure.R             | 1 +
 R/exposure_times.R       | 1 +
 R/extract_climate_data.R | 1 +
 R/get_niche_limits.R     | 1 +
 R/prepare_range.R        | 2 +-
 6 files changed, 10 insertions(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index fbdcf74..fa9f1ba 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
+export(exposure)
+export(exposure_times)
+export(extract_climate_data)
+export(get_niche_limits)
+export(prepare_range)
 importFrom(dplyr,across)
 importFrom(dplyr,case_when)
 importFrom(dplyr,filter)
diff --git a/R/exposure.R b/R/exposure.R
index 393a22e..4aadd7c 100644
--- a/R/exposure.R
+++ b/R/exposure.R
@@ -6,6 +6,7 @@
 #' @param niche Niche limits for each species
 #' @return A data frame with exposure data
 #' @importFrom dplyr filter mutate across relocate case_when
+#' @export
 exposure <- function(data, species_range, climate_data, niche) {
   # Get data for the current species
   spp_data <- species_range[[data]]
diff --git a/R/exposure_times.R b/R/exposure_times.R
index 7239f31..4d87216 100644
--- a/R/exposure_times.R
+++ b/R/exposure_times.R
@@ -6,6 +6,7 @@
 #' @return A tibble with exposure and de-exposure times
 #' @importFrom dplyr filter pull
 #' @importFrom tibble tibble
+#' @export
 exposure_times <- function(data, original.state, consecutive.elements) {
     # Extract species and world_id
     species <- data[1]
diff --git a/R/extract_climate_data.R b/R/extract_climate_data.R
index 15f10e8..0bc3a48 100644
--- a/R/extract_climate_data.R
+++ b/R/extract_climate_data.R
@@ -7,6 +7,7 @@
 #' @importFrom tibble as_tibble
 #' @importFrom terra project rotate
 #' @importFrom exactextractr exact_extract
+#' @export
 extract_climate_data <- function(climate_data, grid) {
     climate <- project(climate_data, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
     climate <- rotate(climate)
diff --git a/R/get_niche_limits.R b/R/get_niche_limits.R
index 74728de..59ac8b9 100644
--- a/R/get_niche_limits.R
+++ b/R/get_niche_limits.R
@@ -6,6 +6,7 @@
 #' @importFrom dplyr filter select
 #' @importFrom tibble tibble
 #' @importFrom stats na.omit sd
+#' @export
 get_niche_limits <- function(species_ranges, climate_df) {
   # Filter climate data for the species ranges
   data <- climate_df %>%
diff --git a/R/prepare_range.R b/R/prepare_range.R
index 2a7cf36..890b40c 100644
--- a/R/prepare_range.R
+++ b/R/prepare_range.R
@@ -7,7 +7,7 @@
 #' @importFrom purrr discard possibly
 #' @importFrom future plan availableCores
 #' @importFrom furrr future_map
-
+#' @export
 prepare_range <- function(range_data, grid) {
     # Filter presence (extant), origin (native and reintroduced), and seasonal (resident and breeding)
     range_filtered <- range_data %>%

From 0a4c76fbcdb3e015c179eb4afe64b1217b587ac5 Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Mon, 18 Nov 2024 15:51:41 -0800
Subject: [PATCH 5/8] Remove source lines from VISS_Sample_Data.R and fix
 linting issues

---
 NAMESPACE                  |   2 +-
 R/exposure_times.R         | 113 ++++++++++++++++++++++---------------
 R/extract_climate_data.R   |  19 ++++---
 R/get_niche_limits.R       |   9 ++-
 R/prepare_range.R          |  56 +++++++++---------
 makefile.R                 |   1 -
 man/exposure_times.Rd      |   6 +-
 scripts/VISS_Sample_Data.R |  12 +---
 8 files changed, 120 insertions(+), 98 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index fa9f1ba..5c39107 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,7 +2,6 @@
 
 export("%>%")
 export(exposure)
-export(exposure_times)
 export(extract_climate_data)
 export(get_niche_limits)
 export(prepare_range)
@@ -13,6 +12,7 @@ importFrom(dplyr,mutate)
 importFrom(dplyr,pull)
 importFrom(dplyr,relocate)
 importFrom(dplyr,select)
+importFrom(dplyr,slice)
 importFrom(exactextractr,exact_extract)
 importFrom(furrr,future_map)
 importFrom(future,availableCores)
diff --git a/R/exposure_times.R b/R/exposure_times.R
index 4d87216..13c2f48 100644
--- a/R/exposure_times.R
+++ b/R/exposure_times.R
@@ -1,69 +1,88 @@
 #' Calculate exposure times for each species
 #'
 #' @param data A row of exposure data
-#' @param original.state Initial exposure state
-#' @param consecutive.elements Minimum consecutive years for state change
+#' @param original_state Initial exposure state
+#' @param consecutive_elements Minimum consecutive years for state change
 #' @return A tibble with exposure and de-exposure times
 #' @importFrom dplyr filter pull
 #' @importFrom tibble tibble
-#' @export
-exposure_times <- function(data, original.state, consecutive.elements) {
-    # Extract species and world_id
-    species <- data[1]
-    world_id <- data[2]
+exposure_times <- function(data, original_state, consecutive_elements) {
+  # Extract species and world_id
+  species <- data[1]
+  world_id <- data[2]
 
-    # Extract year data as numeric vector
-    n <- as.numeric(data[-c(1, 2)])
+  # Extract year data as numeric vector
+  n <- as.numeric(data[-c(1, 2)])
 
-    # Calculate shift sequences
-    rle_x <- data.frame(unclass(rle(n)))
+  # Calculate shift sequences
+  rle_x <- data.frame(unclass(rle(n)))
 
-    # Add year column to represent time steps
-    rle_x$year <- 2015 + cumsum(rle_x$lengths) - rle_x$lengths
+  # Add year column to represent time steps
+  rle_x$year <- 2015 + cumsum(rle_x$lengths) - rle_x$lengths
 
-    # Filter sequences with sufficient consecutive elements
-    rle_x <- rle_x[rle_x$lengths >= consecutive.elements,]
+  # Filter sequences with sufficient consecutive elements
+  rle_x <- rle_x[rle_x$lengths >= consecutive_elements, ]
 
-    # Add a line for the original state to ensure valid transitions
-    rle_x <- rbind(c(1, original.state, 2000), rle_x)
+  # Add a line for the original state to ensure valid transitions
+  rle_x <- rbind(c(1, original_state, 2000), rle_x)
 
-    # Remove unnecessary state repetitions
-    rle_x <- rle_x[c(TRUE, diff(rle_x$values) != 0),]
+  # Remove unnecessary state repetitions
+  rle_x <- rle_x[c(TRUE, diff(rle_x$values) != 0), ]
 
-    # Remove the first line (original state or duplicate)
-    rle_x <- rle_x[-1,]
+  # Remove the first line (original state or duplicate)
+  rle_x <- rle_x[-1, ]
 
-    # Handle cases with no valid exposure sequences
-    if (nrow(rle_x) == 0) {
-        return(tibble(species, world_id, exposure = NA, deexposure = NA, duration = NA))
-    }
+  # Handle cases with no valid exposure sequences
+  if (nrow(rle_x) == 0) {
+    return(tibble(
+      species = species,
+      world_id = world_id,
+      exposure = NA,
+      deexposure = NA,
+      duration = NA
+    ))
+  }
 
-    # Handle cases where all values are 0 (exposure with no de-exposure)
-    if (length(unique(rle_x$values)) == 1 && unique(rle_x$values) == 0) {
-        exposure <- rle_x$year[1]
-        deexposure <- 2101  # Indicates de-exposure did not occur
-        duration <- deexposure - exposure
-        return(tibble(species, world_id, exposure, deexposure, duration))
-    }
+  # Handle cases where all values are 0 (exposure with no de-exposure)
+  if (length(unique(rle_x$values)) == 1 && unique(rle_x$values) == 0) {
+    exposure <- rle_x$year[1]
+    deexposure <- 2101 # Indicates de-exposure did not occur
+    duration <- deexposure - exposure
+    return(tibble(
+      species = species,
+      world_id = world_id,
+      exposure = exposure,
+      deexposure = deexposure,
+      duration = duration
+    ))
+  }
+
+  # Handle cases with both exposure (0) and de-exposure (1)
+  if (length(unique(rle_x$values)) == 2) {
+    exposure <- rle_x %>%
+      filter(values == 0) %>%
+      pull(year)
 
-    # Handle cases with both exposure (0) and de-exposure (1)
-    if (length(unique(rle_x$values)) == 2) {
-        exposure <- rle_x %>%
-            filter(values == 0) %>%
-            pull(year)
+    deexposure <- rle_x %>%
+      filter(values == 1) %>%
+      pull(year)
 
-        deexposure <- rle_x %>%
-            filter(values == 1) %>%
-            pull(year)
+    # If there are more exposures than deexposures,
+    # add a placeholder for deexposure.
+    if (length(exposure) > length(deexposure)) {
+      deexposure[length(exposure)] <- 2101
 
-        # If there are more exposures than deexposures, add a placeholder for deexposure
-        if (length(exposure) > length(deexposure)) {
-            deexposure[length(exposure)] <- 2101
-        }
 
-        # Calculate duration of exposure
-        duration <- deexposure - exposure
+      # Calculate the duration of exposure
+      duration <- deexposure - exposure
 
-        return(tibble(species, world_id, exposure, deexposure, duration))
+      return(tibble(
+        species = species,
+        world_id = world_id,
+        exposure = exposure,
+        deexposure = deexposure,
+        duration = duration
+      ))
     }
+  }
 }
diff --git a/R/extract_climate_data.R b/R/extract_climate_data.R
index 0bc3a48..cd8e0ef 100644
--- a/R/extract_climate_data.R
+++ b/R/extract_climate_data.R
@@ -9,14 +9,17 @@
 #' @importFrom exactextractr exact_extract
 #' @export
 extract_climate_data <- function(climate_data, grid) {
-    climate <- project(climate_data, "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
-    climate <- rotate(climate)
-    climate <- climate - 273.15
+  climate <- project(
+    climate_data,
+    "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
+  )
+  climate <- rotate(climate)
+  climate <- climate - 273.15
 
-    df <- exact_extract(climate, grid, fun = "mean", weights = "area")
-    df <- as_tibble(df) %>%
-        mutate(world_id = grid$world_id) %>%
-        relocate(world_id)
+  df <- exact_extract(climate, grid, fun = "mean", weights = "area")
+  df <- as_tibble(df) %>%
+    mutate(world_id = grid$world_id) %>%
+    relocate(world_id)
 
-    return(df)
+  return(df)
 }
diff --git a/R/get_niche_limits.R b/R/get_niche_limits.R
index 59ac8b9..a04d47b 100644
--- a/R/get_niche_limits.R
+++ b/R/get_niche_limits.R
@@ -52,8 +52,13 @@ get_niche_limits <- function(species_ranges, climate_df) {
     row_min_upper <- row_min_mean + row_min_sd
     row_min_lower <- row_min_mean - row_min_sd
 
-    pre_max <- row_max[which(row_max <= row_max_upper & row_max >= row_max_lower)]
-    pre_min <- row_min[which(row_min <= row_min_upper & row_min >= row_min_lower)]
+    pre_max <- row_max[
+      which(row_max <= row_max_upper & row_max >= row_max_lower)
+    ]
+
+    pre_min <- row_min[
+      which(row_min <= row_min_upper & row_min >= row_min_lower)
+    ]
 
     niche_max <- max(pre_max)
     niche_min <- min(pre_min)
diff --git a/R/prepare_range.R b/R/prepare_range.R
index 890b40c..7a2a943 100644
--- a/R/prepare_range.R
+++ b/R/prepare_range.R
@@ -2,42 +2,44 @@
 #' @param range_data Data frame of species ranges
 #' @param grid Grid data frame for spatial matching
 #' @return A list of matched ranges
-#' @importFrom dplyr filter mutate select pull
+#' @importFrom dplyr filter mutate select pull slice
 #' @importFrom sf st_geometry st_intersects
 #' @importFrom purrr discard possibly
 #' @importFrom future plan availableCores
 #' @importFrom furrr future_map
 #' @export
 prepare_range <- function(range_data, grid) {
-    # Filter presence (extant), origin (native and reintroduced), and seasonal (resident and breeding)
-    range_filtered <- range_data %>%
-        dplyr::filter(
-            presence == 1,
-            origin %in% c(1, 2),
-            seasonal %in% c(1, 2)
-        )
+  # Filter presence (extant), origin (native and reintroduced),
+  # and seasonal (resident and breeding)
+  range_filtered <- range_data %>%
+    dplyr::filter(
+      presence == 1,
+      origin %in% c(1, 2),
+      seasonal %in% c(1, 2)
+    )
 
-    # Enable parallel processing
-    plan("multisession", workers = availableCores() - 1)
+  # Enable parallel processing
+  plan("multisession", workers = availableCores() - 1)
 
-    res <- future_map(
-        st_geometry(range_filtered),
-        possibly(function(x) {
-            y <- st_intersects(x, grid)
-            y <- unlist(y)
-            y <- grid %>%
-                slice(y) %>%
-                pull(world_id)
-            y
-        }, quiet = TRUE),
-        .progress = TRUE
-    )
+  res <- future_map(
+    st_geometry(range_filtered),
+    possibly(function(x) {
+      y <- st_intersects(x, grid)
+      y <- unlist(y)
+      y <- grid %>%
+        slice(y) %>%
+        pull(world_id)
+      y
+    }, quiet = TRUE),
+    .progress = TRUE
+  )
 
-    names(res) <- range_filtered$sci_name
-    res <- discard(res, is.null)
+  names(res) <- range_filtered$sci_name
+  res <- discard(res, is.null)
 
-    # Combine elements with the same name
-    res_final <- tapply(unlist(res, use.names = FALSE), rep(names(res), lengths(res)), FUN = c)
+  # Combine elements with the same name
+  res_final <- tapply(unlist(res, use.names = FALSE),
+                      rep(names(res), lengths(res)), FUN = c)
 
-    return(res_final)
+  return(res_final)
 }
diff --git a/makefile.R b/makefile.R
index f6bdca8..9c502a2 100644
--- a/makefile.R
+++ b/makefile.R
@@ -1,4 +1,3 @@
-
 #### Template of a master script for the project ####
 
 ## Using an advanced tool like {drake} or {targets} is recommended,
diff --git a/man/exposure_times.Rd b/man/exposure_times.Rd
index 9cbf679..6261e2b 100644
--- a/man/exposure_times.Rd
+++ b/man/exposure_times.Rd
@@ -4,14 +4,14 @@
 \alias{exposure_times}
 \title{Calculate exposure times for each species}
 \usage{
-exposure_times(data, original.state, consecutive.elements)
+exposure_times(data, original_state, consecutive_elements)
 }
 \arguments{
 \item{data}{A row of exposure data}
 
-\item{original.state}{Initial exposure state}
+\item{original_state}{Initial exposure state}
 
-\item{consecutive.elements}{Minimum consecutive years for state change}
+\item{consecutive_elements}{Minimum consecutive years for state change}
 }
 \value{
 A tibble with exposure and de-exposure times
diff --git a/scripts/VISS_Sample_Data.R b/scripts/VISS_Sample_Data.R
index fbec661..290ec39 100644
--- a/scripts/VISS_Sample_Data.R
+++ b/scripts/VISS_Sample_Data.R
@@ -3,13 +3,7 @@ library(furrr)
 library(terra)
 library(pbapply)
 library(parallel)
-
-# Source functions from the /R folder
-source("R/prepare_range.R")
-source("R/extract_climate_data.R")
-source("R/get_niche_limits.R")
-source("R/exposure.R")
-source("R/exposure_times.R")
+library(biodiversityhorizons)
 
 # Set the folder containing the files as the working directory
 path <- "data-raw/"
@@ -47,7 +41,7 @@ names(exposure_list) <- names(primates_range_data)
 # Combine exposure data into a single data frame
 exposure_df <- bind_rows(exposure_list) %>%
   mutate(sum = rowSums(select(., starts_with("2")))) %>%
-  filter(sum < 82) %>%  # Exclude species with no exposure
+  filter(sum < 82) %>% # Exclude species with no exposure
   select(-sum)
 
 # Calculate exposure times (using parallel processing)
@@ -62,7 +56,7 @@ res_final <- pbapply(
   cl = cl
 )
 
-stopCluster(cl)  # Stop the parallel cluster
+stopCluster(cl) # Stop the parallel cluster
 
 # Combine results and clean up
 res_final <- bind_rows(res_final) %>%

From 6f97b210976b6650549988d5e534f96cc9162ffe Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Wed, 20 Nov 2024 13:08:43 -0800
Subject: [PATCH 6/8] refactor: update imports, documentation, and package
 metadata

---
 DESCRIPTION                |  1 +
 NAMESPACE                  |  1 +
 R/_globals.R               |  3 ++
 R/exposure_times.R         |  1 +
 R/prepare_range.R          | 16 +++++-----
 scripts/VISS_Sample_Data.R | 63 ++++++++++++++++++++++----------------
 6 files changed, 50 insertions(+), 35 deletions(-)
 create mode 100644 R/_globals.R

diff --git a/DESCRIPTION b/DESCRIPTION
index 66d8fb1..9c2a628 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -17,6 +17,7 @@ Imports:
     terra,
     exactextractr,
     future,
+    pbapply,
     furrr
 Suggests:
     testthat (>= 3.0.0)
diff --git a/NAMESPACE b/NAMESPACE
index 5c39107..0eac0ad 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,7 @@
 
 export("%>%")
 export(exposure)
+export(exposure_times)
 export(extract_climate_data)
 export(get_niche_limits)
 export(prepare_range)
diff --git a/R/_globals.R b/R/_globals.R
new file mode 100644
index 0000000..fb283e7
--- /dev/null
+++ b/R/_globals.R
@@ -0,0 +1,3 @@
+utils::globalVariables(c("presence", "origin",
+                         "seasonal", "world_id", "sci_name",
+                         "values", "year", "lengths"))
diff --git a/R/exposure_times.R b/R/exposure_times.R
index 13c2f48..0efa05b 100644
--- a/R/exposure_times.R
+++ b/R/exposure_times.R
@@ -6,6 +6,7 @@
 #' @return A tibble with exposure and de-exposure times
 #' @importFrom dplyr filter pull
 #' @importFrom tibble tibble
+#' @export
 exposure_times <- function(data, original_state, consecutive_elements) {
   # Extract species and world_id
   species <- data[1]
diff --git a/R/prepare_range.R b/R/prepare_range.R
index 7a2a943..eebc4fa 100644
--- a/R/prepare_range.R
+++ b/R/prepare_range.R
@@ -19,23 +19,23 @@ prepare_range <- function(range_data, grid) {
     )
 
   # Enable parallel processing
-  plan("multisession", workers = availableCores() - 1)
+  future::plan("multisession", workers = future::availableCores() - 1)
 
-  res <- future_map(
-    st_geometry(range_filtered),
-    possibly(function(x) {
-      y <- st_intersects(x, grid)
+  res <- furrr::future_map(
+    sf::st_geometry(range_filtered),
+    purrr::possibly(function(x) {
+      y <- sf::st_intersects(x, grid)
       y <- unlist(y)
       y <- grid %>%
-        slice(y) %>%
-        pull(world_id)
+        dplyr::slice(y) %>%
+        dplyr::pull(world_id)
       y
     }, quiet = TRUE),
     .progress = TRUE
   )
 
   names(res) <- range_filtered$sci_name
-  res <- discard(res, is.null)
+  res <- purrr::discard(res, is.null)
 
   # Combine elements with the same name
   res_final <- tapply(unlist(res, use.names = FALSE),
diff --git a/scripts/VISS_Sample_Data.R b/scripts/VISS_Sample_Data.R
index 290ec39..1e64c2a 100644
--- a/scripts/VISS_Sample_Data.R
+++ b/scripts/VISS_Sample_Data.R
@@ -1,50 +1,56 @@
+# Load necessary libraries
 library(tidyverse)
 library(furrr)
 library(terra)
+library(exactextractr)
 library(pbapply)
+library(sf)
 library(parallel)
-library(biodiversityhorizons)
 
 # Set the folder containing the files as the working directory
 path <- "data-raw/"
 
-# Load input data
+# Load data
 historical_climate <- readRDS(paste0(path, "historical_climaate_data.rds"))
 future_climate <- readRDS(paste0(path, "future_climaate_data.rds"))
 grid <- readRDS(paste0(path, "grid.rds"))
 primates_shp <- readRDS(paste0(path, "primates_shapefiles.rds"))
 
-# Prepare range data
+# Source the functions from the /R directory
+source("R/prepare_range.R")
+source("R/extract_climate_data.R")
+source("R/get_niche_limits.R")
+source("R/exposure.R")
+source("R/exposure_times.R")
+
+# 1. Transform the distribution polygons to match the grid
 primates_range_data <- prepare_range(primates_shp, grid)
 
-# Extract climate data for historical and future datasets
+# 2. Extract climate data using the grid
 historical_climate_df <- extract_climate_data(historical_climate, grid)
 future_climate_df <- extract_climate_data(future_climate, grid)
 
-# Compute thermal niche limits for each species
+# Rename columns
+colnames(historical_climate_df) <- c("world_id", 1850:2014)
+colnames(future_climate_df) <- c("world_id", 2015:2100)
+
+# 3. Compute the thermal limits for each species
 plan("multisession", workers = availableCores() - 1)
-niche_limits <- future_map_dfr(
-  primates_range_data,
-  ~ get_niche_limits(.x, historical_climate_df),
-  .id = "species",
-  .progress = TRUE
-)
+niche_limits <- future_map_dfr(primates_range_data, ~ get_niche_limits(.x, historical_climate_df),
+                               .id = "species", .progress = TRUE)
 
-# Calculate exposure for each species and grid cell
-exposure_list <- future_map(
-  1:length(primates_range_data),
-  ~ exposure(.x, primates_range_data, future_climate_df, niche_limits),
-  .progress = TRUE
-)
+# 4. Calculate exposure
+plan("multisession", workers = availableCores() - 1)
+exposure_list <- future_map(1:length(primates_range_data), ~ exposure(.x, primates_range_data, future_climate_df, niche_limits), .progress = TRUE) #nolint
 names(exposure_list) <- names(primates_range_data)
 
-# Combine exposure data into a single data frame
-exposure_df <- bind_rows(exposure_list) %>%
+# 5. Calculate exposure times
+exposure_df <- exposure_list %>%
+  bind_rows() %>%
   mutate(sum = rowSums(select(., starts_with("2")))) %>%
-  filter(sum < 82) %>% # Exclude species with no exposure
+  filter(sum < 82) %>%  # Select only cells with less than 82 suitable years
   select(-sum)
 
-# Calculate exposure times (using parallel processing)
 cl <- makeCluster(availableCores() - 1)
 clusterEvalQ(cl, library(dplyr))
 clusterExport(cl, "exposure_times")
@@ -52,14 +58,17 @@ clusterExport(cl, "exposure_times")
 res_final <- pbapply(
   X = exposure_df,
   MARGIN = 1,
-  FUN = function(x) exposure_times(data = x, original.state = 1, consecutive.elements = 5),
+  FUN = function(x) exposure_times(data = x, original_state = 1, consecutive_elements = 5),
   cl = cl
 )
 
-stopCluster(cl) # Stop the parallel cluster
-
-# Combine results and clean up
-res_final <- bind_rows(res_final) %>%
+res_final <- res_final %>%
+  bind_rows() %>%
   na.omit()
 
-res_final
+stopCluster(cl)
+
+# Final data frame with exposure times for each species at each grid cell
+print(res_final)
+
+future::plan("sequential")

From ae5886909aea9418d3fea5ea56319161322244d8 Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Wed, 20 Nov 2024 13:30:04 -0800
Subject: [PATCH 7/8] Update globals file and Update .gitignore to exclude
 .tar.gz files

---
 .gitignore   | 1 +
 R/_globals.R | 3 ---
 R/globals.R  | 6 +++---
 3 files changed, 4 insertions(+), 6 deletions(-)
 delete mode 100644 R/_globals.R

diff --git a/.gitignore b/.gitignore
index 9287c8d..7df63eb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -29,6 +29,7 @@ share/python-wheels/
 .installed.cfg
 *.egg
 MANIFEST
+*.tar.gz
 
 # PyInstaller
 #  Usually these files are written by a python script from a template
diff --git a/R/_globals.R b/R/_globals.R
deleted file mode 100644
index fb283e7..0000000
--- a/R/_globals.R
+++ /dev/null
@@ -1,3 +0,0 @@
-utils::globalVariables(c("presence", "origin",
-                         "seasonal", "world_id", "sci_name",
-                         "values", "year", "lengths"))
diff --git a/R/globals.R b/R/globals.R
index 272f9fc..8569d43 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -1,3 +1,3 @@
-utils::globalVariables(c(
-  "world_id", "species", "values", "year", "presence", "origin", "seasonal"
-))
+utils::globalVariables(c("presence", "origin",
+                         "seasonal", "world_id", "sci_name",
+                         "values", "year", "lengths", "species"))

From 775e9f3580682538e8553ff32751131999b2a5d0 Mon Sep 17 00:00:00 2001
From: stlp <ishikakhandelwal02@gmail.com>
Date: Thu, 21 Nov 2024 09:04:43 -0800
Subject: [PATCH 8/8] fix: Resolve warnings and handle imports using usethis
 utilities

---
 DESCRIPTION                      |  3 ++-
 NAMESPACE                        |  4 ++++
 R/biodiversityhorizons-package.R | 23 +++++++++++++++++++++++
 scripts/VISS_Sample_Data.R       | 25 +++++++++++++++----------
 4 files changed, 44 insertions(+), 11 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 9c2a628..fc6939b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -18,7 +18,8 @@ Imports:
     exactextractr,
     future,
     pbapply,
-    furrr
+    furrr,
+    parallel
 Suggests:
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index 0eac0ad..714daa2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -19,6 +19,10 @@ importFrom(furrr,future_map)
 importFrom(future,availableCores)
 importFrom(future,plan)
 importFrom(magrittr,"%>%")
+importFrom(parallel,clusterEvalQ)
+importFrom(parallel,clusterExport)
+importFrom(parallel,makeCluster)
+importFrom(pbapply,pbapply)
 importFrom(purrr,discard)
 importFrom(purrr,possibly)
 importFrom(sf,st_geometry)
diff --git a/R/biodiversityhorizons-package.R b/R/biodiversityhorizons-package.R
index a65cf64..2466c05 100644
--- a/R/biodiversityhorizons-package.R
+++ b/R/biodiversityhorizons-package.R
@@ -2,5 +2,28 @@
 "_PACKAGE"
 
 ## usethis namespace: start
+#' @importFrom dplyr across
+#' @importFrom dplyr filter
+#' @importFrom dplyr mutate
+#' @importFrom dplyr pull
+#' @importFrom dplyr relocate
+#' @importFrom dplyr select
+#' @importFrom exactextractr exact_extract
+#' @importFrom furrr future_map
+#' @importFrom future availableCores
+#' @importFrom future plan
+#' @importFrom magrittr %>%
+#' @importFrom parallel clusterEvalQ
+#' @importFrom parallel clusterExport
+#' @importFrom parallel makeCluster
+#' @importFrom pbapply pbapply
+#' @importFrom purrr discard
+#' @importFrom purrr possibly
+#' @importFrom sf st_geometry
+#' @importFrom sf st_intersects
+#' @importFrom terra project
+#' @importFrom terra rotate
+#' @importFrom tibble as_tibble
+#' @importFrom tibble tibble
 ## usethis namespace: end
 NULL
diff --git a/scripts/VISS_Sample_Data.R b/scripts/VISS_Sample_Data.R
index 1e64c2a..805700a 100644
--- a/scripts/VISS_Sample_Data.R
+++ b/scripts/VISS_Sample_Data.R
@@ -16,12 +16,7 @@ future_climate <- readRDS(paste0(path, "future_climaate_data.rds"))
 grid <- readRDS(paste0(path, "grid.rds"))
 primates_shp <- readRDS(paste0(path, "primates_shapefiles.rds"))
 
-# Source the functions from the /R directory
-source("R/prepare_range.R")
-source("R/extract_climate_data.R")
-source("R/get_niche_limits.R")
-source("R/exposure.R")
-source("R/exposure_times.R")
+devtools::load_all()
 
 # 1. Transform the distribution polygons to match the grid
 primates_range_data <- prepare_range(primates_shp, grid)
@@ -36,8 +31,13 @@ colnames(future_climate_df) <- c("world_id", 2015:2100)
 
 # 3. Compute the thermal limits for each species
 plan("multisession", workers = availableCores() - 1)
-niche_limits <- future_map_dfr(primates_range_data, ~ get_niche_limits(.x, historical_climate_df),
-                               .id = "species", .progress = TRUE)
+
+niche_limits <- future_map_dfr(
+  primates_range_data,
+  ~ get_niche_limits(.x, historical_climate_df),
+  .id = "species",
+  .progress = TRUE
+)
 
 # 4. Calculate exposure
 plan("multisession", workers = availableCores() - 1)
@@ -58,8 +58,13 @@ clusterExport(cl, "exposure_times")
 res_final <- pbapply(
   X = exposure_df,
   MARGIN = 1,
-  FUN = function(x) exposure_times(data = x, original_state = 1, consecutive_elements = 5),
-  cl = cl
+  FUN = function(x) {
+    exposure_times(
+      data = x,
+      original_state = 1,
+      consecutive_elements = 5
+    )
+  }
 )
 
 res_final <- res_final %>%