Authors: The Palaeoverse Development Team
Last updated: 2024-09-16
rmacrostrat
is an R package that allows users to easily
retrieve geologic data from the Macrostrat database and facilitates
analyses of these data within the R environment. This vignette (or
tutorial, if you prefer) is provided to guide you through the
installation process and some of the functionality available within
rmacrostrat
. Specifically, we will focus on making a
geographic map of a specific geologic formation, a commonly required
(and very important) figure for a variety of purposes across the
geosciences.
The rmacrostrat
package can be installed via CRAN, or
its dedicated GitHub repository
if the development version is preferred. To install via CRAN, simply
use:
To install the development version, first install the
devtools
package, and then use install_github
to install rmacrostrat
directly from GitHub.
You can now load rmacrostrat
using the standard
library
function:
Before we get into the good stuff, the development team has a
small request. If you use rmacrostrat
in your
research, please cite the associated publication. This will help us to
continue our work in supporting you to do yours. You can access the
appropriate citation via:
## To cite rmacrostrat in publications, use the following citation:
##
## Jones, L.A., Dean, C.D., Gearty, W., and Allen, B.J. 2024. rmacrostrat: An R package for
## accessing and retrieving data from the Macrostrat geological database. EarthArXiv, 1–29,
## doi: 10.31223/X5XX37
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {rmacrostrat: An R Package for fetching geologic data from the Macrostrat database.},
## author = {Lewis A. Jones and Christopher D. Dean and William Gearty and Bethany J. Allen},
## year = {2024},
## journal = {EarthArXiv},
## pages = {1–29},
## doi = {10.31223/X5XX37},
## }
The Hell Creek is a geologic formation from the late Cretaceous and
early Paleogene of North America, which is found cropping out across
Montana and North and South Dakota, in the United States. It consists of
both freshwater and brackish clays, mudstones, and sandstones, deposited
by a fluvial system on the shoreline of the Western Interior Seaway. The
Hell Creek is famous for its preservation of an extensive array of both
terrestrial and marginal marine fossils from just prior to the K/Pg mass
extinction (66 million years ago), including widely recognizable
dinosaurs such as Edmontosaurus, Triceratops and
Tyrannosaurus rex. In this vignette, we will show you how to
use rmacrostrat
to find the outcrop of the Hell Creek
Formation, and use this to create a map which could be used in a
potential publication. Whilst we have chosen a specific formation, this
approach could be applied to any number of geologic units of interest.
So, let’s dive in!
First, we want to find the strat_name_id
s that are
associated with the Hell Creek. These IDs are linked to unique
stratigraphic names given to one or more geologic units, such as
Macrostrat units or map units, in the Macrostrat database. We can search
for the relevant IDs by using some of the def_*
functions.
These functions provide information about the data stored in the
database. Let’s have a look at the information associated with the name
Hell Creek, using def_strat_names()
:
## strat_name strat_name_long rank strat_name_id concept_id bed bed_id mbr mbr_id
## 1 Hell Creek Hell Creek Formation Fm 71267 8598 0 0
## 2 Hell Creek Hell Creek Member Mbr 71268 8598 0 Hell Creek 71268
## 3 Hell Creek Hell Creek Formation Fm 867 8598 0 0
## fm fm_id subgp subgp_id gp gp_id sgp sgp_id b_age t_age b_period t_period
## 1 Hell Creek 71267 0 Montana 2565 0 100.5 56 Cretaceous Paleogene
## 2 Lance 1090 0 Black Buttes coal 69391 0 100.5 56 Cretaceous Paleogene
## 3 Hell Creek 867 0 0 0 72.1 66 Cretaceous Paleogene
## c_interval t_units ref_id
## 1 0 19
## 2 0 19
## 3 13 1
We can see that there are three different strat_name
s
associated with the name Hell Creek. You’ll also see that all three have
the same concept_id
, showing that these stratigraphic names
are united under a single geologic entity. We can therefore use this
concept_id
to search for all references to the Hell
Creek:
# Search stratigraphic names for all associated with Hell Creek concept_id
hc_def <- def_strat_names(concept_id = 8598)
# View data
hc_def
## strat_name strat_name_long rank strat_name_id concept_id bed bed_id mbr mbr_id
## 1 Hell Creek Hell Creek Formation Fm 71267 8598 0 0
## 2 Hell Creek Hell Creek Member Mbr 71268 8598 0 Hell Creek 71268
## 3 Hell Creek Hell Creek Formation Fm 867 8598 0 0
## fm fm_id subgp subgp_id gp gp_id sgp sgp_id b_age t_age b_period t_period
## 1 Hell Creek 71267 0 Montana 2565 0 100.5 56 Cretaceous Paleogene
## 2 Lance 1090 0 Black Buttes coal 69391 0 100.5 56 Cretaceous Paleogene
## 3 Hell Creek 867 0 0 0 72.1 66 Cretaceous Paleogene
## c_interval t_units ref_id
## 1 0 19
## 2 0 19
## 3 13 1
Yup, these three are the only entries for Hell Creek. But hang on!
Before going ahead and using these to pull our spatial data, you’ll also
notice that one of the strat_name
s here is actually the
Hell Creek “Member”, which is described as belonging to the Lance
Formation. On a hunt in the literature, it turns out that this is a term
used for the Hell Creek before it was designated to the rank of
formation (Archibald, 1982), so it’s okay to include here. But it’s
always worth remembering to check your data for errors prior to using
it!
Let’s go ahead and use the strat_name_id
s associated
with the Hell Creek concept_id
to get some spatial data. To
do this, we’ll need to use the get_map_outcrop()
function,
which allows us to grab outcrop polygons from geologic maps. Let’s do
that now, specifying that we want the output to be returned as an sf, or
‘simple features’, object. An sf object is a collection of simple
features (representations of real world objects) that includes
attributes and geometries in the form of a data.frame
.
# Get spatial outcrop information associated with each Hell Creek strat_name_id
hc <- get_map_outcrop(strat_name_id = hc_def$strat_name_id, sf = TRUE)
Nice. Let’s have a look at the first entry to see what the data looks like.
## Simple feature collection with 1 feature and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -110.0199 ymin: 45.86057 xmax: -110.0168 ymax: 45.86277
## Geodetic CRS: WGS 84
## map_id source_id name strat_name lith descrip
## 1 2092425 25 Hell Creek Formation Hell Creek Formation Major:: {sandstone, mudstone}
## comments macro_units strat_names liths t_int_id t_int_age t_int_name b_int_id b_int_age
## 1 9213 867, 71267, 71268 7, 10 34 66 Late Cretaceous 34 100.5
## b_int_name color t_age b_age best_int_name geometry
## 1 Late Cretaceous #F2FA8C 66 72.1 Maastrichtian MULTIPOLYGON (((-110.0168 4...
That looks like what we need! We have a simple features collection consisting of multi-polygons, with attributes that indicate it contains data for the Hell Creek Formation. Now let’s move on to visualizing the data.
To plot our outcrop map, we’ll have to install and load some packages that let us handle and plot spatial data. Let’s do that now.
# Install rnaturalearthhires to use maps with state boundaries
devtools::install_github("ropensci/rnaturalearthhires")
# Load required packages
library(rnaturalearth)
library(rnaturalearthhires)
library(ggplot2)
library(sf)
First, we can use the packages rnaturalearth
and
rnaturalearthhires
to get the maps and state/province
outlines of the USA, Mexico and Canada to use as a background.
# Load country outlines
us <- ne_states(country = "united states of america", returnclass = "sf")
ca <- ne_states(country = "canada", returnclass = "sf")
mx <- ne_states(country = "mexico", returnclass = "sf")
Now we’re ready to plot our outcrop! We can use the
ggplot2
package to help us do this, in particular the
geom_sf()
function. Let’s plot our background maps and
outcrop now:
# Plot USA, Canada, and Hell Creek outcrop
ggplot() +
geom_sf(data = us) +
geom_sf(data = ca) +
geom_sf(data = hc, color = NA, fill = "#C7622B") +
coord_sf(xlim = c(-112, -97), ylim = c(44, 50), expand = FALSE) +
labs(x = "Longitude (º)", y = "Latitude (º)")
Looking good! But we can take this map a bit further. A common dinosaur found in the Hell Creek is Edmontosaurus, a large duck-billed dinosaur. Let’s combine what we’ve got so far with data from The Paleobiology Database (PBDB) to make a map of Edmontosaurus occurrences in the Hell Creek.
# Access PBDB API for Edmontosaurus occurrences from the Hell Creek
ed <- read.csv(paste0("https://paleobiodb.org/data1.2/occs/list.csv",
"?base_name=Edmontosaurus&strat=Hell%20Creek&show=full"))
# Turn points into spatial features
ed_sf <- st_as_sf(x = ed, coords = c("lng", "lat"))
# Set coordinate system (WGS84) for Edmontosaurus occurrences
st_crs(ed_sf) <- st_crs(4326)
# Plot
ggplot() +
geom_sf(data = us) +
geom_sf(data = ca) +
geom_sf(data = hc, color = NA, fill = "#C7622B") +
geom_sf(data = ed_sf, color = "black", size = 0.7) +
coord_sf(xlim = c(-112, -97), ylim = c(44, 50), expand = FALSE) +
labs(x = "Longitude (º)", y = "Latitude (º)")
Those Edmontosaurus occurrences are now plotted on the map. You may notice that there are some occurrences that don’t overlay the plotted outcrop, despite apparently being from the Hell Creek. This could be due to an error in the PBDB data, potentially relating to the accuracy of the latitude and longitude co-ordinates, or it could be that we’re missing some outcrop data of the Hell Creek. Either way, this shows that plotting occurrences over their respective outcrop can a good way to explore potential data issues.
Aside from those errant occurrences, it’s starting to look pretty
good! But we can add a few more finishing touches to really make it pop.
First, let’s the use package rphylopic
to add a silhouette of Edmontosaurus to the map.
# Get an appropriate silhouette for Edmontosaurus
ed_uuid <- get_uuid(name = "Edmontosaurus", n = 3)[3]
# Get the attribution information for the author of the image
get_attribution(uuid = ed_uuid, text = TRUE)
## Organism silhouettes are from PhyloPic (https://www.phylopic.org/; T. Michael Keesey, 2023) and were added using the rphylopic R package ver. 1.4.0.9000 (Gearty & Jones, 2023). Silhouette was made by Matt Dempsey, 2020 (CC BY 3.0). Silhouette was contributed by Matthew Dempsey.
# Setup data.frame for the silhouette (plot location)
silhouette_df <- data.frame(x = -99, y = 49.5, name = "Edmontosaurus")
# Plot the map
ggplot() +
geom_sf(data = us) +
geom_sf(data = ca) +
geom_sf(data = hc, color = NA, fill = "#C7622B") +
geom_sf(data = ed_sf, color = "black", size = 0.7) +
coord_sf(xlim = c(-112, -97), ylim = c(44, 50), expand = FALSE) +
geom_phylopic(data = silhouette_df, aes(x = x, y = y), uuid = ed_uuid,
size = 0.8, alpha = 1, color = "black") +
labs(x = "Longitude (º)", y = "Latitude (º)")
Nice! Finally, we can add a few other touches to make it a publication-worthy figure, including an inlay map, a scale bar, and a North arrow:
# Load `ggpatial` package to add annotations, and `cowplot` to stitch final map
# together
library(ggspatial)
library(cowplot)
# Make an inlay using the country maps we included earlier
inset <- ggplot() +
geom_sf(data = us) +
geom_sf(data = ca) +
geom_sf(data = mx) +
annotate("rect", xmin = -112, xmax = -97, ymin = 44, ymax = 50,
color = "red", fill = NA) +
xlim(-130, -70) +
ylim(30, 70) +
labs(x = NULL, y = NULL) +
theme_test() +
theme(axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
axis.title = element_blank(),
plot.margin = margin(0, 0, 0, 0, "cm"),
panel.background = element_rect(fill = "lightblue"))
# Add a scale bar, North arrow, and new theme to the original map and save it
ed_hc_map <- ggplot() +
geom_sf(data = us) +
geom_sf(data = ca) +
geom_sf(data = hc,
color = NA,
fill = "#C7622B") +
geom_sf(data = ed_sf,
color = "black",
size = 0.7) +
coord_sf(xlim = c(-112, -97), ylim = c(44, 50), expand = FALSE) +
annotation_north_arrow(location = "br",
pad_y = unit(0.75, "cm"),
height = unit(1, "cm"), width = unit(1, "cm")) +
annotation_scale(location = "br", width_hint = 0.3) +
geom_phylopic(data = silhouette_df, aes(x = x, y = y), uuid = ed_uuid,
size = 0.8, alpha = 1, color = "black") +
labs(x = "Longitude (º)", y = "Latitude (º)") +
theme_bw()
# Plot the two together
ggdraw() +
draw_plot(ed_hc_map) +
draw_plot(inset, 0.08, 0.63, width = 0.2, height = 0.2)