Computational notebook for: Summary of benthic conditions in the Three Bays estuary (Cape Cod, MA) as of 2019

Process

\(~\)

Load data

# ------- 2019 benthic survey data --------------
# load 2019 survey data
data <- read_csv(here("data", "lab_and_field_summary.csv"), col_types = cols())

# reformat columns
data$time <- format(as.POSIXct(data$time), "%H:%M")
data$Secchi_depth_m <- as.numeric(data$Secchi_depth_m)

# add density of organisms
data$density_sq_m <- data$individuals/0.04                        # per square meter
data$density_count_1000s_sq_m <- data$density_sq_m/1000           # 1000s per sq m, for mapping

# load taxa and ecological groups
taxa <- read_csv(here("data", "taxa_and_ecogroups.csv"), skip = 1, col_types = cols())

# create a copy to reformat the data for tables
tabdata <- data

# ------- ancillary data --------------------

# rivers
p <- here::here("data", "ancillary", "3bays_major_rivers.shp")
rivers <- readOGR(dsn = p, layer = "3bays_major_rivers", verbose = FALSE)
rivers <- spTransform(rivers, CRS("+proj=longlat +ellps=GRS80"))

# NEERS data for Waquoit Bay 
needDates <- as.POSIXct(c("09/9/2019 09:00","09/11/2019 13:00"),format="%m/%d/%Y %H:%M",tzoffset=-5)
wb <- read_csv(here("data", "ancillary", "wqbmhwq2019.csv")) %>% 
  mutate(DateTimeStamp = as.POSIXct(DateTimeStamp,format="%m/%d/%Y %H:%M",tzoffset=-5)) %>% 
  filter(between(DateTimeStamp,needDates[1],needDates[2])) %>%
  filter(Sal > 25)   # filter out anomaly

# Tide series exported from NOAA Tides and currents, 6 m interval, meters, MLLW, LST/LDT 
# https://www.tidesandcurrents.noaa.gov/waterlevels.html?id=8449130&units=metric&bdate=20190909&edate=20190912&timezone=LST/LDT&datum=MLLW&interval=6&action=data
tide <- read.csv(here("data", "ancillary", "nantucket_tide_2019_09.csv"))
names(tide) <- c("Date","Time","Predicted","Preliminary","Verified")
tide <- tide %>% mutate(TimeStamp = as.POSIXct(paste(Date,Time),format="%Y/%m/%d %H:%M",tz="America/New_York"))

# tide phase delays for survey stations
tpd <- read.csv(here("data", "ancillary", "tide_phase_delays.csv")) %>%
  mutate(TimeStamp = as.POSIXct(paste(Sample.Date,Sample.Time),format="%m/%d/%Y %H:%M",tz="America/New_York"),
         Tide.TS=TimeStamp-Tide.Phase.Delay*60) # phase delay in minutes

Calculate summary statistics for water quality and sediments data
Calculate benthic infaunal community data metrics

# function to calculate AMBI, modified from Sigovini's M-AMBI function
AMBI <- function(x, eg) { 
  eg <- eg[which(rownames(eg) %in% colnames(x)), , drop = F]
  if(length(which(colnames(x) %in% rownames(eg))) < length(colnames(x))) {
    eg.na <- colnames(x)[-which(colnames(x) %in% rownames(eg))]
    eg <- rbind(eg, rep(NA, length(eg.na)))
    rownames(eg)[(nrow(eg)-length(eg.na)+1):nrow(eg)] <-
      colnames(x)[-which(colnames(x) %in% rownames(eg))]
  }
  eg <- eg[order(rownames(eg)), , drop = F]
  eg.matr <- as.data.frame(matrix(0, nrow = ncol(x), ncol = 6))
  colnames(eg.matr) <- c("EG1", "EG2", "EG3", "EG4", "EG5", "NA")
  rownames(eg.matr) <- rownames(eg)
  options(warn = -1)
  eg.matr[which(eg[, 1] == 1), 1] <- 1
  eg.matr[which(eg[, 1] == 2), 2] <- 1
  eg.matr[which(eg[, 1] == 3), 3] <- 1
  eg.matr[which(eg[, 1] == 4), 4] <- 1
  eg.matr[which(eg[, 1] == 5), 5] <- 1
  eg.matr[which(eg[, 1] == 0), 6] <- 1
  eg.matr[which(is.na(eg[, 1]) == T), 6] <- 1
  options(warn = 0)
  ambi <- as.data.frame(x %*% as.matrix(eg.matr))
  ambi <- cbind("N" = apply(x, 1, sum), ambi, "BC" = rep(NA, nrow(x)))
  for (i in 1:nrow(x)) {
    for (j in 2:6) {
      ambi[i, j] <- ambi[i, j] / (ambi[i, "N"] - ambi[i, "NA"]) * 100
    }
  }
  ambi$"NA" <- ambi$"NA"/ambi$"N"
  for (i in 1:nrow(x)) {
    ambi[i, "BC"] <- round((1.5*ambi[i, "EG2"] + 3*ambi[i, "EG3"] + 4.5*ambi[i, "EG4"] +
                              6*ambi[i, "EG5"]) / 100, 3)
  }
  colnames(ambi)[2:7] <- c("EG1(%)", "EG2(%)", "EG3(%)", "EG4(%)", "EG5(%)", "NA(%)")
  ambi[, 2:7] <- round(ambi[, 2:7], 2)
  return(ambi)
}

# convert taxa counts to community matrix format
taxamat <- taxa %>% select(STATION, SPECIES, COUNT) %>%
  pivot_wider(names_from = SPECIES, values_from = COUNT)%>%
  select(-STATION) %>%
  replace(is.na(.),0)

# create ecological group look-up table
eg <- taxa %>% 
  select(SPECIES, EG_genus)%>%
  filter(!duplicated(SPECIES))
eg <- data.frame(eg[,2], row.names = eg$SPECIES)

# calculate AMBI 
x <- as.matrix(taxamat)[order(as.numeric(rownames(taxamat))), order(colnames(taxamat))]
ambi <- AMBI(x, eg)
BC <- ambi$BC

# calculate other metrics for M-AMBI (and also evenness)
S <- specnumber(taxamat)
H1 <- diversity(taxamat, index = "shannon", MARGIN = 1, base = exp(1))
E <- round(H1/log(S),2)

# exclude sites from M-AMBI calculation when >50% of abundance cannot be assigned to an ecological group
BC[ambi$"NA(%)">0.5] <- NA

Calculate M-AMBI. Note that the supplemental script to calculate M-AMBI must be acquired from Sigovini et al., 2013 and sourced.

# specify order of metrics
ambi_var <- c("BC", "S", "H1")

# create dataframe of metrics in same order
metrics.ex <- data.frame(BC, S, H1)

# specify good and bad metric values (for polyhaline waters of the NE, as is Pelletier et al., 2018)
bad_metrics <- c(6,0,0)
good_metrics <- c(0.72, 44.0, 2.96)

# calculate m-ambi
m <- mambisimpl(m.values = metrics.ex, metrics = ambi_var, trasf = "f", high = good_metrics, bad = bad_metrics)

# select and clean up score results
m <- m[[2]]
m <- m[1:25,1:4]
m[, c(1,3,4)] <- round(m[, c(1,3,4)], 2)
m$`M-AMBI`[m$`M-AMBI`>1] <- 1               # scores above 1 are reset to 1

# assign scores to condition classes
m$condition <- "NA"
m$condition[m$`M-AMBI`<0.2] <- "bad"
m$condition[m$`M-AMBI`>=0.2] <- "poor"
m$condition[m$`M-AMBI`>=0.39] <- "moderate"
m$condition[m$`M-AMBI`>=0.53] <- "good"
m$condition[m$`M-AMBI`>0.77] <- "high"

## save habitat metrics
a <- data %>% select(c(station_id, date, time, lat, lon, depth_m, density_sq_m))
hab <- cbind(a, E, m)
#write_csv(hab, here("data", "habitat_metrics.csv"))

Summarize infaunal community metrics by subembayment and station type

hab <- read_csv(here("data", "habitat_metrics.csv"), col_types = cols())

# join community metrics with station info
a <- left_join(data %>% select(c(subembayment, station_id, station, type)), hab %>% select(c(station_id, density_sq_m:condition)))  %>% 
  arrange(subembayment)

# calculate M-AMBI score by subembayment and station type
b <- a %>% 
  group_by(subembayment, type) %>%
  summarise(avgMAMBI = mean(`M-AMBI`, na.rm=TRUE))

# classify M-AMBI score
b$condition <- "NA"
b$condition[b$avgMAMBI<0.2] <- "bad"
b$condition[b$avgMAMBI>=0.2] <- "poor"
b$condition[b$avgMAMBI>=0.39] <- "moderate"
b$condition[b$avgMAMBI>=0.53] <- "good"
b$condition[b$avgMAMBI>0.77] <- "high"

# summarize by station type
b2 <- a %>% 
  group_by(type) %>%
  summarise(d = mean(density_sq_m), S = mean(S), H1 = mean(H1), E = mean(E), M = mean(`M-AMBI`, na.rm = TRUE))

# join metrics with other data
a <- left_join(data, hab %>% select(c(station_id, S, 'M-AMBI', condition)))

Map the stations

# select columns to pop-up display
a <- a %>% select(-c(station_id, density_sq_m))

# subembayment labels
bays <- c("Warren's Cove", "Prince Cove", "North Bay", "Cotuit Bay", "West Bay", "Nantucket Sound")
lat <- c(41.6457, 41.64156, 41.63261, 41.62023, 41.61290, 41.605)
lon <- c(-70.405, -70.41223, -70.40658, -70.42815, -70.40108, -70.41)
bays <- data.frame(bays, lat, lon, stringsAsFactors = FALSE)

# river labels
rivnames <- c("Little R.", "Marstons Mills R.")
lat <- c(41.63, 41.651)
lon <- c(-70.424, -70.406)
rivnames <- data.frame(rivnames, lat, lon, stringsAsFactors = FALSE)

# create a color palette
pal <- colorFactor(
  palette = c('coral', 'deepskyblue2', 'yellow'),
  domain = a$type)

m = leaflet(a, height = 600, options = leafletOptions(zoomControl = TRUE)) %>% 
  addProviderTiles("CartoDB.PositronNoLabels") %>%
  setView(-70.411, 41.627, zoom = 13) %>%
  addCircleMarkers(lng = ~lon, lat = ~lat, weight = 2, radius = 6, opacity = 1, fillOpacity = 0.5, color = ~pal(type), fillColor = ~pal(type), label = ~station, popup = popupTable(data), labelOptions = labelOptions(noHide = TRUE, offset=c(0,18), textOnly = TRUE, direction = "top", style = list("color" = "black","font-weight" = "bold","font-size" = "9px"))) %>%
  addScaleBar("bottomleft") %>%
  addLegend("topleft", pal = pal, values = ~data$type, title = "Stations", opacity = 1) %>%
  addLabelOnlyMarkers(data = bays, lng = ~lon, lat = ~lat, label = ~bays, labelOptions = labelOptions(noHide = T, textOnly = T, offset=c(0,18), direction = "top", style = list("color" = "gray1","font-style" = "italic","font-size" = "10px"))) %>%
  addPolylines(data=rivers, weight=1, color="gray", opacity=0.35) %>%
  addLabelOnlyMarkers(data = rivnames, lng = ~lon, lat = ~lat, label = ~rivnames, labelOptions = labelOptions(noHide = T, textOnly = T, offset=c(0,18), direction = "top", style = list("color" = "gray45","font-style" = "italic","font-size" = "8px"))) %>%
  addMiniMap(
    tiles = providers$CartoDB.PositronNoLabels,
    position = "topright", 
    width = 80, height = 80, zoomLevelOffset = -7,
    toggleDisplay = FALSE)

m

Map selected observations

# code for sulfidic sites
a$sulfidic = 0
a$sulfidic[grep("sulf", a$comments)] <- 1

b <- a %>% select(c(station, type, lat, lon, Secchi_depth_m, shallow_salinity_ppt, `shallow_DO_%sat`, `fines_%`, `TOC_%`, sulfidic, density_count_1000s_sq_m, S, `M-AMBI`))

# select fixed stations and create column to indicate which stations change over time
fixed <- a %>% filter(type=="fixed") %>%
  mutate(change=0) %>%
  mutate(change=replace(change, station %in% c(7,9), -1)) 

# create long data for mapping
b <- b %>% pivot_longer(Secchi_depth_m:`M-AMBI`, names_to = "parameter", values_to = "value") %>%
  arrange(parameter, station)

# write a mapping function
obsmap <- function(data, fixed, rivers, palette, title, lbins = NULL, scale = FALSE, legend = TRUE, leglab = NULL){
  leaflet(data, options = leafletOptions(zoomControl = FALSE)) %>% 
  addProviderTiles("CartoDB.PositronNoLabels") %>%  
  setView(-70.411, 41.627, zoom = 13) %>%
  addCircleMarkers(lng = ~lon, lat = ~lat, weight = 1, radius = 6, opacity = 1, fillOpacity = 1, color = "gray", fillColor = ~pal(value), label = ~value) %>%
  addCircleMarkers(data = fixed, lng = ~lon, lat = ~lat, weight = 1, radius = 7, opacity = 1, fill = FALSE, color = "black") %>% 
  addPolylines(data=rivers, weight=1, color="gray", opacity=0.35) %>%
  {if(scale == TRUE)addScaleBar(.,"topright") else .} %>%
  {if(legend == TRUE)addLegend(., pal = palette, title = title, opacity = 1, bin = lbins, labels = label, values = ~value, position = "topleft") else addLegend(., pal = pal, title = ttl, opacity = 1, values = ~pal(mapdata$value), label = leglab, position = "topleft", labFormat = function(type, cuts, p) {paste0(leglab)})}
}

# salinity
ttl <- "Salinity (ppt)"
mapdata <- b %>% filter(parameter == "shallow_salinity_ppt")
bins <- c(22,24,26,28,30)
pal <- colorBin(brewer.pal(5, "BuPu")[2:5], domain = c(22,30), bins = bins)
m1 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = FALSE)
#m1

# secchi
ttl <- "Secchi<br>depth (m)"
mapdata <- b %>% filter(parameter == "Secchi_depth_m")
bins <- c(1, 1.5, 2, 2.5)
pal <- colorBin("YlGn", domain = c(1, 2.5), bins = bins)
m2 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = FALSE)
#m2

# do%sat
ttl <- "DO<sub>%sat"
mapdata <- b %>% filter(parameter == "shallow_DO_%sat")
pal <- colorBin("RdPu", bins = c(80, 90, 100, 110, 120, 130))
m3 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = TRUE)
#m3

# fines
ttl <- "Fines (%)"
mapdata <- b %>% filter(parameter == "fines_%")
pal <- colorBin("BuGn", domain = c(0, 100))
m4 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = FALSE)
#m4

# toc
ttl <- "TOC (%)"
mapdata <- b %>% filter(parameter == "TOC_%")
bins <- c(0, 1, 3.5, 5, 9.5)
pal <- colorBin("Oranges", domain = c(0, 9.5), bins = bins)
m5 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = FALSE)
#m5

# sulfidic
ttl <- "Sulfidic odor"
mapdata <- b %>% filter(parameter == "sulfidic" )
bins <- c(0, 0.5, 1)
pal <- colorBin("Greys", domain = c(0,1), bins = bins) 
leglab <- c("No", "Yes")
m6 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = TRUE, legend = FALSE, leglab = leglab)
#m6

# density
ttl <- "Density<br>(1000s/m<sup>2)"
mapdata <- b %>% filter(parameter == "density_count_1000s_sq_m")
bins <- c(0, 10, 20, 30, 40, 50)
pal <- colorBin("Reds", domain = c(0, 50), bins = bins)
m7 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = FALSE)
#m7

# species
ttl <- "Number <br>of species"
mapdata <- b %>% filter(parameter == "density_count_1000s_sq_m")
bins <- c(0, 10, 20, 30, 40, 50)
pal <- colorBin("PuRd", domain = c(0, 50), bins = bins)
m8 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = FALSE)
#m8

# m-ambi
worse <- fixed %>% filter(change<0)

ttl <- "M-AMBI <br> condition"
mapdata <- b %>% filter(parameter == "M-AMBI")
bins <- c(0, 0.2, 0.39, 0.53, 0.77, 1)
leglab <- c("Bad", "Poor", "Moderate", "Good", "High")
pal <- colorBin("RdYlBu", domain = c(0,1), bins = bins)
m9 <- obsmap(mapdata, fixed, rivers, pal, ttl, scale = TRUE, legend = FALSE, leglab = leglab)
m9 <- m9 %>% addCircleMarkers(data = worse, lng = ~lon, lat = ~lat, weight = 1, radius = 10, opacity = 1, fill = FALSE, color = "black")

sync(m1,m2,m3,m4,m5,m6,m7,m8,m9, ncol = 3, sync.cursor = FALSE)

—————————————

M-AMBI scores are binned into condition classes as follows (from Pelletier et al., 2018): Bad (<0.2), Poor (0.2-0.39), Moderate (0.39-0.53), Good (0.53-0.77), High (>0.77).

\(~\)

Comparison of overall habitat condition assessment

comp <- read_csv(here("data", "comparison.csv"), skip = 1, col_types = cols()) %>%
  select("Station ID_1", "Infaunal Indicators", "Habitat condition")

names(comp) <- c("station_id", "indicators", "condition")
comp$indicators[grep("H", comp$indicators)] <- "H or H/MI"
comp$condition[grep('Bad|Poor', comp$condition)] <- "Poor or Bad"

comp <- inner_join(data %>% select(station_id, lat, lon), comp)
comp$indicators <- as.factor(comp$indicators)
comp$condition <- as.factor(comp$condition)
#comp$condition <- ordered(comp$condition, levels = c("Good", "Moderate", "Poor", "Bad"))

colors <- brewer.pal(7, "RdYlBu")[c(1,4,7)]
ttl = "Prior survey <br> (early 2000s)"
pal = colorFactor(colors, levels = comp$indicators)

m10 <- leaflet(comp, options = leafletOptions(zoomControl = FALSE)) %>% 
  addProviderTiles("CartoDB.PositronNoLabels") %>%  
  setView(-70.411, 41.627, zoom = 13) %>%
  addCircleMarkers(lng = ~lon, lat = ~lat, weight = 1, radius = 6, opacity = 1, fillOpacity = 1, color = "gray", fillColor = ~pal(indicators), label = ~indicators) %>%
  addCircleMarkers(data = fixed, lng = ~lon, lat = ~lat, weight = 1, radius = 7, opacity = 1, fill = FALSE, color = "black") %>% 
  addPolylines(data=rivers, weight=1, color="gray", opacity=0.35) %>%
  addLegend(pal = pal, title = ttl, opacity = 1, values = ~indicators, position = "topleft") 

ttl = "Revisited sites <br> (2019)"
pal = colorFactor(colors, levels = comp$condition)

m11 <- leaflet(comp, options = leafletOptions(zoomControl = FALSE)) %>% 
  addProviderTiles("CartoDB.PositronNoLabels") %>%  
  setView(-70.411, 41.627, zoom = 13) %>%
  addCircleMarkers(lng = ~lon, lat = ~lat, weight = 1, radius = 6, opacity = 1, fillOpacity = 1, color = "gray", fillColor = ~pal(condition), label = ~condition) %>%
  addCircleMarkers(data = fixed, lng = ~lon, lat = ~lat, weight = 1, radius = 7, opacity = 1, fill = FALSE, color = "black") %>% 
  addPolylines(data=rivers, weight=1, color="gray", opacity=0.35) %>%
  addScaleBar("topright") %>%
  addLegend(pal = pal, title = ttl, opacity = 1, values = ~condition, position = "topleft") 

sync(m10, m11, ncol = 2, sync.cursor = FALSE)

—————————————

Color scales reflect equivalences that were made between the condition classes used in the prior MEP survey (“Significant Impairment (SI)”, “Moderate Impairment (MI)”, “Healthy habitat (H)”), and the respective groupings of M-AMBI classes (“Bad/Poor”, “Moderate”, “Good/High”).

\(~\)

Session Info

si <- sessioninfo::session_info()
pckgs <- map2(si$packages$package[si$packages$attached==TRUE], 
              si$packages$loadedversion[si$packages$attached==TRUE],
     ~ paste0(.x, " ", .y)) %>% 
  simplify()

OS: Windows 10 x64
Version: R version 4.0.3 (2020-10-10)
Packages: corrplot 0.84, dplyr 1.0.2, forcats 0.5.0, ggplot2 3.3.2, gsw 1.0-5, here 0.1, lattice 0.20-41, leaflet 2.0.3, leafpop 0.0.6, leafsync 0.1.0, mapview 2.9.0, oce 1.2-0, permute 0.9-5, purrr 0.3.4, RColorBrewer 1.1-2, readr 1.4.0, rgdal 1.5-18, sp 1.4-4, stringr 1.4.0, testthat 2.3.2, tibble 3.0.4, tidyr 1.1.2, tidyverse 1.3.0, vegan 2.5-7

Computational notebook for:

Summary of benthic conditions in the Three Bays estuary (Cape Cod, MA)
as of 2019

Laura Erban

2021-03-01

Overview

Process

—————————————

—————————————

Session Info