-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path08-Ocean-Color.qmd
357 lines (313 loc) · 16.7 KB
/
08-Ocean-Color.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
---
title: "Ocean_Color"
author: "[email protected]"
format:
docx:
reference-doc: SAFE-Reference-Doc.docx
---
## Ocean Color
```{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)
library(scales)
reticulate::use_miniconda('r-reticulate')
```
```{r}
#| include: false
# Set report year (RptYr), to make things easier
RptYr <- 2023
# Set path to variable: Ocean_Color
# This is where the data are and where the plots will go
Dir <- here("Ocean_Color")
```
```{r}
#| include: false
### Load data
# Monthly ocean color
# Generated in ~5-yr batches due to memory limitations
Chl_9802 <- read_delim(file = paste(Dir, '/chl_ts_9802.dat', sep = ""), skip = 8, col_names = FALSE)
Chl_0307 <- read_delim(file = paste(Dir, '/chl_ts_0307.dat', sep = ""), skip = 8, col_names = FALSE)
Chl_0812 <- read_delim(file = paste(Dir, '/chl_ts_0812.dat', sep = ""), skip = 8, col_names = FALSE)
Chl_1317 <- read_delim(file = paste(Dir, '/chl_ts_1317.dat', sep = ""), skip = 8, col_names = FALSE)
Chl_1822 <- read_delim(file = paste(Dir, '/chl_ts_1822.dat', sep = ""), skip = 8, col_names = FALSE)
Chl_23 <- read_delim(file = paste(Dir, '/chl_ts_23.dat', sep = ""), skip = 8, col_names = FALSE)
# Concatonate
Chl_full <- rbind(Chl_9802, Chl_0307, Chl_0812, Chl_1317, Chl_1822, Chl_23)
```
```{r}
#| include: false
# Remove seasonal means to calculate anomalies
Chl_climo <- matrix(NA, nrow = length(Chl_full$X2), ncol = 1)
for (m in seq(1,12,1)) {
mo_of_int <- which(month(dmy_hm(Chl_full$X1)) == m)
Chl_climo[mo_of_int,1] <- mean(Chl_full$X2[mo_of_int])
}
Chl_anom_ts <- Chl_full$X2 - Chl_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(Chl_full$X2), 1)
Chl_lm <- lm(Chl_full$X2 ~ n_obs)
Chl_anom_lm <- lm(Chl_anom_ts ~ n_obs)
# Change over time
delta_Chl_lm <- Chl_lm$fitted.values[length(n_obs)] - Chl_lm$fitted.values[1]
delta_Chl_anom_lm <- Chl_anom_lm$fitted.values[length(n_obs)] - Chl_anom_lm$fitted.values[1]
```
```{r}
#| include: false
### Annual values
yrs <- year(dmy_hm(Chl_full$X1))
Chl_full <- bind_cols(Chl_full, yrs)
Chl_full <- rename(Chl_full, Date_Time = X1)
Chl_full <- rename(Chl_full, Chl_mgm3 = X2)
Chl_full <- rename(Chl_full, Year = ...3)
# Add in anomaly to make things easier down the road
Chl_anom_ts <- as_tibble(Chl_anom_ts)
Chl_full <- bind_cols(Chl_full, Chl_anom_ts)
Chl_full <- rename(Chl_full, Anom = V1)
ann_Chl <- Chl_full %>%
group_by(Year) %>%
summarise(Chl_mgm3 = mean(Chl_mgm3, na.rm = TRUE))
ann_anom <- Chl_full %>%
group_by(Year) %>%
summarise(Anom = mean(Anom, na.rm = TRUE))
ann_mean_RptYr <- ann_Chl$Chl_mgm3[which(ann_Chl$Year == RptYr)]
ann_anom_RptYr <- ann_anom$Anom[which(ann_Chl$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(Chl_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(Chl_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(Chl_full$Year == RptYr)
prev_yrs <- which(Chl_full$Year < RptYr)
all_yrs <- which(Chl_full$Year <= RptYr)
monthly_max_RptYr <- max(Chl_full$Chl_mgm3[yr_of_int])
monthly_min_RptYr <- min(Chl_full$Chl_mgm3[yr_of_int])
monthly_max_PrevYrs <- max(Chl_full$Chl_mgm3[prev_yrs])
monthly_min_PrevYrs <- min(Chl_full$Chl_mgm3[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_Chl <- rename(ann_Chl, `mg/m^3` = Chl_mgm3)
write_csv(ann_Chl, file = paste(here(), '/PelagicClimate_', RptYr, '/CHL_', RptYr, '.csv', sep = ""))
ann_anom <- rename(ann_anom, `mg/m^3` = Anom)
write_csv(ann_anom, file = paste(here(), '/PelagicClimate_', RptYr, '/CHLanomaly_', 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(Chl_full$Date_Time))
ID <- rep('Chl', dim(Chl_full)[1])
Units <- rep('mg chl per m3', dim(Chl_full)[1])
Value_lm <- rep(NA, dim(Chl_full)[1]) # No significant trend in chl
Chl_dashboard <- bind_cols(Chl_full$Date_Time,
Chl_full$Year,
Month,
Chl_full$Chl_mgm3,
Chl_full$Anom,
ID,
Value_lm, # Chl_lm$fitted.values,
Chl_anom_lm$fitted.values,
Units)
# Need to figure out how to render this unnecessary
Chl_dashboard <- rename(Chl_dashboard, Date_Time = ...1)
Chl_dashboard <- rename(Chl_dashboard, Year = ...2)
Chl_dashboard <- rename(Chl_dashboard, Month = ...3)
Chl_dashboard <- rename(Chl_dashboard, Value = ...4)
Chl_dashboard <- rename(Chl_dashboard, Anom = ...5)
Chl_dashboard <- rename(Chl_dashboard, ID = ...6)
Chl_dashboard <- rename(Chl_dashboard, Value_lm = ...7)
Chl_dashboard <- rename(Chl_dashboard, Anom_lm = ...8)
Chl_dashboard <- rename(Chl_dashboard, Units = ...9)
write_csv(Chl_dashboard, file = here("Indicator_Dashboard", "Data", paste('Chl_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 <- Chl_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$`mg/m^3`[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(0.07, 0.19), tickvals = list(0.07, 0.08, 0.1, 0.12, 0.14, 0.16, 0.18, 0.19)), #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, '/Chl_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(-0.035, 0.035), tickvals = list(-0.035, -0.025, -0.015, -0.005, 0, 0.005, 0.015, 0.025, 0.035)), #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, '/Chl_anom_ts_', RptYr, '.pdf', sep = ""))
```
Rationale: Phytoplankton are the foundational food source for the fishery. Changes in phytoplankton abundance have been linked to both natural climate variability and anthropogenic climate change. These changes have the potential to impact fish abundance, size, and catch.
\
Status: The mean monthly chlorophyll concentration was `r signif(ann_mean_RptYr, 2)` mg chl m^-^^3^ in `r RptYr`. Monthly mean chlorophyll concentrations ranged from `r signif(monthly_min_RptYr, 2)`–`r signif(monthly_max_RptYr, 2)` mg chl m^-^^3^, which was within the range of values observed during the previous years of the time series (`r signif(monthly_min_PrevYrs, 2)`–`r signif(monthly_max_PrevYrs, 2)` mg chl m^-^^3^). There has been no significant trend in monthly average chlorophyll concentration over the time period, however chlorophyll anomalies have declined by `r abs(signif(delta_Chl_anom_lm, 1))` mg chl m^-^^3^. Chlorophyll concentrations were fairly average across the southern portion of the longline fishing grounds and a little below average north of 30–35°N.
\
Description: Satellite remotely sensed ocean color is used to determine chlorophyll concentrations in the pelagic surface ocean. A time series of mean monthly chlorophyll-a concentrations averaged over the Hawaiʻi longline region is presented. Additionally, spatial climatologies and anomalies are shown. European Space Agency (ESA) Climate Change Initiative (CCI) data are used for this indicator (Sathyendranath et al. 2018).
\
Timeframe: Monthly
\
Region/Location: Hawaii longline region: 5° – 45°N, 180° – 120°W
\
Measurement Platform: Satellite
\
Data available at: <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/esa-cci-chla-monthly-v6-0>, <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/esa-cci-chla-1998-2009-clim-v6-0>, and <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/esa-cci-chla-2023-clim_v6-0>.
\
Sourced from: NOAA OceanWatch (2024b) and Sathyendranath et al. (2018). Graphics produced in part using Stawitz (2023).
\
## Additional Information
Ocean color data are geographically subset and spatially averaged using the pyFerret script `CHLaccess.jnl` which can be found in the Ocean_Color folder. 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 for the anomalies showed that they were evenly distributed, although were more positive (~0.03 max) than negative (~-0.02 min). These residual exhibited periodicity that closely track the time series pattern.
\
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, '/Chl_map_data_', RptYr, '.csv', sep = ""))
#
# # Filter to the given year and the anomaly
# maps_RptYr <- maps |> filter(ID == "Chl")
# maps_anom <- maps |> filter(ID == "Chl_anom")
#
# # Map elements and aesthetics
# seagrass <- nmfs_palette("seagrass")(3)
# pal <- seagrass
# ll_rect_color <- "white" #outline for ll fishing grounds box
# fill_scale <- scale_fill_gradientn(name = "Chl", trans = log10_trans(),
# colors = pal, limits = c(0.01, 10), oob = squish) #log scale
#
# # 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(-0.02),
# color = "black", linetype = 3) # Dotted negative contour
# 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(0.02),
# color = "black", lwd = 0.5) # solid positive contour
#
# p
#
# pdf(paste(Dir, '/Chl_map_', RptYr, '.pdf', sep = ""))
# print(p)
# dev.off()
```