Mytilus edulis size data

First load the required packages:

library(robis)
library(dplyr)
library(DT)
library(stringr)
library(ggplot2)
library(sf)
library(rnaturalearth)

Then use the occurrence() function with mof = TRUE to get all occurrences for Mytilus edulis, including the associated MeasurementOrFact data.

occ <- occurrence("Mytilus edulis", mof = TRUE)

The measurements() function can be used to create a flat measurements table. Use the fields parameter to include occurrence fields in this table.

mof <- measurements(occ, fields = c("decimalLongitude", "decimalLatitude", "collectionCode", "catalogNumber", "eventDate", "date_year", "dataset_id", "samplingProtocol", "samplingEffort"))

mof

## # A tibble: 281,670 x 23
##    decimalLongitude decimalLatitude collectionCode catalogNumber eventDate
##               <dbl>           <dbl> <chr>          <chr>         <chr>    
##  1             18.8            57.7 SHARK_Epibent… EVENT-1244:T… 2006-10-…
##  2             18.8            57.7 SHARK_Epibent… EVENT-1244:T… 2006-10-…
##  3             18.8            57.7 SHARK_Epibent… EVENT-1244:T… 2006-10-…
##  4             18.8            57.7 SHARK_Epibent… EVENT-1244:T… 2006-10-…
##  5             18.8            57.7 SHARK_Epibent… EVENT-1244:T… 2006-10-…
##  6             18.8            57.7 SHARK_Epibent… EVENT-1244:T… 2006-10-…
##  7             16.7            56.9 SHARK Zoobent… Shark zooben… 1990-05-…
##  8             16.7            56.9 SHARK Zoobent… Shark zooben… 1990-05-…
##  9             16.7            56.9 SHARK Zoobent… Shark zooben… 1990-05-…
## 10             16.7            56.9 SHARK Zoobent… Shark zooben… 1990-05-…
## # … with 281,660 more rows, and 18 more variables: date_year <int>,
## #   dataset_id <chr>, samplingProtocol <chr>, samplingEffort <chr>,
## #   measurementDeterminedBy <chr>, measurementAccuracy <chr>,
## #   measurementValue <chr>, measurementRemarks <chr>, measurementValueID <chr>,
## #   level <int>, occurrenceID <chr>, measurementUnit <chr>,
## #   measurementDeterminedDate <chr>, measurementType <chr>,
## #   measurementUnitID <chr>, measurementTypeID <chr>, measurementID <chr>,
## #   measurementMethod <chr>

This dataset includes a large variety of measurement types, let’s create a list:

types <- mof %>%
  group_by(measurementType, measurementTypeID, measurementUnit) %>%
  summarize(records = n(), datasets = length(unique(dataset_id))) %>%
  arrange(desc(records))

datatable(types)

There seem to be only two types of length measurements, both in mm. Let’s select these measurements by their measurementTypeID, and convert the measurement values to numbers.

size <- mof %>%
  filter(str_detect(measurementTypeID, "P01/current/OBSINDLX")) %>%
  mutate(measurementValue = as.numeric(measurementValue))
         
size

## # A tibble: 696 x 23
##    decimalLongitude decimalLatitude collectionCode catalogNumber eventDate
##               <dbl>           <dbl> <chr>          <chr>         <chr>    
##  1             3.32            51.6 RWS_2015_01 B… RWS2015_VOOR… 2015-04-…
##  2             4.02            51.9 Shellfish mon… IMA2016_1982… 2016-06-…
##  3             4.02            51.9 Shellfish mon… IMA2016_1982… 2016-06-…
##  4            -1.20            45.9 benthos_plume… Liens_bentho… 2008-02-…
##  5             4.12            52.1 Shellfish mon… IMA2016_1929… 2016-04-…
##  6             4.12            52.1 Shellfish mon… IMA2016_1929… 2016-04-…
##  7             4.12            52.1 Shellfish mon… IMA2016_1929… 2016-04-…
##  8             4.12            52.1 Shellfish mon… IMA2016_1929… 2016-04-…
##  9             4.12            52.1 Shellfish mon… IMA2016_1929… 2016-04-…
## 10             4.02            51.9 Shellfish mon… IMA2017_1982… 2017-06-…
## # … with 686 more rows, and 18 more variables: date_year <int>,
## #   dataset_id <chr>, samplingProtocol <chr>, samplingEffort <chr>,
## #   measurementDeterminedBy <chr>, measurementAccuracy <chr>,
## #   measurementValue <dbl>, measurementRemarks <chr>, measurementValueID <chr>,
## #   level <int>, occurrenceID <chr>, measurementUnit <chr>,
## #   measurementDeterminedDate <chr>, measurementType <chr>,
## #   measurementUnitID <chr>, measurementTypeID <chr>, measurementID <chr>,
## #   measurementMethod <chr>

Now we can plot the size measurements by collectionCode:

ggplot(size) +
  geom_jitter(aes(collectionCode, measurementValue, color = collectionCode), width = 0.1) +
  ylab("length (mm)")

And finally here are the measurements on a map:

world <- ne_countries(scale = "large", returnclass = "sf")

ggplot() + 
  geom_sf(data = world) +
  geom_point(data = size, aes(decimalLongitude, decimalLatitude, size = measurementValue, color = collectionCode), pch = 21) +
  scale_size(range = c(0.1, 8)) +
  coord_sf(xlim = range(size$decimalLongitude), ylim = range(size$decimalLatitude)) +
  theme_void()