07-Temperature-at-Depth.qmd

---
title: "Temperature_at_Depth"
author: "phoebe.woodworth-jefcoats@noaa.gov"
format:
  docx:
    reference-doc: SAFE-Reference-Doc.docx
---

## Temperature at 200-300m Depth

```{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: Temperature_at_Depth
# This is where the data are and where the plots will go
Dir <- here("Temperature_at_Depth")
```

```{r}
#| include: false
### Load data
# Monthly temperature
TatD_full <- read_delim(file = paste(Dir, '/T_at_200300_ts.dat', sep = ""), skip = 9, col_names = FALSE)
```

```{r}
#| include: false
# Remove seasonal means to calculate anomalies
TatD_climo <- matrix(NA, nrow = length(TatD_full$X2), ncol = 1)
for (m in seq(1,12,1)) {
  mo_of_int <- which(month(dmy_hm(TatD_full$X1)) == m)
  TatD_climo[mo_of_int,1] <- mean(TatD_full$X2[mo_of_int])
}

TatD_anom_ts <- TatD_full$X2 - TatD_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(TatD_full$X2), 1)
TatD_lm <- lm(TatD_full$X2 ~ n_obs)
TatD_anom_lm <- lm(TatD_anom_ts ~ n_obs)

# Change over time
delta_TatD_lm <- TatD_lm$fitted.values[length(n_obs)] - TatD_lm$fitted.values[1]

delta_TatD_anom_lm <- TatD_anom_lm$fitted.values[length(n_obs)] - TatD_anom_lm$fitted.values[1]
```

```{r}
#| include: false
### Annual values
yrs <- year(dmy_hm(TatD_full$X1))
TatD_full <- bind_cols(TatD_full, yrs)
TatD_full <- rename(TatD_full, Date_Time = X1)
TatD_full <- rename(TatD_full, TatD_degC = X2)
TatD_full <- rename(TatD_full, Year = ...3)

# Add in anomaly to make things easier down the road
TatD_anom_ts <- as_tibble(TatD_anom_ts)
TatD_full <- bind_cols(TatD_full, TatD_anom_ts)
TatD_full <- rename(TatD_full, Anom = V1)

ann_TatD <- TatD_full %>%
  group_by(Year) %>%
  summarise(TatD_degC = mean(TatD_degC, na.rm = TRUE))
ann_anom <- TatD_full %>%
  group_by(Year) %>%
  summarise(Anom = mean(Anom, na.rm = TRUE))
ann_mean_RptYr <- ann_TatD$TatD_degC[which(ann_TatD$Year == RptYr)]
ann_anom_RptYr <- ann_anom$Anom[which(ann_TatD$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(TatD_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(TatD_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(TatD_full$Year == RptYr)
prev_yrs <- which(TatD_full$Year < RptYr)
all_yrs <- which(TatD_full$Year <= RptYr)
monthly_max_RptYr <- max(TatD_full$TatD_degC[yr_of_int])
monthly_min_RptYr <- min(TatD_full$TatD_degC[yr_of_int])
monthly_max_PrevYrs <- max(TatD_full$TatD_degC[prev_yrs])
monthly_min_PrevYrs <- min(TatD_full$TatD_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_TatD <- rename(ann_TatD, `Degrees C` = TatD_degC)
write_csv(ann_TatD, file = paste(here(), '/PelagicClimate_', RptYr, '/Tat200to300m_', RptYr, '.csv', sep = ""))

# Anomaly time series not included given that it tracks actual time series so closely.
```

```{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(TatD_full$Date_Time))
ID <- rep('TatD', dim(TatD_full)[1])
Units <- rep('Deg C', dim(TatD_full)[1])
Value_Anom <- rep(NA, dim(TatD_full)[1]) # No anomaly for this variable
Value_Anom_lm <- rep(NA, dim(TatD_full)[1]) # No anomaly for this variable

TatD_dashboard <- bind_cols(TatD_full$Date_Time, 
                           TatD_full$Year, 
                           Month,
                           TatD_full$TatD_degC,
                           Value_Anom, # TatD_full$Anom,
                           ID,
                           TatD_lm$fitted.values,
                           Value_Anom_lm, # TatD_anom_lm$fitted.values,
                           Units)
# Need to figure out how to render this unnecessary
TatD_dashboard <- rename(TatD_dashboard, Date_Time = ...1)
TatD_dashboard <- rename(TatD_dashboard, Year = ...2)
TatD_dashboard <- rename(TatD_dashboard, Month = ...3)
TatD_dashboard <- rename(TatD_dashboard, Value = ...4)
TatD_dashboard <- rename(TatD_dashboard, Anom = ...5)
TatD_dashboard <- rename(TatD_dashboard, ID = ...6)
TatD_dashboard <- rename(TatD_dashboard, Value_lm = ...7)
TatD_dashboard <- rename(TatD_dashboard, Anom_lm = ...8)
TatD_dashboard <- rename(TatD_dashboard, Units = ...9)

write_csv(TatD_dashboard, file = here("Indicator_Dashboard", "Data", paste('TatD_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 <- TatD_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$`Anom`[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(10.8, 11.7), tickvals = list(10.8, 10.9, 11, 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7)), #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, '/TempAtDepth_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'), #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, '/TempAtDepth_anom_ts_', RptYr, '.pdf', sep = ""))
```

Rationale: The temperature at 200–300 m reflects the temperature in the mid-range of depths targeted by the deep-set bigeye tuna fishery. Bigeye have preferred thermal habitat, generally staying within temperatures ranging from 8–14 °C while they are at depth (Howell et al. 2010). Changes in ocean temperature at depth will impact tuna, and in turn, potentially impact their catchability. Understanding the drivers of sub-surface temperature trends and their ecosystem impacts is an area of active research.  
\
Status: In `r RptYr`, 200–300 m temperatures ranged from `r signif(monthly_min_RptYr, 4)`–`r signif(monthly_max_RptYr, 4)` °C with an average value of `r signif(ann_mean_RptYr,4)` °C. These temperatures are within the range of temperatures experienced over the past several decades (`r signif(monthly_min_PrevYrs, 4)`–`r signif(monthly_max_PrevYrs, 4)` °C) and are within the bounds of bigeye tuna’s preferred deep daytime thermal habitat (8–14 °C). Over the period of record (1980–`r RptYr`), 200–300 m temperatures have declined by `r abs(signif(delta_TatD_lm, 1))` °C. The spatial pattern of temperature anomalies was mixed with temperatures at depth around the main Hawaiian Islands roughly 0.5–1.5 °C below average, and temperatures north of about 30°N 0–0.5 °C above average.   
\
Description: Ocean temperature at 200–300 m depth is averaged across the Hawaiʻi-based longline fishing grounds (15° – 45°N, 180° – 120°W). Global Ocean Data Assimilation System (GODAS) data are used. GODAS incorporates global ocean data from moorings, expendable bathythermographs (XBTs), and Argo floats.   
\
Timeframe: Annual, monthly.  
\
Region/Location: Hawaii longline region: 15° – 45°N, 180° – 120°W.  
\
Measurement Platform: _In-situ_ sensors, model.  
\
Sourced from: NOAA (2024d) and APDRC (2024).  Graphics produced in part using Stawitz (2023).

## Additional Information
Temperature-at-depth data are geographically subset and spatially averaged using the pyFerret script `TatDepth_access.jnl` which can be found in the [Temperature_at_Depth](https://github.com/pwoodworth-jefcoats/climate-indicators/tree/main/Temperature_at_Depth) 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 showed that they were evenly distributed, although were more positive (~0.4 max) than negative (~-0.3 min). The residuals for the anomaly model were also fairly evenly distributed in terms of values.  Both sets of 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, '/TempAtDepth_map_data_', RptYr, '.csv', sep = ""))
# 
# # Filter to the given year and the anomaly
# maps_RptYr <- maps |> filter(ID == "TatD")
# maps_anom <- maps |> filter(ID == "TatD_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 = "TatD", 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)), -0.5, 0.5)),
#                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(0.5, ceiling(max(maps_anom$layer)), 0.5)),
#                color = "black", lwd = 0.5) # solid positive contours
# 
# p
# 
# pdf(paste(Dir, '/TempAtDepth_map_', RptYr, '.pdf', sep = ""))
# print(p)
# dev.off()
```