In pbs-assess/PACea: An R package of Pacific ecosystem information to help facilitate an ecosystem approach to fisheries management

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 5.7,
  fig.height = 5
)

Daily average sea surface temperature from buoys

library(pacea)
library(dplyr)
library(tibble)  # Else prints all of a tibble
library(ggplot2)

In pacea we include daily average sea surface temperature (SST) that we have calculated from measurement from r nrow(buoy_metadata) buoys in Canadian Pacific waters, yielding over 200,000 daily means. Data are from Environment and Climate Change Canada, and Fisheries and Oceans Canada. The earliest data are from September 1987, and 14 buoys were still providing data as of May 2023.

Metadata for the buoys is given by

buoy_metadata

which includes each buoy's World Meteorological Organisation's weather station id (wmo_id), its 'common name' (name), location, and depth of the water in which the buoy resides. See ?buoy_metadata for descriptions of all columns. The full names of all buoys are

buoy_metadata$name

The locations of the buoys are given by

knitr::include_graphics(paste0(here::here(),
                               "/man/figures/buoys_map.png"))

(built with code in data-raw/buoys/buoys-map.R), and also shown at https://github.com/IOS-OSD-DPG/Pacific_SST_Monitoring#eccc-buoy-data.

Some of the recent data are also plotted on that website, by Andrea Hilborn, Charles Hannah and Lu Guan, which is automatically updated roughly every week (so look there for a quick glance at recent conditions). For pacea we have, with Andrea, adapted some of that code to include the data and create plotting functions in our package.

The wrangling of data is taken care of within pacea, and includes using protocols to remove certain flagged data, dealing with timezones and pesky daylight savings time changes, and averaging over a day (the original raw data are even higher resolution). We are still refining this to remove outliers -- this has to be somewhat manual.

The sst value are saved in a tibble, which also has class pacea_buoy

buoy_sst
tail(buoy_sst)

with stn_id specifying each buoy as described above in buoy_metadata, and date based on UTC -8 hours (i.e. Pacific Standard Time, not changing due to daylight savings), and sst is the mean SST (deg C) for that day at that station; see ?buoy_sst for full details.

To see the ranges of dates for each buoy, use dplyr functions in the usual way:

buoy_ranges <- buoy_sst %>%
  group_by(stn_id) %>%
  summarise(start = min(date),
            end = max(date))
buoy_ranges

Plot data from a single buoy for all years, the default is for buoy C46205:

plot(buoy_sst)

The red line gives the current year. Since buoy_sst has class pacea_buoy, plot(buoy_sst) function calls our tailored plot.pacea_buoy() function. The gap in the red line at the start of the year indicates missing data (or days that get excluded due to our protocols).

To see data from each buoy in turn:

for(i in 1:nrow(buoy_metadata)){
  plot(buoy_sst,
       stn_id = buoy_metadata[i, ]$stn_id) %>% print()
}

We might adapt Andrea Hilborn's code (saved below in the .Rmd) to produce a panel plot of all buoys at once, to give a figure resembling that at https://github.com/IOS-OSD-DPG/Pacific_SST_Monitoring/blob/main/figures/current/Daily_mean_buoy_overview_2023.png. But for pacea users the individual buoys are probably sufficient.

# Andrea's, and see below:
ggplot(rbind(dat1,dat2), aes(month(date, label=TRUE, abbr=TRUE),
                value, group=factor(year(date)), colour=factor(year(date)))) +

# First stab, to work on if desired.
p <- buoy_sst %>%      # then filter as necessary given function options. Try
                       # and just do all for now.
  arrange(stn_id) %>%      # maybe not needed, or do up front in data-raw
  ggplot() +
  facet_wrap(~stn_id, ncol = 3) +
  geom_ribbon(data = buoy_sst,
              aes(x = date,
                  y = sst,
                  ymin = 0,
                  ymax = 20))
p

                  )) # ymin = SSTP_10perc_day, ymax = SSTP_90perc_day), fill = "grey70", colour = NA, alpha = 0.5) +
  ## geom_path(data = sstclim, aes(x = fakedate, y = SSTP_clim_mean_day), colour = "grey30", linewidth = 0.8) +
  ## # Heat dome - for 2021-2022 plots
  ## # geom_vline(xintercept = as.Date("2021-06-26"), linetype = "dashed") +
  ## geom_path(data = prevyr %>% filter((year(date) == plot_yrs[1]) | is.na(SSTP)), aes(x = fakedate, y = SSTP, colour = col_key), size = 0.6) +

  ## geom_path(data = currentyr %>% filter((year(date) == plot_yrs[2]) | is.na(SSTP)),aes(x = fakedate, y = SSTP), colour = "black", size = 1.4, lineend = "round") +
  ## geom_path(data = currentyr %>% filter(year(date) == plot_yrs[2] | is.na(SSTP)), aes(x = fakedate, y = SSTP), colour = "white", size = 0.5, lineend = "round") +
  ## scale_colour_manual(values = (glasbey_mod[1:(nrow(buoyplot))])) +
  ## scale_colour_identity(breaks = buoyplot$col_key) +
  ## xlab(NULL) +
  ## scale_x_date(labels = scales::date_format("%b"), breaks = "2 months") +
  ## theme(legend.position = "none",
  ##       strip.background = element_rect(colour = "grey70", fill = "grey95"))

# Time series plot ####
s <- prevyr %>%
  arrange(STN_ID, fakedate) %>%          # STN_ID %in% c("C46206","C46146")) %>%
  ggplot() +
  facet_wrap(~name_key, ncol = 3) +
  geom_ribbon(data = sstclim, aes(x = fakedate, ymin = SSTP_10perc_day, ymax = SSTP_90perc_day), fill = "grey70", colour = NA, alpha = 0.5) +
  geom_path(data = sstclim, aes(x = fakedate, y = SSTP_clim_mean_day), colour = "grey30", linewidth = 0.8) +
  # Heat dome - for 2021-2022 plots
  # geom_vline(xintercept = as.Date("2021-06-26"), linetype = "dashed") +
  geom_path(data = prevyr %>% filter((year(date) == plot_yrs[1]) | is.na(SSTP)), aes(x = fakedate, y = SSTP, colour = col_key), size = 0.6) +

  geom_path(data = currentyr %>% filter((year(date) == plot_yrs[2]) | is.na(SSTP)),aes(x = fakedate, y = SSTP), colour = "black", size = 1.4, lineend = "round") +
  geom_path(data = currentyr %>% filter(year(date) == plot_yrs[2] | is.na(SSTP)), aes(x = fakedate, y = SSTP), colour = "white", size = 0.5, lineend = "round") +
  scale_colour_manual(values = (glasbey_mod[1:(nrow(buoyplot))])) +
  scale_colour_identity(breaks = buoyplot$col_key) +
  ylab(expression("Mean Daily SST " ( degree*C))) +
  xlab(NULL) +
  scale_x_date(labels = scales::date_format("%b"), breaks = "2 months") +
  theme(legend.position = "none",
        strip.background = element_rect(colour = "grey70", fill = "grey95"))

pbs-assess/PACea documentation built on April 17, 2025, 11:36 p.m.