Marine Biodiversity Observation Network Pole to Pole of the Americas (MBON Pole to Pole)

Written by E. Montes () and Benjamin Best () on July 6, 2022.

This script pulls dynamic seascape classification data from NOAA’s CoastWatch and creates time series plots and summary tables showing seasacape relative dominance, and maps showing the latest available data for the selected region.

Step 1

First, let’s render a global map with the latest seascape classification

# Renders interactive map showing the latest seascapes globally.
library(rerddap)
library(tidyverse)
library(leaflet)
library(leafem)

seascapeInfo <- info("noaa_aoml_4729_9ee6_ab54", url = "https://cwcgom.aoml.noaa.gov/erddap/")  

# get the lastest monthly seascape classification
get_dates <- function(info){
  info$alldata$time %>%
    filter(attribute_name=="actual_range") %>%
    pull(value) %>%
    str_split(", ", simplify = T) %>%
    as.numeric() %>%
    as.POSIXct(origin = "1970-01-01", tz = "GMT")
}
d <- get_dates(seascapeInfo)[2]

# set lan/lon of site to center seascape map
site_lon <- -40.1
site_lat <- -20.0

# define the size of the domain
box <- 10

# render a map with the latest seascape
leaflet(
  options = leafletOptions(
    crs = leafletCRS(crsClass = "L.CRS.EPSG4326"))) %>%
  # basemap from GBIF in 4326
  addTiles("//tile.gbif.org/4326/omt/{z}/{x}/{y}@1x.png?style=gbif-geyser") %>%
  addWMSTiles(
    baseUrl = 'https://cwcgom.aoml.noaa.gov/erddap/wms/noaa_aoml_4729_9ee6_ab54/request?',
    layers = "noaa_aoml_4729_9ee6_ab54:CLASS",
    options = WMSTileOptions(
      version = "1.3.0", format = "image/png", transparent = T, opacity = 0.7,
      time = format(d,"%Y-%m-%dT00:00:00Z")))  %>%
  addMarkers(lng = site_lon, lat = site_lat)  %>%
  addMouseCoordinates() %>%
  fitBounds(site_lon - box, site_lat - box, site_lon + box, site_lat + box) %>%
  addLegend(
    position="bottomright",
    title = paste0("CLASS<br>", format(d,"%Y-%m-%d")),
    colorNumeric("Spectral", c(1,33), reverse=T), seq(1,33))

Step 2

Now let’s run an example with satellite seascapes using the seascapeR R package

# You may need to install the 'devtools' and 'seascapeR' packages as below:
# install.packages("devtools")
# library(devtools)
# remotes::install_github("marinebon/seascapeR")
# when prompted "Do you want to install from sources the packages which need compilation? (Yes/no/cancel)" select 'no'

library(seascapeR)
# variables
ss_dataset  = "global_monthly" # or "global_8day"
ss_var      = "CLASS"          # or "P"
date_beg    = "2020-01-01"
date_end    = "2022-06-01"

# get sanctuary polygon
lon = -40.; lat = -20.; w = 1.2
ply <- bbox_ply(lon - w, lat - w, lon + w, lat + w)

# get SeaScape dataset information
ss_info <- get_ss_info(dataset = ss_dataset)
ss_info
## <ERDDAP info> noaa_aoml_4729_9ee6_ab54 
##  Base URL: https://cwcgom.aoml.noaa.gov/erddap 
##  Dataset Type: griddap 
##  Dimensions (range):  
##      time: (2003-01-15T12:00:00Z, 2022-05-15T12:00:00Z) 
##      latitude: (-89.975, 89.975) 
##      longitude: (-179.975, 179.975) 
##  Variables:  
##      CLASS: 
##          Units: None 
##      P: 
##          Units: Punits

Step 3

Get SeaScape grids within polyon for date range

# get SeaScape grids within polyon for date range 
grds <- get_ss_grds(
  ss_info, ply, 
  ss_var    = ss_var, 
  date_beg  = date_beg, 
  date_end  = date_end)

# get first grid, a raster layer in the raster stack grds
grd <- raster::raster(grds, 100) # the number here is the index for the time period you would like to see mapped listed in 'tbl'

# map SeaScape grid
map_ss_grd(grd)

Step 4

Summarize SeaScape grids into a time series table and plot

tbl <- sum_ss_grds_to_ts(grds)
tbl
# plot SeaScape time series
plot_ss_ts(tbl, show_legend = "always")