refactor: scripts/VISS_Sample_Data.R

.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

DESCRIPTION

@@ -9,7 +9,17 @@ Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 Imports:
-    magrittr
+    magrittr,
+    dplyr,
+    purrr,
+    sf,
+    tibble,
+    terra,
+    exactextractr,
+    future,
+    pbapply,
+    furrr,
+    parallel
 Suggests:
     testthat (>= 3.0.0)
 Config/testthat/edition: 3

NAMESPACE

@@ -1,4 +1,35 @@
 # 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)
+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)
+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(stats,na.omit)
+importFrom(stats,sd)
+importFrom(terra,project)
+importFrom(terra,rotate)
+importFrom(tibble,as_tibble)
+importFrom(tibble,tibble)

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

R/exposure.R

@@ -0,0 +1,37 @@
+#' 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
+#' @export
+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,89 @@
+#' 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
+#' @export
+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 = 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 = 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)
+
+    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 the duration of exposure
+      duration <- deexposure - exposure
+
+      return(tibble(
+        species = species,
+        world_id = world_id,
+        exposure = exposure,
+        deexposure = deexposure,
+        duration = duration
+      ))
+    }
+  }
+}

R/extract_climate_data.R

@@ -0,0 +1,25 @@
+#' 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
+#' @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
+
+  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,73 @@
+#' 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
+#' @export
+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("presence", "origin",
+                         "seasonal", "world_id", "sci_name",
+                         "values", "year", "lengths", "species"))

R/prepare_range.R

@@ -0,0 +1,45 @@
+#' 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 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)
+    )
+
+  # Enable parallel processing
+  future::plan("multisession", workers = future::availableCores() - 1)
+
+  res <- furrr::future_map(
+    sf::st_geometry(range_filtered),
+    purrr::possibly(function(x) {
+      y <- sf::st_intersects(x, grid)
+      y <- unlist(y)
+      y <- grid %>%
+        dplyr::slice(y) %>%
+        dplyr::pull(world_id)
+      y
+    }, quiet = TRUE),
+    .progress = TRUE
+  )
+
+  names(res) <- range_filtered$sci_name
+  res <- purrr::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)
+}

makefile.R

@@ -1,4 +1,3 @@
-
 #### Template of a master script for the project ####
 
 ## Using an advanced tool like {drake} or {targets} is recommended,