10-Median-Phytoplankton-Size.qmd

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

## Estimated Median Phytoplankton Size

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

```{r}
#| include: false
### Load data
# Monthly estimated median phytoplankton size
# Generated in ~5-yr batches due to memory limitations
med_phyto_9802 <- read_delim(file = paste(Dir, '/medphyto_ts_9802.dat', sep = ""), skip = 8, col_names = FALSE)
med_phyto_0307 <- read_delim(file = paste(Dir, '/medphyto_ts_0307.dat', sep = ""), skip = 8, col_names = FALSE)
med_phyto_0812 <- read_delim(file = paste(Dir, '/medphyto_ts_0812.dat', sep = ""), skip = 8, col_names = FALSE)
med_phyto_1317 <- read_delim(file = paste(Dir, '/medphyto_ts_1317.dat', sep = ""), skip = 8, col_names = FALSE)
med_phyto_1822 <- read_delim(file = paste(Dir, '/medphyto_ts_1822.dat', sep = ""), skip = 8, col_names = FALSE)
med_phyto_23 <- read_delim(file = paste(Dir, '/medphyto_ts_23.dat', sep = ""), skip = 8, col_names = FALSE)


# Concatonate
med_phyto_full <- rbind(med_phyto_9802, 
                        med_phyto_0307, 
                        med_phyto_0812, 
                        med_phyto_1317, 
                        med_phyto_1822,
                        med_phyto_23)
```

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

med_phyto_anom_ts <- med_phyto_full$X2 - medphyto_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(med_phyto_full$X2), 1)
medphyto_lm <- lm(med_phyto_full$X2 ~ n_obs)
medphyto_anom_lm <- lm(med_phyto_anom_ts ~ n_obs)

# Change over time
delta_size_lm <- medphyto_lm$fitted.values[length(n_obs)] - medphyto_lm$fitted.values[1]
delta_size_pct = (medphyto_lm$fitted.values[length(n_obs)] - medphyto_lm$fitted.values[1]) / 
  medphyto_lm$fitted.values[1] * 100

delta_size_anom_lm <- medphyto_anom_lm$fitted.values[length(n_obs)] - 
  medphyto_anom_lm$fitted.values[1]
```

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

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

ann_size <- med_phyto_full %>%
  group_by(Year) %>%
  summarise(ESD_um = mean(ESD_um, na.rm = TRUE))
ann_anom <- med_phyto_full %>%
  group_by(Year) %>%
  summarise(Anom = mean(Anom, na.rm = TRUE))
ann_mean_RptYr <- ann_size$ESD_um[which(ann_size$Year == RptYr)]
ann_anom_RptYr <- ann_anom$Anom[which(ann_size$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(medphyto_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(medphyto_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(med_phyto_full$Year == RptYr)
prev_yrs <- which(med_phyto_full$Year < RptYr)
all_yrs <- which(med_phyto_full$Year <= RptYr)
monthly_max_RptYr <- max(med_phyto_full$ESD_um[yr_of_int])
monthly_min_RptYr <- min(med_phyto_full$ESD_um[yr_of_int])
monthly_max_PrevYrs <- max(med_phyto_full$ESD_um[prev_yrs])
monthly_min_PrevYrs <- min(med_phyto_full$ESD_um[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_size <- rename(ann_size, `ESD in um` = ESD_um)
write_csv(ann_size, file = paste(here(), '/PelagicClimate_', RptYr, '/MedianPhytoSize_', RptYr, '.csv', sep = ""))

ann_anom <- rename(ann_anom, `ESD in um` = Anom)
write_csv(ann_anom, file = paste(here(), '/PelagicClimate_', RptYr, '/MedianPhytoSizeAnomaly_', 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(med_phyto_full$Date_Time))
ID <- rep('MD50', dim(med_phyto_full)[1])
Units <- rep('ESD in um', dim(med_phyto_full)[1])

MD50_dashboard <- bind_cols(med_phyto_full$Date_Time, 
                           med_phyto_full$Year, 
                           Month,
                           med_phyto_full$ESD_um,
                           med_phyto_full$Anom,
                           ID,
                           medphyto_lm$fitted.values,
                           medphyto_anom_lm$fitted.values,
                           Units)
# Need to figure out how to render this unnecessary
MD50_dashboard <- rename(MD50_dashboard, Date_Time = ...1)
MD50_dashboard <- rename(MD50_dashboard, Year = ...2)
MD50_dashboard <- rename(MD50_dashboard, Month = ...3)
MD50_dashboard <- rename(MD50_dashboard, Value = ...4)
MD50_dashboard <- rename(MD50_dashboard, Anom = ...5)
MD50_dashboard <- rename(MD50_dashboard, ID = ...6)
MD50_dashboard <- rename(MD50_dashboard, Value_lm = ...7)
MD50_dashboard <- rename(MD50_dashboard, Anom_lm = ...8)
MD50_dashboard <- rename(MD50_dashboard, Units = ...9)

write_csv(MD50_dashboard, file = here("Indicator_Dashboard", "Data", paste('MD50_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 <- MD50_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$`ESD in um`[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(1.1, 1.9), tickvals = list(1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9)), #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, '/Median_Phyto_Size_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.15, 0.15), tickvals = list(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)), #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, '/Median_Phyto_Size_anom_ts_', RptYr, '.pdf', sep = ""))
```

Rationale: Phytoplankton are the base of the food web and their abundance influences the food available to all higher trophic levels from zooplankton through tuna and billfish. Some studies project that climate change will result in both fewer and smaller phytoplankton. This would reduce the food available to all members of the food web. Understanding trends in phytoplankton abundance and size structure, how they are influenced by oceanographic conditions, and how they influence fish abundance and size structure are areas of active research.  
\
Status: The mean monthly phytoplankton cell size was `r signif(ann_mean_RptYr,3)` $\mu$m Equivalent Spherical Diameter (ESD) in `r RptYr`. Monthly mean cell size ranged from `r signif(monthly_min_RptYr,3)`–`r signif(monthly_max_RptYr,3)` $\mu$m ESD during the year, within the range of values observed over the period of record (`r signif(monthly_min_PrevYrs,3)`–`r signif(monthly_max_PrevYrs,3)` $\mu$m ESD). Over the period of record, there has been weakly significant decline in monthly median phytoplankton size. Over the time series, median phytoplankton size has declined by `r abs(signif(delta_size_lm,2))` $\mu$m ESD, or by `r abs(signif(delta_size_pct,2))`%. The monthly anomaly has declined as well, by `r abs(signif(delta_size_anom_lm,2))` $\mu$m ESD. Average estimated median phytoplankton size was below average across much of the fishing grounds.    
\
Description: Median phytoplankton cell size can be estimated from satellite remotely sensed SST and ocean color (Barnes et al. 2011). A time series of monthly median phytoplankton cell size averaged over the Hawaiʻi longline region is presented, as well as a time series of anomalies. NOAA CoralTemp (see SST indicator) and ESA CCI data (see ocean color indicator) are used to calculate median phytoplankton cell size.  
\
Timeframe: Monthly  
\
Region: Hawaii longline region: 15° – 45°N, 180° – 120°W    
\
Measurement Platform: Satellite  
\
Data available at: <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/md50_exp>, <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/md50_exp-1998-2009-clim>, and <https://oceanwatch.pifsc.noaa.gov/erddap/griddap/md50_exp-2023-clim>.  
\
Sourced from: Barnes et al. (2011) and NOAA OceanWatch (2024c). Graphics produced in part using Stawitz (2023).  
\
## Additional Information
Median phytoplankton size data are geographically subset and spatially averaged using the pyFerret script `MedPhytoSize_access.jnl` which can be found in the [Median_Phytoplankton_Size folder](https://github.com/pwoodworth-jefcoats/climate-indicators/tree/main/Median_Phytoplankton_Size).  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 negative (~-0.39 min) than positive (~0.35 max). The residuals for the anomaly model were also fairly evenly distributed between -0.1 and 0.1.  
\
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, '/Median_Phyto_map_data_', RptYr, '.csv', sep = ""))
# 
# # Filter to the given year and the anomaly
# maps_RptYr <- maps |> filter(ID == "MD50")
# maps_anom <- maps |> filter(ID == "MD50_anom")
# 
# # Map elements and aesthetics
# seagrass <- nmfs_palette("seagrass")(3)
# pal <- seagrass
# ll_rect_color <- "black" #outline for ll fishing grounds box
# fill_scale <- scale_fill_gradientn(name = "MD50", colors = pal, limits = c(0, 8.25))
# 
# # 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(xlim = c(60, 0), ylim = c(15, 45), 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.1, 0.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(0.1, ceiling(max(maps_anom$layer)), 0.1)),
#                color = "black", lwd = 0.5) # solid positive contours
# 
# p
# 
# pdf(paste(Dir, '/Median_Phyto_Size_map_', RptYr, '.pdf', sep = ""))
# print(p)
# dev.off()
```