-
Notifications
You must be signed in to change notification settings - Fork 0
/
buoys.Rmd
172 lines (148 loc) · 7.04 KB
/
buoys.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
---
title: "Buoy Sea Surface Temperature Data"
author: "Andrew Edwards"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Buoy Sea Surface Temperature Data}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
date: "Last rendered on `r format(Sys.time(), '%d %B, %Y')`"
---
<!-- To build either run
rmarkdown::render("buoys.Rmd")
or click the knit button in RStudio -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 5.7,
fig.height = 5
)
```
## Daily average sea surface temperature from buoys
```{r setup}
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
```{r buoymeta}
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
```{r names}
buoy_metadata$name
```
The locations of the buoys are given by
```{r buoysmap, echo=FALSE, out.width = "90%"}
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](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`
```{r buoysst}
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:
```{r buoydates}
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:
```{r plot}
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:
```{r multiloop}
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.
```{r echo=FALSE, eval=FALSE}
# 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"))
```