1. Introduction

This document is intended to act as a primer for the use of the new Neotoma R package, neotoma2, for working with diatom data. Some users may be working with this document as part of a workshop for which there is a Binder instance. The Binder instance will run RStudio in your browser, with all the required packages installed.

If you are using this workflow on its own, or want to use the package directly, the neotoma2 package is available on CRAN. This workflow uses the development version (1.1.1), which you can install by running:

install.packages("neotoma2")
library(neotoma2)

This workshop will also require other packages. To maintain the flow of this document we’ve placed instructions at the end of the document in the section labelled “Installing packages on your own”.

2. Learning Goals

In this tutorial you will learn how to:

  • Search for sites using site names and geographic parameters – Site Searches
  • Filter results by dataset type (diatoms, diatom surface samples, water chemistry) and by whether a record has a chronology – Filter Results
  • Obtain sample information for one of the selected sites – Sample Information
  • Perform basic analysis including connecting your taxon names to GBIF and a stratigraphic diagram – Basic Analysis

2.1. Getting Help with Neotoma

If you’re planning on working with Neotoma, please join us on Slack where we manage a channel specifically for questions about the R package. You may also wish to join our Google Groups mailing list, please contact us to be added.

2.2 Understanding Data Structures in Neotoma

Data in the Neotoma database itself is structured as a set of linked relationships to express the different elements of paleoecological analysis: space and time, raw data, scientific methods and data models. Because paleoecology is such a broad discipline these relationships can be complex, and as such, the database itself is highly structured. If you want to better understand concepts within the database, you can read the Neotoma Database Manual, or take a look at the database schema itself.

In this workshop we want to highlight two key structural concepts:

  1. The way data is structured conceptually within Neotoma (Sites, Collection Units and Datasets).
  2. The way that this structure is adapted within the neotoma2 R package.

2.2.1 Data Structure in the Neotoma Database

Figure. The structure of sites, collection units and datasets within Neotoma. A site contains one or more collection units. Chronologies are associated with collection units. Data of a common type (pollen, diatoms, vertebrate fauna) are assigned to a dataset.
Figure. The structure of sites, collection units and datasets within Neotoma. A site contains one or more collection units. Chronologies are associated with collection units. Data of a common type (pollen, diatoms, vertebrate fauna) are assigned to a dataset.

Data in Neotoma is associated with sites – specific locations with lat/long coordinates.

Within a site, there may be one or more collection units – locations at which samples are physically collected within the site. For example, an archaeological site may have one or more collection units, pits within a broader dig site; a diatom sampling site on a lake may have multiple collection units – core sites within the lake basin. Collection units may have higher resolution GPS locations than the site location, but are considered to be part of the broader site.

Data within a collection unit is collected at various analysis units. Any data sampled within an analysis unit is grouped by the dataset type (charcoal, diatom, dinoflagellate, etc.) and aggregated into a sample. The set of samples for a collection unit of a particular dataset type is then assigned to a dataset.

2.2.2 Data Structures in neotoma2

Figure. Neotoma R Package diagram. Each box represents a data class within the package. Individual boxes show the class object, its name, its properties, and functions that can be applied to those objects. For example, a sites object has a property sites, that is a list. The function plotLeaflet() can be used on a sites object.
Figure. Neotoma R Package diagram. Each box represents a data class within the package. Individual boxes show the class object, its name, its properties, and functions that can be applied to those objects. For example, a sites object has a property sites, that is a list. The function plotLeaflet() can be used on a sites object.

If we look at the UML diagram for the objects in the neotoma2 R package we can see that the data structure generally mimics the structure within the database itself. As we will see in the Site Searches section, we can search for these objects, and begin to manipulate them (in the Simple Analysis section).

It is important to note: within the neotoma2 R package, most objects are sites objects, they just contain more or less data. There are a set of functions that can operate on sites. As we add to sites objects, using get_datasets() or get_downloads(), we are able to use more of these helper functions.

2.2.3 Piping in R

Piping is a technique that simplifies the process of chaining multiple operations on a data object. It involves using either of these operators: |> or %>%. |> is a base R operator while %>% comes from the tidyverse ecosystem in R. In neotoma2 we use %>%.

The pipe operator works as a real-life pipe, which carries water from one location to another. In programming, the output of the function on the left-hand side of the pipe is taken as the initial argument for the function on the right-hand side of the pipe. It helps by making code easier to write and read. Additionally, it reduces the number of intermediate objects created during data processing, which can make code more memory-efficient and faster.

Without using pipes you can use the neotoma2 R package to retrieve a site and then plot it by doing:

# Retrieve the site
plot_site <- neotoma2::get_sites(sitename = "%Swamp%")
# Plot the site
neotoma2::plotLeaflet(object = plot_site)

This would create a variable plot_site that we will not need any more, but it was necessary so that we could pass it to the plotLeaflet function.

With the pipe (%>%) we do not need to create the variable, we can just rewrite our code. Notice that plotLeaflet() doesn’t need the object argument because the response of get_sites(sitename = "%Swamp%") gets passed directly into the function.

2.2.3.1. Code

# get_sites and pipe. The `object` parameter for plotLeaflet will be the
# result of the `get_sites()` function.
get_sites(sitename = "%Swamp%") %>%
  plotLeaflet()

2.2.3.2. Result

3. Site Searches

3.1. get_sites()

There are several ways to find sites in neotoma2, but we think of sites as being spatial objects primarily. They have names, locations, and are found within the context of geopolitical units, but within the API and the package, the site itself does not have associated information about taxa, dataset types or ages. It is simply the container into which we add that information. So, when we search for sites we can search by:

Parameter Description
sitename A valid site name (case insensitive) using % as a wildcard.
siteid A unique numeric site id from the Neotoma Database
loc A bounding box vector, geoJSON or WKT string.
altmin Lower altitude bound for sites.
altmax Upper altitude bound for site locations.
database The constituent database from which the records are pulled.
datasettype The kind of dataset (see get_table("datasettypes"))
datasetid Unique numeric dataset identifier in Neotoma
doi A valid dataset DOI in Neotoma
gpid A unique numeric identifier, or text string identifying a geopolitical unit in Neotoma
keywords Unique sample keywords for records in Neotoma.
contacts A name or numeric id for individuals associuated with sites.
taxa Unique numeric identifiers or taxon names associated with sites.

All sites in Neotoma contain one or more datasets. It’s worth noting that the results of these search parameters may be slightly unexpected. For example, searching for sites by sitename, latitude, or altitude will return all of the datasets for the particular site. Searching for terms such as datasettype, datasetid or taxa will return the site, but the only datasets returned will be those matching the dataset-specific search terms. We’ll see this later.

3.1.1. Site names: sitename="%Lake%"

We may know exactly what site we’re looking for (“Arbutus Lake”), or have an approximate guess for the site name (for example, we know it’s something like “Arbutus Lake” or “Lake Arbutus”, but we’re not sure how it was entered specifically), or we may want to search all sites that have a specific term, for example, Lake.

We use the general format: get_sites(sitename="%Lake%") for searching by name.

PostgreSQL (and the API) uses the percent sign as a wildcard. So "%Lake%" would pick up “Arbutus Lake” for us, and it would also pick up any site where “lake” is part of a longer word, such as “Lakeside” or “Blake”. Note that the search query is also case insensitive, so you could simply write "%LAKE%".

3.1.1.1. Code
lake_sites <- neotoma2::get_sites(sitename = "%Lake%")
plotLeaflet(lake_sites)
3.1.1.2. Result

3.1.2. Location: loc=c()

The neotoma package used a bounding box for locations, structured as a vector of latitude and longitude values: c(xmin, ymin, xmax, ymax). The neotoma2 R package supports both this simple bounding box, but also more complex spatial objects, using the sf package. Using the sf package allows us to more easily work with raster and polygon data in R, and to select sites from more complex spatial objects. The loc parameter works with the simple vector, WKT, geoJSON objects and native sf objects in R.

Here we look for sites using a location. We’re using a rough polygon around Pennsylvania, stretched east to include the New Jersey shore out to Barnegat Bay. To work with this spatial object in R we also transformed the geojson_str element to an object for the sf package. There are many other tools to work with spatial objects in R. Regardless of how you get the data into R, neotoma2 works with almost all objects in the sf package.

geojson_str <- '{"coordinates":
  [[[-80.52, 42.27],
    [-79.76, 42.27],
    [-75.36, 42.00],
    [-74.69, 41.36],
    [-73.95, 40.60],
    [-73.95, 39.45],
    [-75.55, 39.45],
    [-75.79, 39.72],
    [-80.52, 39.72],
    [-80.52, 42.27]]],
  "type":"Polygon"}'

penn_sf <- geojsonsf::geojson_sf(geojson_str)

# Note here we use the `all_data` flag to capture all the sites.
penn_sites <- neotoma2::get_sites(loc = penn_sf, all_data = TRUE)

You can always simply plot() the sites objects, but you will lose some of the geographic context. The plotLeaflet() function returns a leaflet() map, and allows you to further customize it, or add additional spatial data (like our original bounding polygon, penn_sf, which works directly with the R leaflet package):

3.1.2.1. Code
neotoma2::plotLeaflet(penn_sites) %>%
  leaflet::addPolygons(
    map = .,
    data = penn_sf,
    color = "green"
  )
3.1.2.2. Result

3.1.3. Site Helpers

If we look at the data structure diagram for the objects in the neotoma2 R package we can see that there are a set of functions that can operate on sites. As we add to sites objects, using get_datasets() or get_downloads(), we are able to use more of these helper functions.

As it is, we can take advantage of functions like summary() to get a more complete sense of the types of data we have in penn_sites. The following code gives the summary table. We do some R magic here to change the way the data is displayed (turning it into a datatable() object), but the main piece is the summary() call.

3.1.3.1. Code
# Give information about the sites themselves, site names &cetera.
neotoma2::summary(penn_sites)
# Give the unique identifiers for sites, collection units and datasets found at
# those sites.
neotoma2::getids(penn_sites)
3.1.3.2. Result

In this document we list only the first 10 records (there are more, you can use length(penn_sites) to see how many sites you’ve got). The dataset_types column shows what each site holds: mostly diatoms and water chemistry, but also pollen, geochronologic data and other types. We can see that there are no chronologies associated with the site objects (n_chronologies is 0). This is because, at present, we have not pulled in the dataset information we need. In Neotoma, a chronology is associated with a collection unit (and that metadata is pulled by get_datasets() or get_downloads()). All we know from get_sites() are the kinds of datasets we have and the location of the sites that contain the datasets.

3.2. Searching for datasets:

We know that within Neotoma, collection units and datasets are contained within sites. Similarly, a sites object contains collectionunits which contain datasets. From the table above we can see that the sites we’ve looked at contain a mix of dataset types. We could write something like this: table(unlist(strsplit(summary(penn_sites)$dataset_types, ","))) to see the different datasettypes and their counts.

With a sites object we can directly call get_datasets(), to pull in more metadata about the datasets. The get_datasets() method also supports any of the search terms listed above in the Site Search section. At any time we can use datasets() to get more information about any datasets that a sites object may contain. Compare the output of datasets(penn_sites) to the output of a similar call using the following:

3.2.1. Code

penn_datasets <- neotoma2::get_datasets(penn_sites)

datasets(penn_datasets)

# How many datasets of each type are there?
table(datasets(penn_datasets)$datasettype)

3.2.2. Result

You can see that this provides information only about the specific dataset, not the site! For a more complete record we can join site information from summary() to dataset information using datasets() using the getids() function which links sites, and all the collection units and datasets they contain.

3.3. Filter Records

If we choose to pull in information about only certain dataset types, or if there is additional filtering we want to do before we download the data, we can use the filter() function. Neotoma holds several diatom-related dataset types. Run neotoma2::get_table("datasettypes") to see the ones available to you. Here we keep diatom (sediment cores), diatom surface sample (modern samples used to calibrate relationships with the environment) and water chemistry. If we then want only the sediment cores that have a known chronology, we can filter by datasettype and by the presence of an age_range_young, which would indicate that there is a chronology that defines bounds for ages within the record.

3.3.1. Code

diatom_types <- c("diatom", "diatom surface sample", "water chemistry")

penn_diatoms <- penn_datasets %>%
  neotoma2::filter(datasettype %in% diatom_types)

penn_cores <- penn_datasets %>%
  neotoma2::filter(datasettype == "diatom", !is.na(age_range_young))

neotoma2::summary(penn_cores)

# We've removed records, so the new objects should be shorter than the original.
length(penn_diatoms) < length(penn_datasets)
length(penn_cores) < length(penn_diatoms)

3.3.2. Result

We can see now that the data table looks different (comparing it to the table above), and there are fewer total sites. Again, there is no explicit chronology for these records, we need to pull down the complete download for these records, but we begin to get a sense of what kind of data we have.

3.4. Pulling in sample() data

Because sample data adds a lot of overhead (the object that includes the datasets with samples is many times larger than the dataset alone), we try to call get_downloads() after we’ve done our preliminary filtering. After get_datasets() you have enough information to filter based on location, time bounds and dataset type. When we move to get_downloads() we can do more fine-tuned filtering at the analysis unit or taxon level.

We do not need every core for this, so we will work with one site. Look at the table of penn_cores in section 3.3: we pick Crow Hill (siteid 13319), a diatom record in Pennsylvania. We keep only its dataset with neotoma2::filter(), and then download it. The sequence is the same as before, from get_sites() to get_datasets() to filter() and finally get_downloads(), and because we are downloading a single dataset it only takes a few seconds.

crow_hill <- neotoma2::filter(penn_cores, siteid == 13319)
crow_dl <- neotoma2::get_downloads(crow_hill)

Once we’ve downloaded, we now have information for each site about all the associated collection units, the datasets, and, for each dataset, all the samples associated with the datasets. To extract all the samples we can call:

all_samp <- samples(crow_dl)

When we’ve done this, we get a data.frame that is 84 rows long and 39 columns wide. The reason the table is so wide is that we are returning data in a long format. Each row contains all the information you should need to properly interpret it:

##  [1] "age"             "agetype"         "ageolder"        "ageyounger"     
##  [5] "chronologyid"    "chronologyname"  "units"           "value"          
##  [9] "context"         "element"         "taxonid"         "symmetry"       
## [13] "taxongroup"      "elementtype"     "variablename"    "ecologicalgroup"
## [17] "analysisunitid"  "sampleanalyst"   "sampleid"        "depth"          
## [21] "thickness"       "samplename"      "datasetid"       "database"       
## [25] "datasettype"     "age_range_old"   "age_range_young" "age_units"      
## [29] "recdatecreated"  "datasetnotes"    "siteid"          "sitename"       
## [33] "lat"             "long"            "area"            "sitenotes"      
## [37] "description"     "elev"            "collunitid"

For some dataset types, or analyses some of these columns may not be needed, however, for other dataset types they may be critically important. To allow the neotoma2 package to be as useful as possible for the community we’ve included as many as we can.

In a diatom dataset each row is the count of valves (elementtype) of one taxon (variablename) in one sample (sampleid), and value holds the count. Here are the first rows of the columns we need most often:

all_samp %>%
  dplyr::select(sampleid, depth, age, variablename, value, units, elementtype) %>%
  head(10) %>%
  DT::datatable(rownames = FALSE, options = list(scrollX = "100%", dom = "t"))

Crow Hill has 2 samples, both modern (their ages are negative: years after 1950), so this is a snapshot of the diatoms in the lake today rather than a record through time. We will come back to that in section 4.

3.4.1. Extracting Taxa

If you want to know what taxa we have in the record you can use the helper function taxa() on the sites object. The taxa() function gives us, not only the unique taxa, but two additional columns, sites and samples that tell us how many sites the taxa appear in, and how many samples the taxa appear in, to help us better understand how common individual taxa are.

Because we only downloaded one site, the sites column is always 1 here. If we had downloaded several sites it would tell us how widespread each taxon is. The ecologicalgroup column is DIAT for every diatom taxon.

3.4.1.1. Code
neotomatx <- neotoma2::taxa(crow_dl)
3.4.1.2. Results

3.4.2. Understanding Taxonomies in Neotoma

Taxonomies in Neotoma are not as straightforward as we might expect. Diatom identification can be difficult: it depends on the orientation of the valve, the condition of the frustule, the optical quality of the microscope, the analyst’s training and the taxonomic reference they used. Diatom names have also changed a great deal over time, and different labs have identified the same valves at different taxonomic levels, from species to genus, or to a varietal level. You can read more about concepts of taxonomy within Neotoma in the Neotoma Manual’s section on Taxonomic concepts.

We use the unique identifiers (e.g., taxonid, siteid, analysisunitid) throughout the package, since they help us to link between records. The taxonid values returned by the taxa() call can be linked to the taxonid column in the samples() table. This allows us to link the taxa to other tables, as we will do in section 4.1. You may also note that the taxonname is in the field variablename. Individual sample counts are reported in Neotoma as variables. A “variable” may be either a species, something like laboratory measurements, or a non-organic proxy, like water chemistry measurements, and includes the units of measurement and the value.

3.4.2. Extracting Chronologies

Crow Hill has only two modern samples, and so do most of the Pennsylvania diatom records, which is what you want for a snapshot of the lakes today. To see change through time we need a sediment core, with many samples at different depths and ages. Among the diatom cores with a known chronology that we found in section 3.3 (penn_cores), four come from Barnegat Bay, on the New Jersey coast, and they have between 11 and 16 samples each.

Each chronology belongs to a collection unit. chronologies() gives us the chronology used by each dataset, so we can see how each core was dated:

3.4.2.1. Code
core_dl <- penn_cores %>% get_downloads(all_data = TRUE)

chronologies(core_dl) %>%
  as.data.frame() %>%
  dplyr::select(-notes)
3.4.2.2. Results

4. Simple Analytics

4.1. Taxonomy: Connecting to GBIF

The taxon names in a Neotoma download are the names that each analyst used, at the level they could identify (a species, a variety, a genus, or a project code such as Navicula sp. 3 BB DME). Diatom taxonomy is a complex topic, and names change over time, so sometimes we can redirect our names to a different authority.

We use GBIF, the Global Biodiversity Information Facility, because it is free, needs no key, and has an R package (rgbif) on CRAN. The same pattern works for other authorities, such as DiatomBase.

We do not send every taxon. Most of the 317 names in these four cores are rare. We use the common taxa of the four Barnegat Bay cores (at least 3% of the valves in three or more samples of a core), so the list stays small and we are gentle with the GBIF service. First, the valve counts of the cores that have at least five samples, as a proportion of each sample:

4.1.1.1. Code
long_cores <- samples(core_dl) %>%
  dplyr::filter(ecologicalgroup == "DIAT", elementtype == "valve") %>%
  group_by(datasetid) %>%
  dplyr::filter(n_distinct(sampleid) >= 5) %>%
  group_by(siteid, sitename, datasetid, sampleid, age) %>%
  mutate(proportion = value / sum(value)) %>%
  ungroup()
long_cores
4.1.1.2. Results

And then the common taxa, one row per name:

4.1.1.3. Code
common_taxa <- long_cores %>%
  dplyr::filter(proportion > 0.03) %>%
  group_by(siteid, variablename) %>%
  summarise(samples = n_distinct(sampleid), .groups = "drop") %>%
  dplyr::filter(samples >= 3) %>%
  dplyr::distinct(variablename) %>%
  arrange(variablename)
common_taxa
4.1.1.4. Results

That is 32 names, which account for 64% of all the valves counted in these cores.

4.1.2. One call to GBIF

rgbif can match a whole table of names against the GBIF taxonomy in one go with name_backbone_checklist(). We also give it the class, Bacillariophyceae (the diatoms). This matters, because the same genus name can be used for an animal or a plant, and with the class GBIF knows we mean the diatom. We make one call and save the answer with saveRDS(), so we never need to ask again. Please do not put this call in a loop.

4.1.2.1. Code
# The column has to be called `scientificName`, which is what rgbif looks for.
gbif_match <- rgbif::name_backbone_checklist(
  data.frame(
    scientificName = common_taxa$variablename,
    class = "Bacillariophyceae"
  )
)
4.1.2.2. Result

We did this ahead of time, on November 05, 2026. Here is the saved answer:

Remember to save the API’s response so we do not need to ask GBIF again.

saveRDS(
  list(retrieved = Sys.time(), match = gbif_match),
  "data/gbif_match_barnegat.RDS"
)

4.1.3. Merging the two

For each name GBIF tells us the taxon it matched (canonicalName) and at what rank; whether that name is ACCEPTED or a SYNONYM; how it matched (matchType: EXACT, FUZZY when the spelling was slightly different, or HIGHERRANK when it could only find a higher group); a confidence score; and the classification (family, order). The name that we sent comes back in verbatim_name, so we can use it to join the answer to our table with left_join(), which keeps every Neotoma taxon, with or without a match:

4.1.3.1. Code
gbif_small <- gbif_match %>%
  dplyr::select(
    variablename = verbatim_name,
    gbif_name = canonicalName, rank, status,
    match_type = matchType, accepted_name = acceptedScientificName,
    family, order
  )

taxa_gbif <- common_taxa %>%
  left_join(gbif_small, by = "variablename")

# A good check after any join: we should have the same number of rows.
nrow(taxa_gbif) == nrow(common_taxa)
4.1.3.1. Results

Of the 32 names, 27 matched exactly, and 5 could only be matched to a higher group. GBIF considers 6 of the names to be synonyms. Neither of these is an error, but both are worth looking at before you use the result. This is where the taxonomic decisions begin.

4.1.4. Synonyms

Synonyms. The name in Neotoma is not the one that GBIF treats as accepted today. Neither is wrong, and the accepted name is how you would link to other sources:

4.1.4.1. Code
taxa_gbif %>%
  dplyr::filter(status == "SYNONYM") %>%
  dplyr::select(variablename, accepted_name, family)
4.1.4.2. Results

4.1.5. Higher Rank

Matched only at a higher rank. GBIF could not find the species, so it gave us the genus instead. This is what happens with project codes. Be careful: a row that says Navicula is not a species match.

4.1.5.1. Code
taxa_gbif %>%
  dplyr::filter(match_type == "HIGHERRANK") %>%
  dplyr::select(variablename, gbif_name, rank)
4.1.5.2. Results

4.1.6. Harmonizing is up to you

Do you want to merge the synonyms into their accepted names? Fold the varieties into their species? Lump every Navicula with a project code into one group? Different researchers will answer differently, and that is fine. The easiest way to do it is not code: save the merged table as a csv, open it in a spreadsheet, and add a column with the name you want for each taxon.

The new column can then be joined to the samples() table by variablename, in the same way that we joined the GBIF answer, and you can group and sum the counts by your own names.

4.2. Stratigraphic Plotting

Now that we have looked at the taxon names, let’s look at how the taxa change through the depth of a single core. We can use packages like rioja to do stratigraphic plotting for a single record. We pick the Barnegat Bay core at West Creek (siteid 14270), and subset the core_dl object to that site:

4.2.1. Code

plotting_site <- neotoma2::filter(core_dl, siteid == 14270)

4.2.2. Results

rioja needs a “wide” table (one row per sample, one column per taxon), but samples() gives a “long” one. The toWide() function does the conversion, and it lets us choose the taxa, units and elements first. Here we keep only diatoms counted as valves (NISP), we group the samples by age, and we turn the counts into proportions. The first column of the result is age:

4.2.3. Code

counts <- plotting_site %>%
  samples() %>%
  toWide(
    ecologicalgroup = c("DIAT"),
    unit = c("NISP"),
    elementtypes = c("valve"),
    groupby = "age",
    operation = "prop"
  )

4.2.4. Results

There are hundreds of taxa in this core, and most of them are rare. We keep the age column, plus the taxa that make up at least 3% of the valves in three or more samples (the same rule we used for the common taxa in section 4.1, here for one core). (The ages here include negative values, after 1950, so we do not test the age column itself.)

4.2.5. Code

common <- colSums(counts[, -1] > 0.03, na.rm = TRUE) >= 3
counts <- counts[, c(TRUE, common)]

4.2.6. Result

I’ve also added a CONISS plot to the edge of the plot, to show how the new wide data frame works with distance metric functions. First we run the constrained clustering:

clust <- rioja::chclust(dist(sqrt(counts[, -1])), method = "coniss")
clust
## 
## Call:
## rioja::chclust(d = dist(sqrt(counts[, -1])), method = "coniss")
## 
## Cluster method   : coniss 
## Distance         : euclidean 
## Number of objects: 16

To plot the data we use rioja’s strat.plot(), converting proportions to percentages and sorting the taxa using weighted averaging scores (wa.order). Base R plots are drawn one layer on top of the other, so the zone lines and the site name are added in the same chunk:

4.2.7. Code

plot <- rioja::strat.plot(counts[, -1] * 100,
  yvar = counts$age,
  title = "",
  ylabel = "Calibrated Years BP (negative = after 1950)",
  xlabel = "Diatoms (% of valves counted)",
  y.rev = TRUE,
  clust = clust,
  wa.order = "topleft",
  scale.percent = TRUE,
  yTop = 0.62, # leave room for the taxon names
  cex.xlabel = 0.8
)

# Mark the two main zones of the clustering and name the site:
rioja::addClustZone(plot, clust, 2, col = "red")
mtext(plotting_site[[1]]$sitename,
  side = 3, line = -1.5, adj = 0.02, font = 2, outer = TRUE
)

4.2.8. Results

5. Conclusion

So, we’ve done a lot in this example. We’ve (1) searched for sites using site names and geographic parameters, (2) filtered results by dataset type and by whether a record has a chronology, (3) downloaded the sample information for one of the sites we found, and (4) performed basic analysis including a connection to GBIF to look up our diatom names (and a way to harmonize them yourself) and a stratigraphic diagram. Hopefully you can use these examples as templates for your own future work, or as a building block for something new and cool!

6. Installing packages on your own

We use several packages in this document, including leaflet, sf and others. We load the packages using the pacman package, which will automatically install the packages if they do not currently exist in your set of packages.

options(warn = -1)
pacman::p_load(
  neotoma2, dplyr, ggplot2, sf, geojsonsf, leaflet, DT, readr,
  rioja
)
# rgbif is not attached, because it masks neotoma2::datasets().
if (!requireNamespace("rgbif", quietly = TRUE)) install.packages("rgbif")

This workflow uses the development version of neotoma2 (version 1.1.1.9000 or later), which can search for several dataset types at once. You can install it with:

install.packages("remotes")
remotes::install_github("NeotomaDB/neotoma2", ref = "dev")

Note that R is sensitive to the order in which packages are loaded. Using neotoma2:: tells R explicitly that you want to use the neotoma2 package to run a particular function. So, for a function like filter(), which exists in other packages such as dplyr, you may see an error that looks like:

Error in UseMethod("filter") : 
  no applicable method for 'filter' applied to an object of class "sites"