refactor: scripts/VISS_Sample_Data.R

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

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)

R/exposure.R

@@ -0,0 +1,36 @@
+#' 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 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]]
+  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)
+}

R/exposure_times.R

@@ -0,0 +1,68 @@
+#' 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]
+    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))
+    }
+}

R/extract_climate_data.R

@@ -0,0 +1,21 @@
+#' Extract 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)
+    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)
+}

R/get_niche_limits.R

@@ -0,0 +1,67 @@
+#' Calculate thermal niche limits for each species
+#'
+#' @param species_ranges List of grid cell IDs 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 %>%
+    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))
+}

R/globals.R

@@ -0,0 +1,3 @@
+utils::globalVariables(c(
+  "world_id", "species", "values", "year", "presence", "origin", "seasonal"
+))

R/prepare_range.R

@@ -0,0 +1,43 @@
+#' Prepare range data to match the grid
+#' @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 %>%
+        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)
+}