Connecting MERMAID ecological data with covariates example

Overview

This analysis explores the relationship between hard coral cover and degree heating weeks (DHW) in East Africa. We use the mermaidrcovariates package to access environmental data for MERMAID survey locations.

Degree Heating Weeks (DHW) measure accumulated heat stress rather than temperature at a single moment. They sum how far sea surface temperature has risen above the highest monthly mean expected at that location, across the preceding 12 weeks, so one DHW is equivalent to one week at 1 degree C above that threshold. NOAA Coral Reef Watch treats 4 DHW as the point at which significant bleaching becomes likely (Bleaching Alert Level 1), and 8 DHW as the point at which severe bleaching and mortality become likely (Alert Level 2). Those two values are used as thresholds throughout this document.

East Africa is represented here by surveys from Kenya and Tanzania.

Key questions:

• How does hard coral cover vary in East Africa?
• What is the exposure of different reefs to Degree Heating Weeks (DHW)?

Setup

library(here)

# Pin the project root explicitly. _quarto.yml sits in analysis/, which makes
# that folder look like a project root, so during a render here::here() would
# otherwise resolve to analysis/ instead of the repository root - and the
# cached data in data/ would not be found.
here::i_am("analysis/coral-cover-covariates-DHW-East-Africa.qmd")
library(tidyverse)
library(mermaidr)
library(plotly)
library(DT)
library(janitor)

# Install mermaidr.covariates (if needed)
# remotes::install_github("data-mermaid/mermaidr-covariates")

library(mermaidrcovariates)

# Set theme for ggplot
theme_set(theme_minimal(base_size = 12))

# Helper for the coloured status messages used throughout this document.
# knitr::asis_output() emits raw HTML on its own, so the chunk does NOT need
# `results: asis` and ordinary printed output (tibbles, data frames) in the
# same chunk still gets formatted normally.
status_msg <- function(..., color = "#006400") {
  knitr::asis_output(
    paste0('<span style="color: ', color, ';">', paste(...), '</span>\n\n')
  )
}

# Show file paths relative to the project root, so local absolute paths
# (C:/Users/.../) don't end up in the rendered HTML.
rel_path <- function(path) {
  sub(paste0(here::here(), "/"), "", path, fixed = TRUE)
}

# Set this to TRUE to force re-download of covariates
force_download <- FALSE

# --- Shared thresholds and colours -------------------------------------------
# Defined once here so every figure in this document uses the same values.

# NOAA Coral Reef Watch bleaching alert levels (degree heating weeks)
alert_1 <- 4   # Alert Level 1 - significant bleaching likely
alert_2 <- 8   # Alert Level 2 - severe bleaching and mortality likely

# Hard coral cover cut-points (percent) used to colour the cover histogram
cover_low  <- 10
cover_high <- 30

# Red / amber / green used throughout
col_red   <- '#d13823'
col_amber <- '#f3a224'
col_green <- '#277d1d'

# Named DHW bands, used for the map legend. Labels are built from the
# thresholds above so they cannot drift out of step with them.
dhw_bands <- setNames(
  c(col_green, col_amber, col_red),
  c(paste0("< ", alert_1, " DHW"),
    paste0(alert_1, " - ", alert_2, " DHW"),
    paste0("> ", alert_2, " DHW"))
)

Data Export

MERMAID Data

Export all publicly available summary sample events from MERMAID from Kenya and Tanzania.

Both the MERMAID export and the covariate extraction are cached to the project’s data/ folder the first time this document is rendered, and read from there on subsequent renders. Delete the relevant .rds file to force a fresh download - worth doing periodically for current data, though rebuilding the covariate extraction takes around 15 minutes.

# Get all public benthic data at the sample event level (all protocols)
mermaid_summ_ses_path <- here::here("data", "mermaid_summ_ses.rds")

if (file.exists(mermaid_summ_ses_path) && !force_download) {
  mermaid_summ_ses <- readRDS(mermaid_summ_ses_path)
  status_msg("Loaded cached summary sample events from",
             rel_path(mermaid_summ_ses_path),
             color = "#4682B4")
} else {
  mermaid_summ_ses <- mermaid_get_summary_sampleevents()
  dir.create(dirname(mermaid_summ_ses_path), showWarnings = FALSE, recursive = TRUE)
  saveRDS(mermaid_summ_ses, mermaid_summ_ses_path)
}

Loaded cached summary sample events from data/mermaid_summ_ses.rds

#file.exists("data/mermaid_summ_ses.rds")
#file.size("data/mermaid_summ_ses.rds")

mermaid_summ_ses <- mermaid_summ_ses %>% 
  filter(country %in% c("Kenya", "Tanzania"))

# Check how many sample events we have
status_msg("Total sample events:", nrow(mermaid_summ_ses))

Total sample events: 1291

Extract Hard Coral Cover

Hard coral cover is stored in the percent cover columns. We’ll focus on total hard coral cover (from all benthic protocols that provide coral cover) and show a histogram of the data.

The dotted lines at 10% and 30% mark the cut-points used to colour the bars: below 10% in red, 10-30% in amber, and above 30% in green.

# Extract relevant columns
coral_data <- mermaid_summ_ses %>% 
  rowwise() %>% 
  mutate(hard_coral_cover =
           mean(c(`benthicpit_percent_cover_benthic_category_avg_Hard coral`,
                  `benthiclit_percent_cover_benthic_category_avg_Hard coral`,
                  `benthicpqt_percent_cover_benthic_category_avg_Hard coral`,
                  quadrat_benthic_percent_percent_hard_avg_avg), na.rm = T)) %>% 
  ungroup() %>% 
  filter(!is.na(hard_coral_cover)) %>% 
  select(
    project,
    country,
    site,
    latitude,
    longitude,
    sample_date,
    hard_coral_cover
  ) %>% 
  # Convert date to Date format
  mutate(sample_date = as.Date(sample_date))

# Check how many sample events have hard coral cover estimates 
status_msg("Total sample events with hard coral data:", nrow(coral_data))

Total sample events with hard coral data: 789

# Create the histogram to get bin data to assign colors
hist_data <- hist(coral_data$hard_coral_cover,
                  breaks = seq(0, 100, by = 2),
                  plot = FALSE)

# Colour each bar by the cover value at its midpoint, not by its position in
# the sequence, so the colours stay tied to cover_low / cover_high even if the
# bin width above is changed.
bin_mid <- head(hist_data$breaks, -1) + diff(hist_data$breaks) / 2

bin_colors <- case_when(
  bin_mid < cover_low  ~ col_red,
  bin_mid < cover_high ~ col_amber,
  TRUE                 ~ col_green
)

# Create the histogram using plotly
hardCoralAggHist <-
  plot_ly(x = coral_data$hard_coral_cover,
          type = 'histogram',
          xbins = list(start = 0, size = 2, end = 100),
          marker = list(color = bin_colors), height = 450,
          hovertemplate = "Bin: %{x}%<br>%{y} surveys<extra></extra>") %>%
  config(displayModeBar = TRUE,
         displaylogo = FALSE,
         modeBarButtonsToRemove = c('zoom','pan', 'select', 'zoomIn', 'zoomOut',
                                    'autoScale', 'resetScale', 'lasso2d',
                                    'hoverClosestCartesian',
                                    'hoverCompareCartesian')) %>% 
  layout(bargap = 0.1,
         # Shapes and annotations below use "paper" coordinates on one axis:
         # 0 to 1 spans the plot area, and values above 1 sit above it (which
         # is how the title and subtitle are placed outside the plot).
         shapes = list(
           list(type = "line", 
                x0 = cover_low, x1 = cover_low, y0 = 0, y1 = 1, yref = "paper", 
                line = list(color = "black", dash = "dot")),
           list(type = "line", 
                x0 = cover_high, x1 = cover_high, y0 = 0, y1 = 1, yref = "paper", 
                line = list(color = "black", dash = "dot"))
         ),
         xaxis = list(title = "Hard coral cover (%)",
                      linecolor = "black",
                      linewidth = 2,
                      tickvals = seq(0, 100, by = 10),  # Set x-axis tick values
                      ticktext = seq(0, 100, by = 10)),
         yaxis = list(title = "Number of surveys",
                      linecolor = "black",   # Set the y-axis line color to black
                      linewidth = 2),
         annotations = list(
           list(x = 0, y = 1.15, text = "HARD CORAL COVER", showarrow = FALSE, 
                xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
                font = list(size = 20)),
           list(x = 0, y = 1.08,
                text = paste0(nrow(coral_data), " Surveys"),
                showarrow = FALSE, 
                xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
                font = list(size = 12))
         ),
         margin = list(t = 50, b = 75)) # Increase top margin to create more space for title and subtitle

# Visualize the plot
hardCoralAggHist

Get Covariates

Now we’ll use mermaidrcovariates to get DHW data.

List covariates

The first function of interest is list_covariates, which can be used to find out which covariates are available.

# Extract id and title from list_covariates() output
covariates <- list_covariates()

covariates_df <- covariates %>% 
  select(id, title, description)

# Create interactive table
datatable(
  covariates_df,
  options = list(
    pageLength = 10,
    scrollX = TRUE,
    autoWidth = TRUE
  ),
  caption = "Available MERMAID Covariates",
  rownames = FALSE,
  filter = "top"  # Adds search boxes at the top of each column
)

Get zonal statistics for DHW (raster)

The next function of interest is get_zonal_statistics, which can be used to extract raster data, like DHW, from the covariates. In this case, we calculate the mean DHW for 30 days previous to the sample date, as an example, but it can be used to extract data across any number of days prior to the survey using the n_days argument in the function.

Values are averaged within a 1 km radius of each survey location (radius = 1000, spatial_stats = "mean"). Since the underlying satellite product is gridded at 5 km, that radius in practice samples the pixel or pixels immediately covering the site rather than smoothing across neighbouring reefs.

# How the surveys break down by country
coral_data %>% 
  tabyl(country) %>% 
  adorn_pct_formatting(digits = 1) %>% 
  rename(Country = country, Surveys = n, `Percent of surveys` = percent) %>% 
  knitr::kable(align = c("l", "r", "r"))

Country	Surveys	Percent of surveys
Kenya	326	41.3%
Tanzania	463	58.7%

# get_zonal_statistics() needs the covariate's full title, so look it up by id
# rather than typing it out (pull() returns a plain string, not a tibble)
dhw_covariate <- covariates_df %>% 
  filter(id == "daily_dhw") %>% 
  pull(title)

status_msg("Using covariate:", dhw_covariate, color = "#4682B4")

Using covariate: Daily Global 5km Satellite Coral Bleaching Degree Heating Week

#set path, use if/else to read or cache the API call
covars_path <- here::here("data", "coral_data_dhw_30days.rds")

# Start timer
start_time <- Sys.time()

if (file.exists(covars_path) && !force_download) {
  coral_data_dhw <- readRDS(covars_path)
  status_msg("Loaded cached DHW data from",
             rel_path(covars_path),
             color = "#4682B4")
} else {
  coral_data_dhw <- coral_data %>%
    mermaidrcovariates::get_zonal_statistics(
      covariate = dhw_covariate,
      n_days = 30,
      radius = 1000,
      spatial_stats = "mean")
  
  dir.create(dirname(covars_path), showWarnings = FALSE, recursive = TRUE)
  
  saveRDS(coral_data_dhw, covars_path)
}

Loaded cached DHW data from data/coral_data_dhw_30days.rds

# Calculate elapsed time
end_time <- Sys.time()
elapsed_time <- end_time - start_time

# Print the result
status_msg("Data retrieval completed in:",
           round(elapsed_time, 2), attr(elapsed_time, "units"))

Data retrieval completed in: 0.02 secs

get_zonal_statistics() returns one row per day per survey, so a 30-day window gives a short time series for each sample event rather than a single number. summarise_zonal_statistics collapses that series to one value per survey. We take the maximum rather than the mean, because bleaching is driven by peak accumulated stress - averaging across the month would dilute a short, severe event into an unremarkable number.

The histogram below shows only those surveys that exceeded 4 DHW, rather than all surveys, so that the distribution of genuinely stressful conditions is visible.

coral_data_dhw_summary <- coral_data_dhw %>%
  select(site, sample_date, zonal_statistics) %>%
  unnest(zonal_statistics)

#names(coral_data_dhw_summary)


max_dhw <- coral_data_dhw %>%
  summarise_zonal_statistics("max")

#max_dhw

max_dhw_summary <- max_dhw %>%
  select(site, sample_date, summary_zonal_statistics) %>%
  unnest(summary_zonal_statistics)

#max_dhw_summary
#names(max_dhw_summary)

dhw_4_or_more <- max_dhw_summary %>% 
  filter(value > alert_1) %>% 
  arrange(-value)

#nrow(dhw_4_or_more) #74 sites with DHW>4 in 30 days prior to survey

# Create the histogram to get bin data to assign colors
# Bin the data manually (single source of truth for bin edges)
bin_width <- 0.5
bin_start <- 0
bin_end <- ceiling(max(dhw_4_or_more$value))

breaks <- seq(bin_start, bin_end, by = bin_width)

dhw_binned <- dhw_4_or_more %>%
  mutate(bin = cut(value,
                   breaks = breaks,
                   right = FALSE,
                   include.lowest = TRUE)) %>%
  count(bin, .drop = FALSE) %>%
  mutate(
    bin_left = breaks[-length(breaks)],
    bin_mid  = bin_left + bin_width / 2,
    bin_color = if_else(bin_mid < alert_2, col_amber, col_red)
  )

# Create the histogram using plotly (as a bar chart on pre-binned data)
dhwAggHist <-
  plot_ly(data = dhw_binned, height = 450) %>%
  add_bars(x = ~bin_left + bin_width / 2,  # plot at bin midpoint
           y = ~n,
           width = bin_width * 0.9,        # bar width in x-axis units (trace attribute)
           marker = list(color = ~bin_color),
           hovertemplate = "Bin: %{x} DHW<br>%{y} surveys<extra></extra>") %>%
  config(displayModeBar = TRUE,
         displaylogo = FALSE,
         modeBarButtonsToRemove = c('zoom','pan', 'select', 'zoomIn', 'zoomOut',
                                    'autoScale', 'resetScale', 'lasso2d',
                                    'hoverClosestCartesian',
                                    'hoverCompareCartesian')) %>% 
  layout(bargap = 0.1,
         shapes = list(
           list(type = "line", 
                x0 = alert_2, x1 = alert_2, y0 = 0, y1 = 1, yref = "paper", 
                line = list(color = "black", dash = "dot"))
         ),
         xaxis = list(title = "Degree Heating Weeks (DHW)",
                      linecolor = "black",
                      linewidth = 2),
         yaxis = list(title = "Number of surveys",
                      linecolor = "black",
                      linewidth = 2),
         annotations = list(
           list(x = 0, y = 1.15, text = "DEGREE HEATING WEEKS", showarrow = FALSE, 
                xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
                font = list(size = 20)),
           list(x = 0, y = 1.08,
                text = paste0(nrow(dhw_4_or_more), " Surveys"),
                showarrow = FALSE, 
                xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
                font = list(size = 12))
         ),
         margin = list(t = 50, b = 75))

dhwAggHist

Exposure to DHW by coral cover

Joining coral cover to heat stress needs some care. A single sample event can appear more than once - once per benthic protocol, for instance - so both tables are reduced to one row per site and date before joining, averaging where duplicates remain.

Each point below is one survey. The shaded bands and dotted lines mark the same 4 and 8 DHW thresholds described above: green below 4, amber from 4 to 8, and red above 8.

#First make sure there are no duplicates in the max_dhw_summary data to join
max_dhw_summary_join <- max_dhw_summary |>
  filter(n_dates != 0) |> #get rid of any without data
  select(site, sample_date, covariate, value) |> 
  distinct() |> 
  group_by(site, sample_date, covariate) |> 
  summarise(value = mean(value)) |> 
  ungroup()

#Also remove duplicates in the coral data
coral_data_join <- coral_data |> 
  distinct() |> 
  group_by(project, country, site, latitude, longitude, sample_date) |> 
  summarise(hard_coral_cover = mean(hard_coral_cover)) |> 
  ungroup()

coral_data_with_covariate <- coral_data_join %>% 
  left_join(max_dhw_summary_join)

# Determine y-axis max to size the top band appropriately
y_max <- max(coral_data_with_covariate$value, na.rm = TRUE) * 1.05

# Create the scatter plot using plotly
coralDhwScatter <-
  plot_ly(data = coral_data_with_covariate,
          x = ~hard_coral_cover,
          y = ~value,
          type = 'scatter',
          mode = 'markers',
          marker = list(color = "black",
                        size = 8,
                        opacity = 0.6),
          height = 450,
          text = ~paste0(site, ", ", country, "<br>", format(sample_date, "%b %d, %Y")),
          hovertemplate = "%{text}<br>Hard coral cover: %{x:.1f}%<br>DHW: %{y:.1f}<extra></extra>") %>%
  config(displayModeBar = TRUE,
         displaylogo = FALSE,
         modeBarButtonsToRemove = c('zoom','pan', 'select', 'zoomIn', 'zoomOut',
                                    'autoScale', 'resetScale', 'lasso2d',
                                    'hoverClosestCartesian',
                                    'hoverCompareCartesian')) %>% 
  layout(shapes = list(
           # Background bands, drawn first so they sit behind the points.
           # x0/x1 with xref = "paper" makes each band span the full plot
           # width; y0/y1 are in real DHW units.
           list(type = "rect",
                x0 = 0, x1 = 1, xref = "paper",
                y0 = 0, y1 = alert_1,
                fillcolor = col_green, opacity = 0.12, line = list(width = 0),
                layer = "below"),
           list(type = "rect",
                x0 = 0, x1 = 1, xref = "paper",
                y0 = alert_1, y1 = alert_2,
                fillcolor = col_amber, opacity = 0.15, line = list(width = 0),
                layer = "below"),
           list(type = "rect",
                x0 = 0, x1 = 1, xref = "paper",
                y0 = alert_2, y1 = y_max,
                fillcolor = col_red, opacity = 0.15, line = list(width = 0),
                layer = "below"),
           # Threshold lines
           list(type = "line", 
                x0 = 0, x1 = 1, xref = "paper", y0 = alert_1, y1 = alert_1, 
                line = list(color = "black", dash = "dot")),
           list(type = "line", 
                x0 = 0, x1 = 1, xref = "paper", y0 = alert_2, y1 = alert_2, 
                line = list(color = "black", dash = "dot"))
         ),
         xaxis = list(title = "Hard coral cover (%)",
                      linecolor = "black",
                      linewidth = 2,
                      range = c(0, 100)),
         yaxis = list(title = "DHW exposure (past 30 days)",
                      linecolor = "black",
                      linewidth = 2,
                      range = c(0, y_max)),
         annotations = list(
           list(x = 0, y = 1.15, text = "CORAL COVER VS. HEAT STRESS", showarrow = FALSE, 
                xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
                font = list(size = 20)),
           list(x = 0, y = 1.08,
                text = paste0(nrow(coral_data_with_covariate), " Surveys"),
                showarrow = FALSE, 
                xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
                font = list(size = 12))
         ),
         margin = list(t = 50, b = 75))

coralDhwScatter

Coral cover and heat stress by year

Splitting the same data by year separates chronic patterns from individual heat stress events. Only years in which at least one survey exceeded 4 DHW are shown - years with no recorded heat stress are omitted rather than drawn as empty panels.

# Add year column (in case it wasn't already created upstream)
coral_data_with_covariate <- coral_data_with_covariate %>%
  mutate(year = lubridate::year(sample_date))

# Recreate y_max and bands so this chunk doesn't depend on earlier chunks
y_max <- max(coral_data_with_covariate$value, na.rm = TRUE) * 1.05

bands <- tibble(
  ymin = c(0, alert_1, alert_2),
  ymax = c(alert_1, alert_2, y_max),
  fill = c(col_green, col_amber, col_red)
)

# Identify years where at least one survey had DHW > 4
years_with_heat_stress <- coral_data_with_covariate %>%
  filter(value > alert_1) %>%
  distinct(year) %>%
  pull(year)

# Filter to only those years (keeping all points within them)
coral_data_heat_years <- coral_data_with_covariate %>%
  filter(year %in% years_with_heat_stress)

coralDhwFacetGG <-
  ggplot(coral_data_heat_years, aes(x = hard_coral_cover, y = value)) +
  geom_rect(data = bands, inherit.aes = FALSE,
            aes(xmin = -Inf, xmax = Inf, ymin = ymin, ymax = ymax, fill = fill),
            alpha = 0.15) +
  scale_fill_identity() +
  geom_hline(yintercept = c(alert_1, alert_2), linetype = "dotted", color = "black") +
  geom_point(color = "black", alpha = 0.6, size = 1.8) +
  facet_wrap(~ year, ncol = 3) +
  coord_cartesian(xlim = c(0, 100), ylim = c(0, y_max)) +
  labs(title = "CORAL COVER VS. HEAT STRESS, BY YEAR",
       subtitle = paste0(nrow(coral_data_heat_years), " Surveys (years with DHW > 4 only)"),
       x = "Hard coral cover (%)",
       y = "DHW exposure (past 30 days)") +
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    strip.background = element_rect(fill = "grey90", color = NA),
    strip.text = element_text(face = "bold"),
    plot.title = element_text(size = 20, face = "bold"),
    plot.subtitle = element_text(size = 12, color = "grey30"),
    axis.line = element_line(color = "black", linewidth = 0.6)
  )

coralDhwFacetGG

ggsave(
  filename = here::here("plots", "coral_cover_heat_stress_by_year.png"),
  plot     = coralDhwFacetGG,
  width    = 9,
  height   = 6,
  dpi      = 300,
  bg       = "white"
)

Survey sites and heat stress exposure

Finally, the same surveys mapped in space, using the DHW banding from the figures above. The dropdown filters to a single survey year, since many sites were surveyed repeatedly and would otherwise plot on top of each other.

#This block creates a map of the sites, colored by the number of DHW

# dhw_bands, alert_1 and alert_2 are defined once in the setup chunk
map_data <- coral_data_with_covariate %>%
  mutate(year = lubridate::year(sample_date)) %>%   # in case it is not already there
  filter(!is.na(value), !is.na(latitude), !is.na(longitude)) %>%
  mutate(
    dhw_band   = cut(value,
                     breaks = c(-Inf, alert_1, alert_2, Inf),
                     labels = names(dhw_bands),
                     right  = FALSE),
    band_color = unname(dhw_bands[as.character(dhw_band)]),
    hover_text = paste0(site, ", ", country,
                        "<br>", format(sample_date, "%b %d, %Y"),
                        "<br>Max DHW: ", round(value, 1),
                        "<br>Hard coral cover: ", round(hard_coral_cover, 1), "%")
  )

map_years <- sort(unique(map_data$year))
n_legend  <- length(dhw_bands)
n_years   <- length(map_years)

dhwSiteMap <- plot_ly(height = 600)

# Data-free traces that exist only to hold the colour legend in place,
# so the key stays visible whichever year is selected
for (b in names(dhw_bands)) {
  dhwSiteMap <- dhwSiteMap %>%
    add_trace(type = "scattermapbox", mode = "markers",
              lat = NA, lon = NA,
              marker = list(size = 10, color = unname(dhw_bands[b])),
              name = b, showlegend = TRUE, hoverinfo = "skip")
}

# One trace per year, each point coloured by its DHW band
for (y in map_years) {
  d <- map_data %>% filter(year == y)
  dhwSiteMap <- dhwSiteMap %>%
    add_trace(type = "scattermapbox", mode = "markers",
              lat = d$latitude, lon = d$longitude,
              marker = list(size = 9, opacity = 0.8, color = d$band_color),
              text = d$hover_text,
              hovertemplate = "%{text}<extra></extra>",
              name = as.character(y), showlegend = FALSE)
}

# Dropdown entries: "All years" first, then one per year.
# Each sets the visible flag for every trace, in order.
year_buttons <- list(
  list(method = "restyle", label = "All years",
       args = list("visible", rep(TRUE, n_legend + n_years)))
)

for (i in seq_along(map_years)) {
  vis <- c(rep(TRUE, n_legend), rep(FALSE, n_years))
  vis[n_legend + i] <- TRUE
  year_buttons <- c(year_buttons,
                    list(list(method = "restyle",
                              label = as.character(map_years[i]),
                              args  = list("visible", vis))))
}

dhwSiteMap <- dhwSiteMap %>%
  config(displayModeBar = TRUE,
         displaylogo = FALSE,
         modeBarButtonsToRemove = c('select', 'lasso2d',
                                    'hoverClosestCartesian',
                                    'hoverCompareCartesian')) %>%
  layout(
    mapbox = list(
      style  = "open-street-map",   # free style, no Mapbox token needed
      center = list(lat = mean(map_data$latitude, na.rm = TRUE),
                    lon = mean(map_data$longitude, na.rm = TRUE)),
      zoom   = 5
    ),
    legend = list(title = list(text = "<b>Max DHW</b>"),
                  x = 0.01, y = 0.99,
                  bgcolor = "rgba(255,255,255,0.8)"),
    updatemenus = list(
      list(type = "dropdown", direction = "down", showactive = TRUE,
           x = 1, xanchor = "right", y = 1.10, yanchor = "top",
           buttons = year_buttons)
    ),
    annotations = list(
      list(x = 0, y = 1.14, text = "SURVEY SITES BY HEAT STRESS", showarrow = FALSE,
           xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
           font = list(size = 20)),
      list(x = 0, y = 1.06,
           text = paste0(nrow(map_data), " Surveys, ", n_years, " years"),
           showarrow = FALSE,
           xref = 'paper', yref = 'paper', xanchor = 'left', yanchor = 'top',
           font = list(size = 12))
    ),
    margin = list(t = 95, b = 10, l = 10, r = 10))

dhwSiteMap

Data Access

This analysis uses:

• MERMAID data: Publicly available benthic sample events for Kenya and Tanzania, accessed via the mermaidr package
• Environmental covariates: Accessed via the mermaidrcovariates package
  • DHW: Daily Global 5km Satellite Coral Bleaching Degree Heating Week, maximum over the 30 days prior to each survey, within a 1 km radius

Session Info

sessionInfo()

R version 4.6.1 (2026-06-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mermaidrcovariates_1.0.2 janitor_2.2.1            DT_0.34.0               
 [4] plotly_4.12.1            mermaidr_1.3.1           lubridate_1.9.5         
 [7] forcats_1.0.1            stringr_1.6.0            dplyr_1.2.1             
[10] purrr_1.2.2              readr_2.2.0              tidyr_1.3.2             
[13] tibble_3.3.1             ggplot2_4.0.3            tidyverse_2.0.0         
[16] here_1.0.2              

loaded via a namespace (and not attached):
 [1] gtable_0.3.6        bslib_0.12.0        xfun_0.60          
 [4] htmlwidgets_1.6.4   tzdb_0.5.0          vctrs_0.7.3        
 [7] tools_4.6.1         crosstalk_1.2.2     generics_0.1.4     
[10] curl_8.0.0          proxy_0.4-29        pkgconfig_2.0.3    
[13] KernSmooth_2.23-26  data.table_1.18.6.1 RColorBrewer_1.1-3 
[16] S7_0.2.2            lifecycle_1.0.5     compiler_4.6.1     
[19] farver_2.1.2        textshaping_1.0.5   snakecase_0.11.1   
[22] httpuv_1.6.17       sass_0.4.10         htmltools_0.5.9    
[25] rstac_1.0.1         class_7.3-23        yaml_2.3.12        
[28] jquerylib_0.1.4     later_1.4.8         pillar_1.11.1      
[31] crayon_1.5.3        classInt_0.4-11     cachem_1.1.0       
[34] tidyselect_1.2.1    digest_0.6.39       stringi_1.8.9      
[37] sf_1.1-2            labeling_0.4.3      rprojroot_2.1.1    
[40] fastmap_1.2.0       grid_4.6.1          cli_3.6.6          
[43] magrittr_2.0.5      e1071_1.7-17        withr_3.0.3        
[46] scales_1.4.0        promises_1.5.0      timechange_0.4.0   
[49] rmarkdown_2.31      httr_1.4.8          jpeg_0.1-11        
[52] otel_0.2.0          ragg_1.5.2          png_0.1-9          
[55] hms_1.1.4           evaluate_1.0.5      knitr_1.51         
[58] viridisLite_0.4.3   rlang_1.3.0         Rcpp_1.1.2         
[61] glue_1.8.1          DBI_1.3.0           rstudioapi_0.19.0  
[64] jsonlite_2.0.0      R6_2.6.1            systemfonts_1.3.2  
[67] units_1.0-1