Natural history records and the case of the missing absence.
Target-group background sampling for CART species distribution models
Overview: The Pseudo-Absence Problem
Global databases like GBIF provide “presence-only” data. A botanist records where they found a plant present, but they rarely upload coordinates for where the plant was absent. However, machine learning algorithms like Classification and Regression Trees (CART) require both presences and absences to learn the boundaries of a species’ niche.
Pseudo-absences are often drawn at random across the map in hopes of catching the broad pattern. This introduces errors that complicate the interpretation of the model. For example, random points might land in the middle of a military testing range where no botanist has ever looked. Our model might wrongly conclude the plant can’t live there due to climate, when in reality, it’s just unsurveyed.
The Target-Group Background (TGB) Solution
Luckily botanists do often note other plants growing nearby. We can use this attribute of occurrence records as a biological basis for selecting absences. To solve this, we define our sample domain for absences using associated taxa. For an area with many occurrence records for species that commonly grows alongside our target species, but didn’t record our target, we can have some confidence that our target is likely to absent there.
Distributions are rarely in equilibrium with the climate, so further consideration about absences is needed here. For example, some plant distrubtions respond quickly to interglacial periods while others are slow and may still be adjusting to the current inter-glacial climate. A
In this tutorial we will:
1. Download occurrences for Brittlebush (Encelia farinosa).
2. Extract and summarize the top 10 associated taxa.
3. Define a 100km sample domain buffer (simplified to a WKT envelope).
4. Download the associated taxa within that domain to serve as our “Background.”
5. Spatially isolate the pixels that contain background species, but lack E. farinosa.
6. Run a CART SDM using these true presences and validated absences.
1. Target Species Acquisition
First, we acquire the occurrence records for our single target species and convert them to a spatial vector.
# ==============================================================================
# Script 1: target_species.r
# Purpose: Download and clean Encelia farinosa occurrences
# ==============================================================================
library(terra)
library(rgbif)
library(dplyr)
library(tidyr)
library(stringr)
# 1.1 Download E. farinosa
target_key \<- name_backbone("Encelia farinosa")\$usageKey
gbif_target \<- occ_download(
pred("taxonKey", target_key),
pred("hasCoordinate", TRUE),
pred("hasGeospatialIssue", FALSE),
pred_gte("year", 2000),
format = "SIMPLE_CSV"
)
occ_download_wait(gbif_target)
target_data \<- occ_download_import(occ_download_get(gbif_target, path = "data"))
# 1.2 Clean target data
target_clean \<- target_data %\>%
filter(!is.na(decimalLongitude), !is.na(decimalLatitude))
# Convert to spatial vector
target_vect \<- vect(target_clean, geom = c("decimalLongitude", "decimalLatitude"), crs = "EPSG:4326")2. Summarizing Associated Taxa
GBIF data utilizes the Darwin Core format, which often includes a column named associatedTaxa. We will parse this column, count the species, and identify the top 10 most common companions to E. farinosa.
# ==============================================================================
# Script 2: associated_taxa.r
# Purpose: Extract and rank the most common co-occurring species
# ==============================================================================
# 2.1 Parse the associatedTaxa column
# The column is often a messy text string separated by pipes (|) or commas.
# We will separate the strings, clean the whitespace, and tally them up.
assoc_summary <- target_clean %>%
filter(!is.na(associatedTaxa), associatedTaxa != "") %>%
separate_rows(associatedTaxa, sep = "[|,;]") %>% # split by common delimiters
mutate(associatedTaxa = str_squish(associatedTaxa)) %>% # remove extra whitespace
count(associatedTaxa, sort = TRUE) %>%
filter(n > 5) # Filter out rare/typo taxa
# 2.2 View the top 10
head(assoc_summary, 10)
# Save summary to table
write.csv(assoc_summary, "outputs/associated_taxa_summary.csv", row.names = FALSE)
# 2.3 Store the names of the top 10 (excluding generic/unidentified terms if any)
# Note: You may need to manually verify these names are clean binomials
top_10_names <- head(assoc_summary$associatedTaxa, 10)3. Defining the Sample Domain (100km Buffer)
We only want to draw absences from the general geographic region where our target species exists. We will buffer our target presences by 100km, and then simplify that shape into a single bounding box (extent) so we can easily pass it to GBIF as WKT.
# ==============================================================================
# Script 3: sample_domain.r
# Purpose: Create a 100km geographic domain around known presences
# ==============================================================================
# Buffer the occurrences by 100km (100,000 meters)
# Note: buffering in geographic coordinates (WGS84) requires defining meters
domain_buffer <- buffer(target_vect, width = 100000)
# To simplify the GBIF request, we convert this complex buffer into a simple bounding box
domain_extent <- ext(domain_buffer)
# Convert the bounding box extent into a spatial polygon vector
domain_poly <- as.polygons(domain_extent, crs = "EPSG:4326")
# Convert to WKT format for the GBIF API
domain_wkt <- geom(domain_poly, wkt = TRUE)4. Downloading the Missing Absences
Now we query GBIF for any observations of our top 10 associated taxa that fall within our 100km WKT bounding box.
# ==============================================================================
# Script 4: acquire_background.r
# Purpose: Download occurrences of associated taxa within the sample domain
# ==============================================================================
# Get taxon keys for our top 10 associated species
# We use purrr/lapply to get keys for the whole list
assoc_keys <- sapply(top_10_names, function(x) name_backbone(x)$usageKey)
# Download the background data
gbif_bg <- occ_download(
pred_in("taxonKey", assoc_keys),
pred("hasCoordinate", TRUE),
pred("hasGeospatialIssue", FALSE),
pred_within(domain_wkt),
format = "SIMPLE_CSV"
)
occ_download_wait(gbif_bg)
bg_data <- occ_download_import(occ_download_get(gbif_bg, path = "data"))
# Clean and convert to spatial vector
bg_clean <- bg_data %>% filter(!is.na(decimalLongitude), !is.na(decimalLatitude))
bg_vect <- vect(bg_clean, geom = c("decimalLongitude", "decimalLatitude"), crs = "EPSG:4326")5. Identifying the Final Absence Pixels
We now have presences and background data. We will create a grid. Any pixel containing E. farinosa is a 1 (Presence). Any pixel containing a background species, but no E. farinosa, becomes a 0 (Absence).
# ==============================================================================
# Script 5: pixel_classification.r
# Purpose: Isolate pixels with associated species but lacking the target species
# ==============================================================================
# 5.1 Create a template grid at ~10km resolution covering our domain
grid_template <- rast(domain_poly, resolution = 0.1, crs = "EPSG:4326")
# 5.2 Rasterize Presences (Assign value of 1)
# length function counts how many points fall in the cell. We convert any count > 0 to 1.
pres_raster <- rasterize(target_vect, grid_template, fun = length)
pres_raster <- ifel(pres_raster > 0, 1, NA)
# 5.3 Rasterize Background (Assign value of 0)
bg_raster <- rasterize(bg_vect, grid_template, fun = length)
bg_raster <- ifel(bg_raster > 0, 0, NA)
# 5.4 Combine and classify
# We use terra::cover. It takes values from the first raster (presences).
# If a pixel is NA in the first, it takes the value from the second (absences).
sdm_grid <- cover(pres_raster, bg_raster)
# Let's extract these pixel coordinates to build our final modeling dataframe
sdm_points <- as.points(sdm_grid)
sdm_df <- geom(sdm_points)[, c("x", "y")] %>% as.data.frame()
sdm_df$occurrence <- sdm_points[[1]] # Extract the 1s and 0s
cat("Total Presences (1):", sum(sdm_df$occurrence == 1), "\n")
cat("Total Validated Absences (0):", sum(sdm_df$occurrence == 0), "\n")6. CART Model Development
With our presences and validated absences ready, we can extract environmental variables and run our decision tree.
# ==============================================================================
# Script 6: cart_sdm.r
# Purpose: Extract climate data and build the CART model
# ==============================================================================
library(geodata)
library(rpart)
library(rpart.plot)
# 6.1 Get climate data and extract to our points
bioclim <- worldclim_global(var = "bio", res = 5, path = "data")
sdm_vect <- vect(sdm_df, geom = c("x", "y"), crs = "EPSG:4326")
clim_extract <- extract(bioclim, sdm_vect, ID = FALSE)
# Combine occurrence status (1/0) with climate predictors
model_data <- data.frame(occurrence = as.factor(sdm_df$occurrence), clim_extract)
model_data <- na.omit(model_data)
# 6.2 Build and Plot the CART Model
# We use all predictors (.) to see which ones the tree selects
cart_model <- rpart(occurrence ~ ., data = model_data, method = "class")
# Save the decision tree plot
tiff("outputs/encelia_farinosa_sdm_tree.tif", width = 2400, height = 1800, res = 300)
rpart.plot(cart_model,
type = 2,
extra = 104,
box.palette = "RdYlGn", # Red for absence, Green for presence
main = "SDM Decision Tree: E. farinosa")
dev.off()
# 6.3 Spatial Prediction
# Predict across the climate raster to see the potential distribution
predicted_map <- predict(bioclim, cart_model, type = "class")
# Crop map to our sample domain for plotting
predicted_domain <- crop(predicted_map, domain_poly)
predicted_domain <- mask(predicted_domain, domain_poly)
plot(predicted_domain, main = "Predicted E. farinosa Distribution (1 = Presence, 0 = Absence)")
writeRaster(predicted_domain, "outputs/e_farinosa_predicted_map.tif", overwrite = TRUE)