Code 
library(tidyverse)
library(here)
library(patchwork)

# epc data by subsegment
load(file = here('data-clean/epcwq_clean.RData'))

epc <- epcwq3 |>
  filter(param == 'Chla') |>
  mutate(
    date = ymd(date),
    yr = year(date),
    mo = month(date)
  ) |> 
  select(subsegment = subseg, mo, yr, obs = value, date) |> 
  summarise(
    obs = mean(obs, na.rm = TRUE),
    .by = c('subsegment', 'mo', 'yr', 'date')
  )
  
# casm data by subsegment
casm <- list(
    "prd486" = read.csv(here("data-raw/casm486.csv")),
    "prdhist" = read.csv(here("data-raw/casmhistproj.csv"))
  ) |> 
  enframe('run', 'data') |> 
  mutate(
    data = map(data, function(x) x |> 
                 pivot_longer(cols = -c(1:2), names_to = 'yr', values_to = 'prd') |> 
                 rename(mo = Month) |> 
                 mutate(
                   yr = gsub('X', '', yr) |> as.numeric(),
                   mo = factor(mo, levels = c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'), 
                               labels = c(1:12)) |> as.numeric(),
                   date = as.Date(paste0(yr, '-', mo, '-01'))
                 )
    )
  ) |> 
  unnest('data') |> 
  pivot_wider(names_from = 'run', values_from = 'prd')

##
# combine epc and casm data

# monthly
mo <- casm |> 
  filter(yr <= 2021) |> 
  inner_join(epc, by = c('yr', 'mo', 'date', 'subsegment')) |> 
  mutate(
    subsegment = factor(subsegment, levels = c('NW', 'NE', 'CW', 'CE', 'SW', 'SE'))
  )

# annual
ann <- mo |> 
  summarise(
    obs = mean(obs, na.rm = TRUE),
    prd486 = mean(prd486, na.rm = TRUE),
    prdhist = mean(prdhist, na.rm = TRUE),
    .by = c('subsegment', 'yr')
  )

Report plot

Code 
toplo <- casm |> 
  summarise(
    prd486 = mean(prd486, na.rm = T),
    prdhist = mean(prdhist, na.rm = T),
    .by = yr
  )
ggplot(toplo, aes(x = yr)) + 
  geom_line(aes(y = prdhist, color = 'Predicted historical and projected')) +
  geom_line(aes(y = prd486, color = 'Predicted 486 tons N/yr or less')) + 
  theme_minimal() +
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  labs(
    x = NULL,
    y = 'Chlorophyll-a (µg/L)',
    colour = NULL
  )

Annual all

Code 
annall <- ann |> 
  summarise(
    obs = mean(obs, na.rm = TRUE),
    prd486 = mean(prd486, na.rm = TRUE),
    prdhist = mean(prdhist, na.rm = TRUE), 
    .by = c('yr')
  )

p1 <- ggplot(annall, aes(x = yr)) + 
  geom_line(aes(y = prdhist, color = 'Predicted')) +
  geom_line(aes(y = obs, color = 'Observed')) + 
  theme_minimal() +
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  labs(
    x = NULL,
    y = 'Chlorophyll-a (µg/L)',
    colour = NULL
  )
p2 <- ggplot(annall, aes(x = obs, y = prdhist)) + 
  geom_point() + 
  geom_abline(aes(linetype = '1:1', slope = 1, intercept = 0)) +
  geom_smooth(se = F, method = 'lm', formula = y ~ x, aes(linetype = 'reg')) +
  theme_minimal() +
  scale_linetype_manual(values = c('1:1' = 'dashed', 'reg' = 'solid')) +
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  # coord_equal() + 
  labs(
    x = 'Observed',
    y = 'Predicted',
    linetype = NULL
  )

p1 + p2 + plot_layout(ncol = 2, widths = c(1, 0.5))

Code 
mod <- lm(prdhist ~ obs, annall)
r2 <- round(summary(mod)$r.squared, 2)

R2 = 0.48

Monthly all

Code 
moall <- mo |> 
  summarise(
    obs = mean(obs, na.rm = TRUE),
    prd486 = mean(prd486, na.rm = TRUE),
    prdhist = mean(prdhist, na.rm = TRUE), 
    .by = c('yr', 'mo', 'date')
  )

p1 <- ggplot(moall, aes(x = date)) + 
  geom_line(aes(y = prdhist, color = 'Predicted')) +
  geom_line(aes(y = obs, color = 'Observed')) + 
  theme_minimal() +
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  labs(
    x = NULL,
    y = 'Chlorophyll-a (µg/L)',
    colour = NULL
  )
p2 <- ggplot(moall, aes(x = obs, y = prdhist)) + 
  geom_point() + 
  geom_abline(aes(linetype = '1:1', slope = 1, intercept = 0)) +
  geom_smooth(se = F, method = 'lm', formula = y ~ x, aes(linetype = 'reg')) +
  theme_minimal() +
  scale_linetype_manual(values = c('1:1' = 'dashed', 'reg' = 'solid')) +
  # coord_equal() + 
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  labs(
    x = 'Observed',
    y = 'Predicted',
    linetype = NULL
  )

p1 + p2 + plot_layout(ncol = 2, widths = c(1, 0.5))

Code 
mod <- lm(prdhist ~ obs, moall)
r2 <- round(summary(mod)$r.squared, 2)

R2 = 0.16

Annual by subsegment

Code 
p1 <- ggplot(ann, aes(x = yr)) + 
  geom_line(aes(y = prdhist, color = 'Predicted')) +
  geom_line(aes(y = obs, color = 'Observed')) + 
  theme_minimal() +
  facet_wrap(~subsegment, ncol = 1) + 
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  labs(
    x = NULL,
    y = 'Chlorophyll-a (µg/L)',
    colour = NULL
  )
p2 <- ggplot(ann, aes(x = obs, y = prdhist)) + 
  geom_point() + 
  geom_abline(aes(linetype = '1:1', slope = 1, intercept = 0)) +
  geom_smooth(se = F, method = 'lm', formula = y ~ x, aes(linetype = 'reg')) +
  facet_wrap(~subsegment, ncol = 1) + 
  theme_minimal() +
  scale_linetype_manual(values = c('1:1' = 'dashed', 'reg' = 'solid')) +
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  # coord_equal() + 
  labs(
    x = 'Observed',
    y = 'Predicted',
    linetype = NULL
  )

p1 + p2 + plot_layout(ncol = 2, widths = c(1, 0.5))

Code 
mods <- ann |> 
  summarise(
    R2 = round(summary(lm(prdhist ~ obs))$r.squared, 2),
    .by = subsegment
  )
knitr::kable(mods)
subsegment R2
NW 0.02
NE 0.10
CW 0.74
CE 0.16
SW 0.63
SE 0.17

Monthly by subsegment

Code 
p1 <- ggplot(mo, aes(x = date)) + 
  geom_line(aes(y = prdhist, color = 'Predicted')) +
  geom_line(aes(y = obs, color = 'Observed')) + 
  theme_minimal() +
  facet_wrap(~subsegment, ncol = 1) + 
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  labs(
    x = NULL,
    y = 'Chlorophyll-a (µg/L)',
    colour = NULL
  )
p2 <- ggplot(mo, aes(x = obs, y = prdhist)) + 
  geom_point() + 
  geom_abline(aes(linetype = '1:1', slope = 1, intercept = 0)) +
  geom_smooth(se = F, method = 'lm', formula = y ~ x, aes(linetype = 'reg')) +
  facet_wrap(~subsegment, ncol = 1) + 
  theme_minimal() +
  scale_linetype_manual(values = c('1:1' = 'dashed', 'reg' = 'solid')) +
  theme(
    legend.position = 'top',
    panel.grid.minor = element_blank()
  ) +
  # coord_equal() + 
  labs(
    x = 'Observed',
    y = 'Predicted',
    linetype = NULL
  )

p1 + p2 + plot_layout(ncol = 2, widths = c(1, 0.5))

Code 
mods <- mo |> 
  summarise(
    R2 = round(summary(lm(prdhist ~ obs))$r.squared, 2),
    .by = subsegment
  )
knitr::kable(mods)
subsegment R2
NW 0.17
NE 0.03
CW 0.08
CE 0.15
SW 0.04
SE 0.03