Natural History Occurence records in R

a simple pipeline to acquire and prepare records for spatial analysis

Author

J Gibson, Arizona Western College

Published

September 14, 2026

Overview of GBIF occurrence data and R

Digitized Natural History records are basically the best development of the 21st Century. You have access to over 3.75 billion occurrence records across the planet, and this grows by about 420 million records a year. You can contribute to this data set and perform cutting edge analysis on it. R allows access to this data and combines it with the data management, statistical and spatial capabilities of the R program. With this, every student can contribute to the radical transformation of biodiversity science, sensu Sutherland et al. 2026

Global Biodiversity Information Facility

GBIF Aggregates natural history collections from data portals, including physical collections like those published through the Big Bee Library or SCAN & human observations like those published through iNaturalist.

The R Project for Statistical Computing

The R program provides a reliable, repeatable, and shareable way to access GIBF. R is a statistical programming language that also handles spatial data.

R studio

The RStudio environment provides a standardized coding environment across various operating systems. The RStudio environment combines a script editor with an R console. A script editor allows you to develop R code and save it. The R console allows you to execute the code on input and output datasets.

Establishing a stable code pipeline from GBIF to R

A sequence of processing is called a pipeline. We need a stable pipeline for accessing and preparing occurrence records. This tutorial walks through a simple pipeline you can follow along and adapt to your own taxonomic and geographic interest. You will need to have R studio.

Basic R tutorials

Look to the source: https://education.rstudio.com/learn/beginner/

The process pipeline

The process has four main components. In each component we will develop a script.

1. setup_project.r

  • 1.1 define directories

  • 1.2 install packages

  • 1.3 set R environment

2. search_parameters.r

  • 2.1 geographic parameters

  • 2.2 taxonomic parameters

  • 2.3 other parameters

3. download_clean.r

  • 3.1 download data

  • 3.2 clean the occurrence data

  • 3.3 summarize occurrence records as a table

4. occurrence_maps.r

  • 4.1 make occurrence records spatial

  • 4.2 get base map

  • 4.3 create occurrence map as a figure

  • 4.4 save map as a .tif

1. Project setup (project_setup.r)

Purpose: This script will define the directory structure for saving the data that will be downloaded and processed, acquire necessary packages, and set the R environment for connecting to GBIF.

Concepts:

  • R can read and write directory structures (e.g. folders) as well as files (e.g. occurrence data).

  • Begin each script with a header.

  • Objects can be used as arguments in a function.

Open RStudio and open a new R Script

  • In Rstudio select File -> New File -> R Script

Begin the script with a header that has your name and the purpose of the script. Copy, paste, and adapt this code at the top of your blank new script.

# ==============================================================================
# Script: project_setup.r
# Author: Your Name
# Purpose: Setup a project directory structure and acquire packages
# ==============================================================================
#
#

1.1 Define Directories

  • Next we define the directory structure.

    • create objects of the text string of the directory path

    • use those objects as the argument for the dir.create() function

  • Add the following code to your script

    • adapt it to the root directory where you want to create your project folder.

    • NOTE: make sure to use forward slashes / in your pathnames, not back slashes

# ==============================================================================
# 1.1 Define directory structure, check if it exists and create if not
# ==============================================================================
## Define directory structure
root_dir <- "C:/example" #NOTE: make sure you seperate folders with the forward slash / 
setwd(root_dir) 
project_dir <- paste(root_dir,"project",sep="/") # change the word project to the name you want your folder to be
occurrence_dir <- paste(project_dir,"occurrence",sep="/") # folder for tabular occurence data 
maps_dir <- paste(project_dir,"maps",sep="/") # folder for spatial data 
scripts_dir <- paste(project_dir,"scripts",sep="/") # folder for project scripts
  • Now run code to create the the directories if they do not exist.
## create directories. This will fail if the file already exists.

dir.create(project_dir)
dir.create(occurrence_dir)
dir.create(maps_dir)
dir.create(scripts_dir)
  • Save your script now.
    • In Rstudio select File -> Save As…
    • Navigate to your new script folder and save as projectSetup.r

1.2 Install R Packages

We will need a set of packages with functions for working with GBIF and spatial data.

  • Add the following code to your script and run. You do not need to modify this code.
# ==============================================================================
# 1.2 Check and install package dependencies
# ==============================================================================
required_packages <- c("rgbif", "dplyr", "purrr", "CoordinateCleaner", "geodata", "terra","maptiles")
missing_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]

if(length(missing_packages)) {
  install.packages(missing_packages)
  message("Missing packages installed successfully.")
} else {
  message("All packages are already installed!")
}

  • You will be prompted to select a CRAN. This is where you will download from.

  • Save your script! RStudio does not automatically save your work.

1.3 Set R environment for GBIF API permissions

Set your GBIF credentials in the R environment.

  • If you have not already, make an account with GBIF. This is necassary for the data_download() function in the rgbif library.

  • Add, adapt, and run the following code:

# ==============================================================================
# 1.3 set R environment for GBIF API permissions
# ==============================================================================
## set up gbif creds
#install.packages("usethis")
usethis::edit_r_environ()

GBIF_USER="username"
GBIF_PWD="password"
GBIF_EMAIL="email"
  • This is the end of the script. Make sure to save it.

2. Search Parameters (search_parameters.r)

Purpose: This script will define search parameters that will be used to query GBIF. We will set geographic and taxonomic parameters, and consider other common search parameters.

Concepts:

  • GBIF is accessed through the library rgbif and uses the same search terms that are used in the GBIF GUI.

  • Load packages with library() to use them in an R session.

  • The terra package handles spatial data

  • The vect() function creates or reads spatial data in vector format, such as .shp or .kml

  • The rast() function creates or reads spatial data in raster format, such as .tif or .grd

  • The help() function opens the documentation for libraries.

In RStudio select File -> New File -> R Script. Save as study_area.r in your scripts folder. Add and adapt the following code to the new script.

# ==============================================================================
# Script: search_parameters.r
# Author: Your Name
# Purpose: Define the gegraphic and taxonomic scope of your work
# ==============================================================================

## load libraries
library(terra)
library(rgbif)

## directory structure
#copy your directory structure from the previous script and paste here
#
#

2.1 Geographic parameters

Geography can be used as a search parameter in two main ways.

  • Define a geometry like a circular rectangular area.

    • In this case you pass the geometry as an argument.

  • Use one of the existing geography fields in the GBIF database, such as administrative unit or continent.

    • In this case you pass the search term as an argument.

In this example we define a study area as a circular area around a point.

# ==============================================================================
# 2.1 Define a Study Area 
# ==============================================================================
# Point coordinates for center of circle
latitude <- 32.69
longitude <- -114.62

# Define radii of circle in meters

radii <- 100000

# Format your coordinates into a matrix
point_coords <- matrix(c(longitude, latitude), ncol = 2)

# use coordinates to create a spatial vector with WGS84 projection
pt <- vect(point_coords, crs = "EPSG:4326")

# use the buffer function to create a circle
circle <- buffer(pt, width = radii)

# Check out the vector by submitting the name of the vector object
circle

# Plot the vector
plot(circle)

# Write out the shapefile in the ESRI shp format
output_shp <- "circle.shp" # name of output .shp file

writeVector(
  x         = circle, 
  filename  = paste(maps_dir, output_shp, sep="/"),
  filetype  = "ESRI Shapefile", 
  overwrite = TRUE
)

# Check to see output files were written
list.files(maps_dir)

2.2 Taxonomic paramaters

Before downloading data, we need to find how GBIF classifies our target taxon. We use the name_backbone() function to look up the unique taxonKey. This can be done at many taxonomic levels, but not all levels. Families can be searched, but not epifamilies, for example.

# ==============================================================================
# 2.2 Taxonomic Parameters
# ==============================================================================
# You can search for a taxon, like all bees in the family Andrenidae
target_taxon <- name_backbone("Andrenidae")

# Extract the unique usageKey 
taxon_key <- target_taxon$usageKey

cat("The GBIF Taxon Key for your taxon is:", taxon_key, "\n")

# Or you can search for multiple taxa, like all bee families
bee_families <- c("Apidae","Halictidae","Andrenidae","Megachilidae","Colletidae","Siricidae")

bee_keys <- sapply(bee_families, function(x) name_backbone(x)$usageKey)

cat("The GBIF Taxon Keys for your taxa are:", bee_keys, "\n")

2.3 Other parameters

There are many other ways to refine your occurrence search prior to downloading. You can refer to the official documentation by running help(download_predicate_dsl) in your console. Common filters include:

  • Administrative Units: You can use specific codes such as ISO-2 letter country codes (e.g., pred(“country”, “US”)) or GADM codes.

  • Continents: Broad geographic regions (e.g., pred(“continent”, “NORTH_AMERICA”)).

  • Basis of Observation: Filter by how the data was gathered (e.g., pred(“basisOfRecord”, “HUMAN_OBSERVATION”) or “PRESERVED_SPECIMEN”).

  • Time of Observation: Target specific times (e.g., pred(“year”, “2020”), pred_gte(“year”, 2000), or pred(“month”, “4”)).

  • Publishing Institution / Observer: Target specific institutions (publishingOrg) or the person who recorded it (recordedBy).

# ==============================================================================
# 2.3 Other Parameters
# ==============================================================================
# Example of defining a minimum year to exclude very old historical data
min_year <- 2000

3. Download Data (download_clean.r)

Purpose: This script will downloaded GBIF data and clean it for further analysis.

Concepts:

  • The rgbif library provides functions for downloading data.

  • The CoordinateCleaner library provides functions for cleaning common errors from the coordinates of natural history occurrence records.

  • The dplyr library provides functions for summarizing tabular data.

In RStudio select File -> New File -> R Script. Save as download_clean.r in your scripts folder.

# ==============================================================================
# Script: download_clean.r
# Author: Your Name
# Purpose: Download, clean, and summarize GBIF records
# ==============================================================================
library(terra)
library(rgbif)
library(CoordinateCleaner)
library(dplyr)

# (Copy your directory structure from your project_setup.r script and paste here)
# (Copy your taxon_key and any other search parameters from search_parameters.r)

# Read in study area shapefile
circle <- vect(paste(maps_dir, "circle.shp", sep="/"))

3.1 Download data

The occ_download() function allows larger downloads and is recommended for robust data pipelines. This function submits a job to the GBIF servers. We can use the terra package to assign our circle shape to a format that GBIF will accept.

# ==============================================================================
# 3.1 Download GBIF Data
# ==============================================================================
# Convert terra SpatVector to GBIF-Compliant WKT Polygon
# rgbif requires WKT formatted counter-clockwise (anti-clockwise) for exterior rings
# WKT polygon string:
# Extract the coordinate matrix from the terra circle
coords <- terra::geom(circle)

# Reverse the row order to change from clockwise to counter-clockwise
coords_rev <- coords[nrow(coords):1, ]

# Manually format the reversed coordinates into a GBIF-approved WKT string
wkt_poly <- paste0("POLYGON((", paste(coords_rev[, "x"], coords_rev[, "y"], sep = " ", collapse = ", "), "))")

# Note: You must set environment variables GBIF_USER, GBIF_PWD, and GBIF_EMAIL
# or pass them directly to occ_download()

gbif_job <- rgbif::occ_download(
  pred_in("taxonKey", bee_keys), # from search_parameters
  pred("hasCoordinate", TRUE),
  pred("hasGeospatialIssue", FALSE),
  pred_within(wkt_poly),      # Geometry constraint: inside the buffer polygon!
  format = "SIMPLE_CSV"
)

cat("GBIF Download Requested! Download Key:", gbif_job[1], "\n")


# Wait for completion (this will pause R until GBIF finishes processing)
occ_download_wait(gbif_job)

# Fetch zip archive and import as a data frame
res_file <- occ_download_get(gbif_job, path = occurrence_dir)
df_occ   <- occ_download_import(res_file)

3.2 Clean the occurrence data

Raw database occurrences often have errors with coordinates recorded as country centroids, institutions, or in the ocean. We will remove these obvious errors.

# ==============================================================================
# 3.2 Clean the data
# ==============================================================================
# Flag common spatial errors using CoordinateCleaner
clean_occ <- clean_coordinates(
  x = df_occ,
  lon = "decimalLongitude",
  lat = "decimalLatitude",
  species = "species",
  tests = c("centroids", "equal", "gbif", "institutions", "zeros", "seas")
)

# Filter out the flagged records
df_cleaned <- clean_occ %>% 
  filter(.summary == TRUE)
  
# Save the cleaned tabular data
write.csv(df_cleaned, paste(occurrence_dir, "cleaned_occurrences.csv", sep="/"), row.names = FALSE)

3.3 Summarize occurrence records

We can summarize our data in several ways. We can simply look at the number of records in a certain attribute, like family or genus. We can also use multiple attributes to summarize our records, like year of the sample, the basis of observation and genera.

# ==============================================================================
# 3.3 Summarize Records
# ==============================================================================
# Summarize records by a single variable
table(df_cleaned$family)

# Count records by a combination of certain attributes, like genus, year, and basis of record
summary_table <- df_cleaned %>%
  group_by(genus, year, basisOfRecord) %>%
  summarize(record_count = n(), .groups = 'drop')

# Count records by a combination of certain attributes, like genus, year, and basis of record
summary_table <- df_cleaned %>%
  group_by(genus, year, basisOfRecord) %>%
  summarize(record_count = n(), .groups = 'drop')

print(summary_table)

4. Map occurrences (occurrence_maps.r)

Purpose: Turn tabular data into spatial data and visualize it.

Concepts:

  • The maptiles library provides access to many base layers like Open Street Maps. These are for display purposes to make figures and do not contain the actual data.

  • The geodata library provides access to many environmental datasets in gridded format. These are for analysis purposes but can be plotted with symbology to create figures.

In RStudio select File -> New File -> R Script. Save as occurrence_maps.r in your scripts folder.

# ==============================================================================
# Script: occurrence_maps.r
# Author: Your Name
# Purpose: Map cleaned occurrence records
# ==============================================================================
library(terra)
library(maptiles)

# (Copy your directory structure here)

4.1 Make occurrence records spatial

Convert the cleaned dataframe into a SpatVector object so terra can map it.

# ==============================================================================
# 4.1 Create Spatial Vector from DataFrame
# ==============================================================================
# Load the cleaned csv
df_cleaned <- read.csv(paste(occurrence_dir, "cleaned_occurrences.csv", sep="/"))
circle <- vect(paste(maps_dir, "circle.shp", sep="/"))
pt <- centroids(circle) # recreate center point

# Convert to terra SpatVector
occ_vect <- vect(df_cleaned, geom = c("decimalLongitude", "decimalLatitude"), crs = "EPSG:4326")

4.2 Get base map

Next, we acquire a base map to plot our study area on using the maptiles library.

# ==============================================================================
# 4.2 Acquire a base map
# ==============================================================================
# Download OpenStreetMap tiles cropped to our study area
osm <- get_tiles(circle, 
  provider = "OpenStreetMap", 
  zoom     = 8,   
  crop     = TRUE
)

4.3 Create occurrence map as a figure

Now we plot the basemap, the study area boundaries, and occurrence records.

# ==============================================================================
# 4.3 Plot Basemap with Study Area and Occurrences
# ==============================================================================
# Plot the RGB OpenStreetMap tiles
plotRGB(osm, main = "Bee occurrences in Study Area")

# Overlay the study area circle 
plot(circle, add = TRUE, border = "blue", lwd = 2)

# Overlay center point
plot(pt, add = TRUE, col = "red", pch = 19, cex = 1.5)

# Overlay the cleaned occurrences
plot(occ_vect, add = TRUE, col = "darkgreen", pch = 20, cex = 0.8)

4.4 Save map as a .tif We can write our mapped plot directly out of R into an image file format.

# ==============================================================================
# 4.4 Save the Map as a TIF
# ==============================================================================
# Open a tiff graphics device 
tiff(paste(maps_dir, "occurrence_map.tif", sep="/"), width = 800, height = 800, res = 150)

# Re-run the plotting commands
plotRGB(osm, main = "Andrenidae occurrences in Study Area")
plot(circle, add = TRUE, border = "blue", lwd = 2)
plot(pt, add = TRUE, col = "red", pch = 19, cex = 1.5)
plot(occ_vect, add = TRUE, col = "darkgreen", pch = 20, cex = 0.8)

# Close the device to finalize the save
dev.off()