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

Written by E. Montes () and Eduardo Klein () on July 6, 2022.

This code pulls data from NOAA’s ERDDAP servers and creates time series plots of sea surface temperature (SST) and chlorophyll-a concentration (CHL), and maps showing the latest available data for the selected region.

Step 1

First, let’s load required libraries

library(readr)
library(rerddap)
library(lubridate)
library(dplyr)
library(flexdashboard)
library(reshape2)
library(leaflet)
library(ggplot2)
library(vegan)
library(xts)
library(dygraphs)
library(plotly)
library(mapdata)

library(RColorBrewer)
palette(brewer.pal(8, "Set2"))

Step 2

Query sea surface temperature (SST) data from ERDDAP. These are monthly SST means at ~ 1km pixel resolution starting in 2002.

## remove all spaces from string
NoSpaces = function(x){
  return(gsub(" ", "", x))
}

## set site coordinates and time for SST extraction. 
SSTSiteName = "Costa Das Algas"   ## for the resulting file name
SSTcoords.lon = -40.05  ## this is decimal longitude
SSTcoords.lat = -20.0  ## this is decimal latitude

SSTstartDate = "2017-01-01"  ## define start date of your time series 

## set dataset source (monthly SST)
SSTsource = info("jplMURSST41mday")

##
## Get sst 
SST <- griddap(SSTsource, 
              time=c(SSTstartDate, "last"),
              longitude = c(SSTcoords.lon,SSTcoords.lon),
              latitude = c(SSTcoords.lat,SSTcoords.lat),
              fields = "sst",
              fmt = "csv")

SST = SST[,c(1,4)]
names(SST) = c("time", "SST")

## convert time to a Data object
SST$time = as.Date(ymd_hms(SST$time))

Step 3

Plot SST time series

SST.xts = as.xts(SST$SST, SST$time)
dygraph(SST.xts, 
        ylab = "Sea Surface Temperature (Deg C)") %>% 
  dySeries("V1", label ="SST (Deg C)", color = "steelblue") %>%
  dyHighlight(highlightCircleSize = 5, 
              highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>% 
  dyOptions(fillGraph = FALSE, fillAlpha = 0.4) %>% 
  dyRangeSelector(dateWindow = c(min(SST$time), max(SST$time)))

Step 4

Save SST time series data as a CSV table

write_csv(SST, path = paste0(NoSpaces(SSTSiteName), "_SST.csv"))

Step 5

Create a map of the latest SST data. The blue dot indicates the location of the site where the time series was extracted in the code chunck above.

# define geographic domain
lat_lims<- c(-25., -15.)  ## change the first and second value for desired minimum and maximum latitude.
lon_lims <- c(-45., -35.)  ## change the first and second value for desired minimum and maximum longitude.

sstInfo <- info('jplMURSST41mday')
# get latest composite sst
GHRSST <- griddap(sstInfo, latitude = lat_lims, longitude = lon_lims, time = c('last','last'), fields = 'sst')

mycolor <- colors$temperature
w <- map_data("worldHires", ylim = lat_lims, xlim = lon_lims)
ggplot(data = GHRSST$data, aes(x = lon, y = lat, fill = sst)) + 
  geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
  geom_raster(interpolate = FALSE) +
  geom_point(x = SSTcoords.lon, y = SSTcoords.lat, colour = "blue", size = 3) +
  scale_fill_gradientn(colours = mycolor, na.value = NA) +
  theme_bw() + ylab("latitude") + xlab("longitude") +
  coord_fixed(1.1, xlim = lon_lims,  ylim = lat_lims) + ggtitle("Latest SST")

Step 6

Now, let’s query surface chlorophyll-a concentration (CHL; proxy for phytoplankton biomass) data from ERDDAP. These are weekly CHL means from the MODIS Aqua sensor at ~ 4km pixel resolution starting in 2003.

## remove all spaces from string
NoSpaces = function(x){
  return(gsub(" ", "", x))
}

## set site coordinates and time for CHL extraction
CHLSiteName = "Costa Das Algas"   ## for the resulting file name
CHLcoords.lon = -40.05 ## this is decimal longitude
CHLcoords.lat = -20.0 ## this is decimal latitude

CHLstartDate = "2017-01-01"

## set dataset source
CHLsource = info("erdMH1chla8day")

##
## Get CHL 
CHL <- griddap(CHLsource, 
               time=c(CHLstartDate, "last"),
               longitude = c(CHLcoords.lon,CHLcoords.lon),
               latitude = c(CHLcoords.lat,CHLcoords.lat),
               fields = "chlorophyll", fmt = "csv")

CHL = CHL[,c(1,4)]
names(CHL) = c("time", "CHL")
CHL = na.omit(CHL)

## convert time to a Data object
CHL$time = as.Date(ymd_hms(CHL$time))

Step 7

Plot CHL time series

CHL.xts = as.xts(CHL$CHL, CHL$time)
dygraph(CHL.xts, 
        ylab = "Chlorophyll a (mg m-3)") %>% 
  dySeries("V1", label ="CHL", color = "steelblue") %>%
  dyHighlight(highlightCircleSize = 5, 
              highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>% 
  dyOptions(fillGraph = FALSE, fillAlpha = 0.4) %>% 
  dyRangeSelector(dateWindow = c(min(CHL$time), max(CHL$time)))

Step 8

Save CHL time series data

write_csv(CHL, path = paste0(NoSpaces(CHLSiteName), "_CHL.csv"))

Step 9

Create a map of the latest CHL data. The blue dot indicates the location of the site where the time series was extracted in the code chunck above.

require("rerddap")
require("ggplot2")
require("mapdata")

# define geographic domain
lat_lims<- c(-25., -15.)  ## change the first and second value for desired minimum and maximum latitude.
lon_lims <- c(-45., -35.)  ## change the first and second value for desired minimum and maximum longitude.

chlaInfo <- info('nesdisVHNSQchlaMonthly')
CHLA <- griddap(chlaInfo, latitude = lat_lims, longitude = lon_lims, time = c('last','last'), fields = 'chlor_a')

# Map monthly chl (VIIRS)
mycolor <- colors$chlorophyll
w <- map_data("worldHires", ylim = lat_lims, xlim = lon_lims)
ggplot(data = CHLA$data, aes(x = lon, y = lat, fill = log(chlor_a))) + 
  geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") +
  geom_raster(interpolate = FALSE) +
  geom_point(x = CHLcoords.lon, y = CHLcoords.lat, colour = "blue", size = 3) +
  scale_fill_gradientn(colours = mycolor, na.value = NA) +
  theme_bw() + ylab("latitude") + xlab("longitude") +
  coord_fixed(1.3, xlim = lon_lims,  ylim = lat_lims) + ggtitle("Last month")