-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path06-Sea-Surface-Temperature.qmd
360 lines (315 loc) · 17.2 KB
/
06-Sea-Surface-Temperature.qmd
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
---
title: "Sea_Surface_Temperature"
author: "[email protected]"
format:
docx:
reference-doc: SAFE-Reference-Doc.docx
---
## Sea Surface Temperature
```{r}
#| include: false
### Load libraries
library(tidyverse)
library(lubridate)
library(here)
library(stringr)
# remotes::install_github("nmfs-fish-tools/nmfspalette")
library(nmfspalette)
library(plotly)
library(reticulate)
reticulate::use_miniconda('r-reticulate')
```
```{r}
#| include: false
# Set report year (RptYr), to make things easier
RptYr <- 2023
# Set path to variable: Sea_Surface_Temperature
# This is where the data are and where the plots will go
Dir <- here("Sea_Surface_Temperature")
```
```{r}
#| include: false
### Load data
# Monthly SST
# Generated in ~5-yr batches due to memory limitations
SST_8590 <- read_delim(file = paste(Dir, '/sst_ts_8590.dat', sep = ""), skip = 8, col_names = FALSE)
SST_9195 <- read_delim(file = paste(Dir, '/sst_ts_9195.dat', sep = ""), skip = 8, col_names = FALSE)
SST_9600 <- read_delim(file = paste(Dir, '/sst_ts_9600.dat', sep = ""), skip = 8, col_names = FALSE)
SST_0105 <- read_delim(file = paste(Dir, '/sst_ts_0105.dat', sep = ""), skip = 8, col_names = FALSE)
SST_0610 <- read_delim(file = paste(Dir, '/sst_ts_0610.dat', sep = ""), skip = 8, col_names = FALSE)
SST_1115 <- read_delim(file = paste(Dir, '/sst_ts_1115.dat', sep = ""), skip = 8, col_names = FALSE)
SST_1620 <- read_delim(file = paste(Dir, '/sst_ts_1620.dat', sep = ""), skip = 8, col_names = FALSE)
SST_2123 <- read_delim(file = paste(Dir, '/sst_ts_2123.dat', sep = ""), skip = 8, col_names = FALSE)
# Concatonate
SST_full <- rbind(SST_8590, SST_9195, SST_9600, SST_0105, SST_0610, SST_1115, SST_1620, SST_2123)
```
```{r}
#| include: false
# Remove seasonal means to calculate anomalies
SST_climo <- matrix(NA, nrow = length(SST_full$X2), ncol = 1)
for (m in seq(1,12,1)) {
mo_of_int <- which(month(dmy_hm(SST_full$X1)) == m)
SST_climo[mo_of_int,1] <- mean(SST_full$X2[mo_of_int])
}
SST_anom_ts <- SST_full$X2 - SST_climo
```
```{r}
#| include: false
### Linear fit
# Note that this assumes that observations are equally spaced in time, which they're not
n_obs <- seq(1, length(SST_full$X2), 1)
SST_lm <- lm(SST_full$X2 ~ n_obs)
SST_anom_lm <- lm(SST_anom_ts ~ n_obs)
# Change over time
delta_SST_lm <- SST_lm$fitted.values[length(n_obs)] - SST_lm$fitted.values[1]
delta_SST_anom_lm <- SST_anom_lm$fitted.values[length(n_obs)] -
SST_anom_lm$fitted.values[1]
```
```{r}
#| include: false
### Annual values
yrs <- year(dmy_hm(SST_full$X1))
SST_full <- bind_cols(SST_full, yrs)
SST_full <- rename(SST_full, Date_Time = X1)
SST_full <- rename(SST_full, SST_degC = X2)
SST_full <- rename(SST_full, Year = ...3)
# Add in anomaly to make things easier down the road
SST_anom_ts <- as_tibble(SST_anom_ts)
SST_full <- bind_cols(SST_full, SST_anom_ts)
SST_full <- rename(SST_full, Anom = V1)
ann_SST <- SST_full %>%
group_by(Year) %>%
summarise(SST_degC = mean(SST_degC, na.rm = TRUE))
ann_anom <- SST_full %>%
group_by(Year) %>%
summarise(Anom = mean(Anom, na.rm = TRUE))
ann_mean_RptYr <- ann_SST$SST_degC[which(ann_SST$Year == RptYr)]
ann_anom_RptYr <- ann_anom$Anom[which(ann_SST$Year == RptYr)]
```
```{r}
#| echo: false
# Note that the above needs to be 'echo' and not 'include' so that the error checks print.
# This section includes some error checks to prompt fixing the text
if (any(summary(SST_lm)$coefficients[,4] > 0.05)) {
print('The linear fit to monthly values is not signficant. Remove text related to the linear trend.')
}
if (any(summary(SST_anom_lm)$coefficients[,4] > 0.05)) {
print('The linear fit to anomaly values is not signficant. Remove text related to the linear trend.')
}
# Pull out the MONTHLY values we need for the report & plots
yr_of_int <- which(SST_full$Year == RptYr)
prev_yrs <- which(SST_full$Year < RptYr)
all_yrs <- which(SST_full$Year <= RptYr)
monthly_max_RptYr <- max(SST_full$SST_degC[yr_of_int])
monthly_min_RptYr <- min(SST_full$SST_degC[yr_of_int])
monthly_max_PrevYrs <- max(SST_full$SST_degC[prev_yrs])
monthly_min_PrevYrs <- min(SST_full$SST_degC[prev_yrs])
if (monthly_max_RptYr > monthly_max_PrevYrs) {
print('The greatest monthly value was during the report year. Revise text to reflect this.')
}
if (monthly_min_RptYr < monthly_min_PrevYrs) {
print('The lowest monthly value was during the report year. Revise text to reflect this.')
}
```
```{r}
#| include: false
# Write csv for portal
# Note that output csvs go in their own folder
ann_SST <- rename(ann_SST, `Degrees C` = SST_degC)
write_csv(ann_SST, file = paste(here(), '/PelagicClimate_', RptYr, '/SST_', RptYr, '.csv', sep = ""))
ann_anom <- rename(ann_anom, `Degrees C` = Anom)
write_csv(ann_anom, file = paste(here(), '/PelagicClimate_', RptYr, '/SSTanomaly_', RptYr, '.csv', sep = ""))
```
```{r}
#| include: false
# Write csv for dashboard
# Note that dashboard output has its own folder
# Add columns for month, variable ID, and units
Month <- month(dmy_hm(SST_full$Date_Time))
ID <- rep('SST', dim(SST_full)[1])
Units <- rep('Deg C', dim(SST_full)[1])
SST_dashboard <- bind_cols(SST_full$Date_Time,
SST_full$Year,
Month,
SST_full$SST_degC,
SST_full$Anom,
ID,
SST_lm$fitted.values,
SST_anom_lm$fitted.values,
Units)
# Need to figure out how to render this unnecessary
SST_dashboard <- rename(SST_dashboard, Date_Time = ...1)
SST_dashboard <- rename(SST_dashboard, Year = ...2)
SST_dashboard <- rename(SST_dashboard, Month = ...3)
SST_dashboard <- rename(SST_dashboard, Value = ...4)
SST_dashboard <- rename(SST_dashboard, Anom = ...5)
SST_dashboard <- rename(SST_dashboard, ID = ...6)
SST_dashboard <- rename(SST_dashboard, Value_lm = ...7)
SST_dashboard <- rename(SST_dashboard, Anom_lm = ...8)
SST_dashboard <- rename(SST_dashboard, Units = ...9)
write_csv(SST_dashboard, file = here("Indicator_Dashboard", "Data", paste('SST_Dashboard_Data_', RptYr, '.csv', sep = "")))
```
```{r}
#| include: false
# Borrowing code from the dashboard for this chunk
# so that figures look the same across products
indicator_data <- SST_dashboard |>
filter(Year <= RptYr)
# Create color palette for easy reference
oceans <- nmfs_palette("oceans")(3) # 1 = report_year, 3 = previous years
crustacean <- nmfs_palette("crustacean")(4) # 1 = linear trend
coral <- nmfs_palette("coral")(3) # 3 = annual average for Rpt Yr
ann_grey <- "#D0D0D0" # annual means; in NMFS branding guide but not in package
waves <- nmfs_palette("waves")(3) # annual means; in NMFS branding guide but not in package
seagrass <- nmfs_palette("seagrass")(3)
pal <- c(oceans[3], coral[2], waves[2], coral[3], crustacean[2])
# Formatting
plot_title_font <- list(size = 14)
plot_height <- 350 #in pixels
# Calculate annual means
ann_vals <- indicator_data |>
group_by(Year) |>
summarise(Value = mean(Value, na.rm = TRUE))
# Identify the current year, to overlay on plot
given_yr <- indicator_data |>
filter(Year == RptYr)
given_yr_ann <- bind_cols(rep(ann_vals$Value[dim(ann_vals)[1]]),
given_yr$Date_Time,
ann_anom$`Degrees C`[which(ann_anom$Year == RptYr)])
given_yr_ann <- rename(given_yr_ann, Value = ...1)
given_yr_ann <- rename(given_yr_ann, Date_Time = ...2)
given_yr_ann <- rename(given_yr_ann, Anom = ...3)
p1 <- plot_ly(indicator_data, x = dmy_hm(indicator_data$Date_Time), y = ~Value,
type = "scatter", mode = "lines", line = list(color = pal[1]),
name = ~ID[1], height = plot_height) |>
add_trace(indicator_data, x = dmy_hm(indicator_data$Date_Time), y = ~Value_lm,
type = "scatter", mode = "lines", line = list(color = pal[5]),
name = "Long-term Trend") |>
add_trace(ann_vals, x = ymd_hm(paste(ann_vals$Year, '0601 00:00', sep = "")), y = ann_vals$Value,
type = "scatter", mode = "lines", line = list(color = pal[3]),
name = "Annual Mean") |>
add_trace(given_yr, x = dmy_hm(given_yr$Date_Time), y = given_yr$Value,
type = "scatter", mode = "lines", line = list(color = pal[2]),
name = ~ID[1], height = plot_height) |>
add_trace(given_yr_ann, x = dmy_hm(given_yr_ann$Date_Time), y = given_yr_ann$Value,
type = "scatter", mode = "lines", line = list(color = pal[4]),
name = ~ID[1], height = plot_height)
#apply same layout parameters for all plots
#custom x axis (min, every decade, report year)
all_years <- unique(indicator_data$Year)
first_date <- as.character(parse_date_time(indicator_data$Date_Time[1], orders = "d-b-Y H:M"))
date_axis <- c(first_date,
all_years[which(all_years %% 10 == 0)],
RptYr)
p1 <- p1 |> layout(title = list(text = "Indicator Time Series", x = 0.01, font = plot_title_font), #add title
xaxis = list(type = "date", tickformat = "%Y", tickmode = "array", tickvals = date_axis, tickangle = 90),
#xaxis = list(tick0 = min(indicator_data$Date_Time), dtick = "M24"),
yaxis = list(title = indicator_data$Units[1], hoverformat = '.3f', range = list(17, 25), tickvals = list(17, 18, 19, 20, 21, 22, 23, 24, 25)), #add units and round values in hover display; not sure of a better way to set axis limits and ticks...
paper_bgcolor = 'transparent', plot_bgcolor = 'transparent', #transparent background
hovermode = "x unified", #show data for all traces simultaneously
hoverlabel = list(namelength = -1)) #don't cutoff hoverinfo due to length
# return plot
save_image(p1, paste(Dir, '/SST_ts_', RptYr, '.pdf', sep = ""))
### Anomaly plot
# Calculate annual means
ann_vals <- indicator_data |>
group_by(Year) |>
summarise(Anom = mean(Anom, na.rm = TRUE))
p2 <- plot_ly(indicator_data, x = dmy_hm(indicator_data$Date_Time), y = ~Anom, height = plot_height,
type = "scatter", mode = "lines", line = list(color = pal[1]),
name = "Monthly Anomaly") |>
add_trace(indicator_data, x = dmy_hm(indicator_data$Date_Time), y = ~Anom_lm,
type = "scatter", mode = "lines", line = list(color = pal[5]),
name = "Long-term Trend") |>
add_trace(ann_vals, x = ymd_hm(paste(ann_vals$Year, '0601 00:00', sep = "")), y = ann_vals$Anom,
type = "scatter", mode = "lines", line = list(color = pal[3]),
name = "Annual Mean") |>
add_trace(given_yr, x = dmy_hm(given_yr$Date_Time), y = given_yr$Anom,
type = "scatter", mode = "lines", line = list(color = pal[2]),
name = ~ID[1], height = plot_height) |>
add_trace(given_yr_ann, x = dmy_hm(given_yr_ann$Date_Time), y = given_yr_ann$Anom,
type = "scatter", mode = "lines", line = list(color = pal[4]),
name = ~ID[1], height = plot_height) |>
layout(xaxis = list(type = "date", tickformat = "%Y", tickmode = "array", tickvals = date_axis, tickangle = 90),
yaxis = list(title = indicator_data$Units[1], hoverformat = '.3f', range = list(-1.5, 1.5), tickvals = list(-1.5, -1, -0.5, 0, 0.5, 1.0, 1.5)), #add units and round values in hover display
title = list(text = "Anomaly Time Series", x = 0.01, font = plot_title_font), #add title
paper_bgcolor = 'transparent', plot_bgcolor = 'transparent', #transparent background
hovermode = "x unified", #show data for all traces simultaneously
hoverlabel = list(namelength = -1)) #don't cutoff hoverinfo due to length
# return plot
save_image(p2, paste(Dir, '/SST_anom_ts_', RptYr, '.pdf', sep = ""))
```
Rationale: Sea surface temperature is one of the most directly observable existing measures for tracking increasing ocean temperatures. SST varies in response to natural climate cycles such as the El Niño – Southern Oscillation (ENSO) and is rising as a result of anthropogenic climate change. Both short-term variability and long-term trends in SST impact the marine ecosystem. Understanding the mechanisms through which organisms are impacted and the time scales of these impacts is an area of active research.
\
Status: Annual mean SST was `r signif(ann_mean_RptYr,3)` ºC in `r RptYr`. Over the period of record, SST across the longline fishing grounds has increased by `r signif(delta_SST_lm, 1)` ºC and the monthly SST anomaly increased by `r signif(delta_SST_anom_lm, 1)` ºC, both at a rate of roughly 0.03 ºC yr^-^^1^. Monthly SST values in `r RptYr` ranged from `r signif(monthly_min_RptYr, 3)`–`r signif(monthly_max_RptYr, 3)` ºC, within the range of temperatures experienced over the past several decades (`r signif(monthly_min_PrevYrs, 3)`–`r signif(monthly_max_PrevYrs, 3)` ºC). Overall, SST was above the long-term average across most of the Hawaiʻi longline region in `r RptYr`. The exception to this was a patch of slightly cooler waters in the southeastern corner of the fishing grounds where very little fishing takes place.
\
Description: Satellite remotely sensed monthly sea surface temperature (SST) is averaged across the Hawaiʻi-based longline fishing grounds (15° – 45°N, 180° – 120°W). A time series of monthly mean SST averaged over the Hawaiʻi longline region is presented. Additionally, spatial climatologies and anomalies are shown. CoralTemp data are used to calculate this indicator.
\
Timeframe: Monthly.
\
Region/Location: Hawaiʻi longline region: 15° – 45°N, 180° – 120°W.
\
Measurement Platform: Satellite.
\
Data available at: <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_monthly>, <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_1985-2009-clim>, and <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v3_1_2023-clim>.
\
Sourced from: NOAA OceanWatch (2024a). Graphics produced in part using Stawitz (2023).
\
## Additional Information
SST data are geographically subset and spatially averaged using the pyFerret script `SSTaccess.jnl` which can be found in the [Sea_Surface_Temperature folder](https://github.com/pwoodworth-jefcoats/climate-indicators/tree/main/Sea_Surface_Temperature). You could do this step in R, but I use [PyFerret](https://ferret.pmel.noaa.gov/Ferret/) because it's freely available software developed by NOAA specifically for working with large gridded datasets.
\
A plot of the residuals from the linear model showed that they were evenly distributed, although were more positive (~4.0 max) than negative (~-3.0 min). The residuals for the anomaly model were also fairly evenly distributed in terms of values, however they did appear to exhibit some periodicity.
\
To prepare the spatial data for mapping, you'll need to run `map_data_for_dashboard.R`, which can be found in the [PreProcessing](https://github.com/pwoodworth-jefcoats/climate-indicators/tree/main/Indicator_Dashboard/PreProcessing) folder of the [Indicator_Dashboard](https://github.com/pwoodworth-jefcoats/climate-indicators/tree/main/Indicator_Dashboard) folder in this repository.
```{r}
#| include: FALSE
# # After running the chunks above and prepping the map data, you can uncomment and run this chunk.
#
# # Load basemap data
# land_proj <- readRDS(here('Indicator_Dashboard', 'Data', 'rnatearth_land.RData'))
# coast_proj <- readRDS(here('Indicator_Dashboard', 'Data', 'rnatearth_coast.RData'))
#
# # Load data
# maps <- read.csv(paste(Dir, '/SST_map_data_', RptYr, '.csv', sep = ""))
#
# # Filter to the given year and the anomaly
# maps_RptYr <- maps |> filter(ID == "SST")
# maps_anom <- maps |> filter(ID == "SST_anom")
#
# # Map elements and aesthetics
# waves <- nmfs_palette("waves")(3)
# pal <- rev(waves)
# ll_rect_color <- "white" #outline for ll fishing grounds box
# fill_scale <- scale_fill_gradientn(name = "SST", colors = pal, limits = c(0, 31))
#
# # Create map
# p <- ggplot() +
# geom_raster(data = maps_RptYr, mapping = aes(x = x, y = y, fill = layer)) + #data as raster
# geom_rect(aes(xmin = 0, xmax = 60, ymin = 15, ymax = 45), color = ll_rect_color, fill = NA, linewidth = 0.5) + #ll fishing grounds box
# annotate("text", x = 30, y = 47, label = "Longline fishing grounds", size = 3.2, color = ll_rect_color) +
# fill_scale + #indicator-dependent color and scale from above
# geom_sf(data = land_proj, fill = "#A5A5A5", color = "#A5A5A5", linewidth = 0.5) + #base map background
# geom_sf(data = coast_proj) + #base map outline
# coord_sf(expand = F) + #don't expand past x/y limits
# xlab("") + ylab("") +
# theme_bw() + theme(panel.grid = element_line(color = "black", linewidth = 0.1),
# legend.background = element_rect(fill = 'transparent'))
# # Add anomaly
# p <- p +
# geom_contour(data = maps_anom, mapping = aes(x = x, y = y, z = layer),
# breaks = c(seq(floor(min(maps_anom$layer)), -1, 1)),
# color = "black", linetype = 3) # Dotted negative contours
# p <- p +
# geom_contour(data = maps_anom, mapping = aes(x = x, y = y, z = layer), breaks = c(0),
# color = "black", lwd = 1) # heavy zero line
# p <- p +
# geom_contour(data = maps_anom, mapping = aes(x = x, y = y, z = layer),
# breaks = c(seq(1, ceiling(max(maps_anom$layer)), 1)),
# color = "black", lwd = 0.5) # solid positive contours
#
# p
#
# pdf(paste(Dir, '/SST_map_', RptYr, '.pdf', sep = ""))
# print(p)
# dev.off()
```