GBIF biodiversity example

This short example shows how the packages rgbif, taxa, and metacoder can be used to retrieve biodiversity data, filter it by taxonomic characteristics, and visualize it. Each of these packages do much more than this; see the links at the end of this example for more information.

Install packages

install.packages("metacoder")
install.packages("taxa")
install.packages("rgbif")

Load packages

library(rgbif)
library(taxa)
library(metacoder)

Geting biodiversity data

We can use the rgbif package to get a table of recorded plant occurrences in Oregon from the Global Biodiversity Information Facility (GBIF). Each row is a record of a plant being observed.

occ <- occ_data(scientificName = "Plantae", stateProvince = "Oregon")
occ$data
## # A tibble: 500 x 72
##       key scientificName decimalLatitude decimalLongitude issues datasetKey publishingOrgKey
##     <int> <chr>                    <dbl>            <dbl> <chr>  <chr>      <chr>           
##  1 2.01e9 Quercus suber…            44.6            -123. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  2 1.99e9 Polystichum m…            45.5            -123. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  3 1.99e9 Chrysolepis c…            44.4            -122. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  4 1.99e9 Lysichiton am…            46.1            -124. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  5 1.99e9 Cardamine hir…            45.5            -123. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  6 1.99e9 Thuja plicata…            45.4            -123. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  7 2.01e9 Quercus garry…            44.9            -123. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  8 1.99e9 Arum italicum…            45.5            -122. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
##  9 1.99e9 Umbellularia …            43.5            -124. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
## 10 1.99e9 Oenanthe sarm…            44.6            -123. cdrou… 50c9509d-… 28eb1a3f-1c15-4…
## # … with 490 more rows, and 65 more variables: networkKeys <list>, installationKey <chr>,
## #   publishingCountry <chr>, protocol <chr>, lastCrawled <chr>, lastParsed <chr>, crawlId <int>,
## #   basisOfRecord <chr>, taxonKey <int>, kingdomKey <int>, …

Converting to the taxmap format

Note how the taxonomy is encoded in a set of columns:

occ$data[26:31]
## # A tibble: 500 x 6
##    kingdom phylum       order        family          genus        species                 
##    <chr>   <chr>        <chr>        <chr>           <chr>        <chr>                   
##  1 Plantae Tracheophyta Fagales      Fagaceae        Quercus      Quercus suber           
##  2 Plantae Tracheophyta Polypodiales Dryopteridaceae Polystichum  Polystichum munitum     
##  3 Plantae Tracheophyta Fagales      Fagaceae        Chrysolepis  Chrysolepis chrysophylla
##  4 Plantae Tracheophyta Alismatales  Araceae         Lysichiton   Lysichiton americanus   
##  5 Plantae Tracheophyta Brassicales  Brassicaceae    Cardamine    Cardamine hirsuta       
##  6 Plantae Tracheophyta Pinales      Cupressaceae    Thuja        Thuja plicata           
##  7 Plantae Tracheophyta Fagales      Fagaceae        Quercus      Quercus garryana        
##  8 Plantae Tracheophyta Alismatales  Araceae         Arum         Arum italicum           
##  9 Plantae Tracheophyta Laurales     Lauraceae       Umbellularia Umbellularia californica
## 10 Plantae Tracheophyta Apiales      Apiaceae        Oenanthe     Oenanthe sarmentosa     
## # … with 490 more rows

We can use these columns to infer a taxonomy and organize the data according to it. The function parse_tax_data from the taxa package can convert data that includes a taxonomy into the taxmap format.

obj <- parse_tax_data(occ$data, class_cols = 26:31, named_by_rank = TRUE)
obj
## <Taxmap>
##   454 taxa: ab. Plantae, ac. Tracheophyta ... rm. Clematis vitalba
##   454 edges: NA->ab, ab->ac, ab->ad, ab->ae, ab->af ... im->rj, gs->rk, fq->rl, ki->rm
##   1 data sets:
##     tax_data:
##       # A tibble: 500 x 73
##         taxon_id    key scientificName decimalLatitude decimalLongitude issues datasetKey
##         <chr>     <int> <chr>                    <dbl>            <dbl> <chr>  <chr>     
##       1 kj       2.01e9 Quercus suber…            44.6            -123. cdrou… 50c9509d-…
##       2 kk       1.99e9 Polystichum m…            45.5            -123. cdrou… 50c9509d-…
##       3 kl       1.99e9 Chrysolepis c…            44.4            -122. cdrou… 50c9509d-…
##       # … with 497 more rows, and 66 more variables: publishingOrgKey <chr>,
##       #   networkKeys <list>, installationKey <chr>, publishingCountry <chr>,
##       #   protocol <chr>, lastCrawled <chr>, lastParsed <chr>, crawlId <int>,
##       #   basisOfRecord <chr>, taxonKey <int>, …
##   0 functions:

Note how the original table is now part of the object and a unique taxon ID has been added to each row.

Plotting and filtering taxonomic data

Using the taxmap format keeps the taxonomic and non-taxonomic parts of the data in sync, allowing for filtering data based on taxon characteristics and vice versa. It also allows you to use the heat_tree function from the metacoder package to plot per-taxon data. In the example below, we plot the number of occurrence for each taxon after filtering out:

  • any taxa at the rank of species
  • any taxa not a subtaxa of Tracheophyta (vascular plants)
  • any order that does not have at least 10 subtaxa
obj %>% 
  filter_taxa(taxon_ranks != "species") %>% 
  filter_taxa(taxon_names == "Tracheophyta", subtaxa = TRUE) %>%
  filter_taxa(taxon_ranks == "order", n_subtaxa > 10, subtaxa = TRUE, supertaxa = TRUE) %>% 
  heat_tree(node_label = taxon_names,
            node_color = n_obs,
            node_size = n_obs,
            node_color_axis_label = "# occurrences")