Data preperation for Melanoplus niche modelling

Download occurrence data

# ==============================================================================
# Script: 01_data_prep.R
# Purpose: Pre-download and project data for the Melanoplus Shiny Apps
# ==============================================================================
library(rgbif)
library(terra)
library(geodata)
library(dplyr)
library(tidyr)
library(CoordinateCleaner)
library(jsonlite)
library(stringr)
library(httr)

# ==============================================================================
# 0. Setup Directories
# ==============================================================================
shi.dir <- "" #set your working directory here
setwd(shi.dir)
dir.create("data", showWarnings = FALSE)
dir.create("raw", showWarnings = FALSE)

# ==========================================
# 1. GBIF DOWNLOAD (Genus: Melanoplus)
# ==========================================
Melanoplus_key <- name_backbone("Melanoplus", rank = "genus")$usageKey
print(paste("Using Melanoplus Key:", Melanoplus_key))

download_request <- occ_download(
  pred("taxonKey", Melanoplus_key),
  pred("hasCoordinate", TRUE)
)

occ_download_wait(download_request)

raw_df <- occ_download_import(d)

Clean occurrence data

# ==========================================
# 2. METADATA FILTERING & COORDINATE CLEANING & TAXON GROUPING
# ==========================================
df_filtered <- raw_df %>%
  select(genus, species, decimalLongitude, decimalLatitude, basisOfRecord, year) %>%
  filter(!is.na(genus), !is.na(species), species != "") %>%
  filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
  filter(basisOfRecord != "FOSSIL_SPECIMEN")

df_filtered$decimalLongitude <- as.numeric(df_filtered$decimalLongitude)
df_filtered$decimalLatitude <- as.numeric(df_filtered$decimalLatitude)

df_clean <- clean_coordinates(
  x = df_filtered,
  lon = "decimalLongitude",
  lat = "decimalLatitude",
  species = "species",
  tests = c("capitals", "centroids", "equal", "gbif", "institutions", "zeros", "seas")
) %>% filter(.summary == TRUE)

# 1. Fetch the Top Level (The 32 Groups)
url <- "https://sfg.taxonworks.org/api/v1/otus/816542/inventory/taxonomy.json?project_token=3oerVKf82_196cIECvHYNg&max_descendants_depth=2"
response <- GET(url)
data <- fromJSON(content(response, "text", encoding = "UTF-8"))

data <- fromJSON("Melanoplus_groups.txt")

# 2. Extract the Group and Species Mappings
mapping <- data.frame(species = character(), species_group = character(), stringsAsFactors = FALSE)
groups_df <- data$descendants

for (i in 1:nrow(groups_df)) {
  # Clean HTML tags from the group name
  group_name_raw <- str_replace_all(groups_df$name[i], "<[^>]+>", "")
  
  if (str_detect(group_name_raw, "supersp\\.")) {
    clean_group_name <- paste(str_extract(group_name_raw, "(?<=supersp\\. )[a-z]+"), "group")
    
    # Extract the nested dataframe of species for this specific group
    children_df <- groups_df$descendants[[i]]
    
    if (is.data.frame(children_df) && nrow(children_df) > 0) {
      for (j in 1:nrow(children_df)) {
        sp_name_raw <- str_replace_all(children_df$name[j], "<[^>]+>", "")
        
        # Keep only "Melanoplus [species]" (drop the author/year)
        clean_sp_name <- str_extract(sp_name_raw, "^Melanoplus [a-z]+")
        
        mapping <- rbind(mapping, data.frame(species = clean_sp_name, species_group = clean_group_name))
      }
    }
  }
}

# Remove any duplicates just in case
mapping <- unique(mapping)
print(head(mapping)) # Let's print the first few rows to verify it worked!

# 3. Load your 300 GBIF species and merge
gbif_species <- read.csv("Melanopulus_species.csv") 

final_csv <- gbif_species %>%
  left_join(mapping, by = c("x" = "species")) %>%
  rename(species = x) %>%
  # If a species isn't in a TaxonWorks group, label it "Ungrouped"
  mutate(species_group = ifelse(is.na(species_group), "Ungrouped", species_group)) %>%
  # Drop the extra unnamed column if it exists
  select(species, species_group)

# 4. Save it for the Shiny app
write.csv(final_csv, "data/Melanopulus_species.csv", row.names = FALSE)
cat("Successfully mapped", nrow(final_csv), "species to their BugGuide groups!\n")

Download and process climate data

# ==========================================
# 3. BIOCLIM DATA, CONTINENT MASK & CUSTOM ALBERS
# ==========================================
# Download 30s (~1km) global bioclim data
bioclim_global <- worldclim_global(var = "bio", res = 0.5, path = tempdir())

# 1. Isolate the North American Continent
world_borders <- geodata::world(path = tempdir())
# Subset to the primary NA landmass (add Central American countries here if needed)
na_borders <- world_borders[world_borders$NAME_0 %in% c("Canada", "United States", "Mexico"), ]

cat("Cropping and masking climate data to the North American landmass...\n")
# Crop sets the bounding box; Mask converts everything outside the borders (oceans) to NA
bioclim_na <- crop(bioclim_global, na_borders)
bioclim_na <- mask(bioclim_na, na_borders)
names(bioclim_na) <- paste0("BIO", 1:19)

# 2. Define Custom Albers Equal Area Conic (Central Meridian = -100)
# This centers the map over the Great Plains instead of further east
custom_albers <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-100 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

# 3. Project Layer-by-Layer to bypass GDAL RAM limits
cat("Projecting climate raster to custom Albers layer by layer (this will take a while)...\n")
proj_layers <- list()

for (i in 1:nlyr(bioclim_na)) {
  var_name <- names(bioclim_na)[i]
  cat(sprintf("Processing %s (%d of 19)...\n", var_name, i))
  
  # Project just this single layer and save it to a temp file
  temp_file <- file.path(tempdir(), paste0(var_name, "_albers.tif"))
  
  layer_proj <- project(
    x = bioclim_na[[i]], 
    y = custom_albers, 
    res = 1000, 
    filename = temp_file, 
    overwrite = TRUE
  )
  
  proj_layers[[i]] <- layer_proj
}

# Stack the 19 temporary layers back into a single multi-band raster
bioclim_albers <- rast(proj_layers)
bioclim_albers_trimmed <- trim(bioclim_albers)

## trim of the aleutions and hawaii
# 1. Extract the current bounding box dimensions
current_ext <- ext(bioclim_albers_trimmed)

# 2. Create a new extent, keeping East, South, and North the same, but moving West to -3,400,000
new_ext <- ext(-3400000, current_ext[2], current_ext[3], current_ext[4])

# 3. Crop the raster using the new bounding box
bioclim_albers_final <- crop(bioclim_albers_trimmed, new_ext)

# 4. Plot to verify the new, perfectly tight continental frame
plot(bioclim_albers_final[[1]], maxcell = 5000000, main = "Final Continental Albers (BIO1)")

# 5. Overwrite the main climate raster with this optimized version
writeRaster(bioclim_albers_final, "data/bioclim_na_albers_1km.tif", overwrite = TRUE)

# 6. Create and save a 100km template for richness
template_100km <- aggregate(bioclim_albers_final[[1]], fact = 100)
writeRaster(template_100km, "data/template_100km.tif", overwrite = TRUE)

#  Clean up the individual temp files to free up SSD space
file.remove(unlist(lapply(proj_layers, sources)))


# ==========================================
# 3a. ***PATCH*** TOPOGRAPHIC ROUGHNESS PATCH
# ==========================================
library(terra)
library(geodata)

# 1. Download global elevation data (same 30-arc-second resolution)
cat("Downloading elevation data...\n")
elev_global <- elevation_global(res = 0.5, path = tempdir())

# 2. Crop and mask to the existing North American borders
elev_na <- crop(elev_global, na_borders)
elev_na <- mask(elev_na, na_borders)

# 3. Compute Topographic Complexity 
# TRI (Terrain Ruggedness Index) is the mean difference between a central pixel and its 8 surrounding cells.
cat("Calculating Terrain Ruggedness Index (TRI)...\n")
tri_na <- terrain(elev_na, v = "TRI")
names(tri_na) <- "TRI"

# 4. Project TRI to your Custom Albers
cat("Projecting TRI to Custom Albers...\n")
tri_albers <- project(
  x = tri_na, 
  y = custom_albers, 
  res = 1000, 
  filename = file.path(tempdir(), "tri_albers_temp.tif"),
  overwrite = TRUE
)

# 5. Load your beautifully trimmed 19-layer Bioclim raster
bioclim_albers_final <- rast("data/bioclim_na_albers_1km.tif")

# 6. Crop the newly projected TRI to match your exact western cutoff (-3400000)
tri_albers_matched <- crop(tri_albers, ext(bioclim_albers_final))

# 7. Stack them together into a 20-layer super-raster
env_stack <- c(bioclim_albers_final, tri_albers_matched)

# 8. Overwrite the main file with the new 20-layer version
writeRaster(env_stack, "data/bioclim_na_albers_1km.tif", overwrite = TRUE)
cat("Success! TRI is now Layer 20.\n")

Aggregate occurrence records & extract variables

# ==========================================
# 4. SPATIAL AGGREGATION & EXTRACTION
# ==========================================
# Convert GBIF points to a spatial vector and project them to our new custom Albers
pts <- vect(as.data.frame(df_clean), geom = c("decimalLongitude", "decimalLatitude"), crs = "EPSG:4326")
pts_albers <- project(pts, custom_albers)

# Extract cell IDs from the projected Albers raster
df_clean$cell <- cellFromXY(env_stack, crds(pts_albers))
df_clean <- df_clean %>% filter(!is.na(cell))

# Aggregate occurrences by grid cell and species to reduce sampling bias
df_agg <- df_clean %>%
  group_by(cell, genus, species) %>%
  summarise(occurrence_count = n(), .groups = "drop")

# Retrieve the central projected X/Y coordinates for each populated cell
cell_coords <- xyFromCell(env_stack, df_agg$cell)
df_agg$x <- cell_coords[, "x"]
df_agg$y <- cell_coords[, "y"]

# Extract the 19 Bioclim variables for those specific cells
cat("Extracting climate variables for occurrences...\n")
clim_vals <- terra::extract(env_stack, df_agg$cell)

# Bind the climate data (dropping the redundant 'ID' column from extract)
final_data <- cbind(df_agg, clim_vals)

# ==========================================
# 5. ROBUST SPECIES FILTERING
# ==========================================
# Filter for species with at least 50 occurrences to ensure model stability
species_counts <- final_data %>%
  group_by(species) %>%
  summarize(total_n = sum(occurrence_count)) %>%
  arrange(desc(total_n))

robust_species <- species_counts %>%
  filter(total_n >= 50) %>%
  pull(species)

final_data_pruned <- final_data %>%
  filter(species %in% robust_species)

write.csv(final_data_pruned, "data/melanoplus_aggregated.csv", row.names = FALSE)
cat("Data prep complete. Saved", length(unique(final_data_pruned$species)), "robust Melanoplus species.\n")

Create species groups

# ==========================================
# 6. CREATE SPECIES GROUPS
# ==========================================
library(dplyr)
library(stringr)
library(tidyr)

# 1. Read and parse the TaxonWorks JSON text
cat("Parsing BugGuide/TaxonWorks groups...\n")
lines <- readLines("raw/Melanoplus_groups.txt", warn = FALSE)
raw_text <- paste(lines, collapse = " ")
chunks <- str_split(raw_text, "supersp\\.")[[1]]

mapping <- data.frame(species = character(), species_group = character(), stringsAsFactors = FALSE)

for (i in 2:length(chunks)) {
  chunk <- chunks[i]
  group_name_match <- str_extract(chunk, "^\\s*([a-z]+)")
  
  if (!is.na(group_name_match)) {
    clean_group_name <- paste0(trimws(group_name_match), " group")
    sp_matches <- str_extract_all(chunk, "Melanoplus\\s+[a-z]+")[[1]]
    
    if (length(sp_matches) > 0) {
      chunk_df <- data.frame(
        species = unique(sp_matches),
        species_group = clean_group_name,
        stringsAsFactors = FALSE
      )
      mapping <- rbind(mapping, chunk_df)
    }
  }
}

# Ensure no duplicate species exist in the mapping dictionary
mapping <- distinct(mapping, species, .keep_all = TRUE)

# 2. Load GBIF species list, merge with mapping, and save
cat("Merging groups with GBIF species list...\n")
gbif_species <- read.csv("data/Melanopulus_species.csv")

# Safely rename 'x' to 'species' if it's the raw file
if ("x" %in% colnames(gbif_species)) {
  gbif_species <- dplyr::rename(gbif_species, species = x)
}

final_csv <- gbif_species %>%
  # 1. Force the dataframe to ONLY keep the 'species' column, destroying any old groups
  dplyr::select(species) %>%
  dplyr::distinct(species, .keep_all = TRUE) %>%           
  
  # 2. Now the join has zero risk of .x/.y collisions
  dplyr::left_join(mapping, by = "species") %>%            
  
  # 3. Apply the "Ungrouped" label to anything missing a group
  dplyr::mutate(species_group = ifelse(is.na(species_group), "Ungrouped", species_group)) %>%
  
  # 4. Drop NA species rows
  tidyr::drop_na(species)                           

# 3. Save the clean mapping file for the Shiny App
write.csv(final_csv, "data/Melanopulus_species.csv", row.names = FALSE)

cat("\nSuccessfully mapped", nrow(final_csv), "species.\n")
cat("Ungrouped Status (FALSE = Successfully Grouped):\n")
print(table(final_csv$species_group == "Ungrouped"))