bivariate_map <- 
pco2_product_coarse_annual_regression %>%
  filter(year == 2023,
         name %in% c("fgco2", "temperature")) %>%
  select(product, name, lon, lat, resid) %>% 
  pivot_wider(names_from = name,
              values_from = resid) %>% 
  drop_na()

map +
  geom_tile(data = bivariate_map,
            aes(lon, lat, fill = temperature)) +
  scale_fill_continuous_divergingx(palette = "RdBu", rev = TRUE) +
  facet_wrap(~ product, ncol = 2)

map +
  geom_tile(data = bivariate_map,
            aes(lon, lat, fill = fgco2)) +
  scale_fill_continuous_divergingx(palette = "RdBu", rev = TRUE) +
  facet_wrap(~ product, ncol = 2)

dim_set <- 4

bivariate_map <-
  bi_class(
    bivariate_map,
    x = temperature,
    y = fgco2,
    dim = dim_set,
    style = "quantile"
  )

bi_breaks <-
  bi_class_breaks(
    bivariate_map,
    x = temperature,
    y = fgco2,
    dim = dim_set,
    style = "quantile",
    dig_lab = 1
  )

map +
  geom_tile(data = bivariate_map,
            aes(lon, lat, fill = bi_class)) +
  bi_scale_fill(pal = "DkViolet2", dim = dim_set) +
  theme(legend.position = "none") +
  facet_wrap(~ product, ncol = 2)

bi_legend(
  pal = "DkViolet2",
  xlab = labels_breaks("temperature")$i_legend_title,
  ylab = labels_breaks("fgco2")$i_legend_title,
  dim = dim_set,
  pad_width = 2,
  breaks = bi_breaks,
  arrows = FALSE
) +
  theme(
    axis.title.x = element_markdown(),
    axis.title.y = element_markdown(),
    axis.ticks = element_blank(),
    axis.text = element_text(size = 10)
  )

# cowplot::plot_grid(
#   bi_map, bi_legend,
#   rel_heights = c(5, 3),
#   ncol = 1
# )

pco2_product_coarse_annual_regression %>%
  filter(year == 2023) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid)) +
      labs(title =  paste(2023, "anomaly")) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$resid, .01), quantile(.x$resid, .99)),
        oob = squish
      )+
      theme(legend.title = element_markdown()) +
      facet_wrap(~ product, ncol = 2)
  )

pco2_product_coarse_annual_regression_ensemble <-
  pco2_product_coarse_annual_regression %>% 
  filter(year == 2023,
         product %in% pco2_product_list) %>%
  group_by(name, lon, lat) %>%
  summarize(
    resid_sd = sd(resid),
    resid_range = max(resid) - min(resid),
    resid_mean = mean(resid),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n == length(pco2_product_list)) %>% 
  select(-n)

pco2_product_coarse_annual_regression_ensemble %>%
  mutate(product = "Ensemble mean") %>% 
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid_mean)) +
      geom_point(data = .x %>% filter(abs(resid_mean) < resid_sd),
                aes(lon, lat, shape = "Ensemble mean\n< StDev")) +
      labs(title =  paste(2023, "anomaly")) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      )+
      scale_shape_manual(values = 4, name = "") +
      theme(legend.title = element_markdown()) +
      facet_wrap(~ product, ncol = 2)
  )

pco2_product_coarse_annual_regression_ensemble <-
left_join(
    pco2_product_coarse_annual_regression_ensemble,
    pco2_product_coarse_annual_regression %>% 
      filter(year == 2023,
             product %in% pco2_product_list)
  ) %>%
  mutate(`Anomaly offset` = resid - resid_mean) %>% 
  select(name, lon, lat, product, `Anomaly offset`)

pco2_product_coarse_annual_regression_ensemble_baseline <-
  pco2_product_coarse_annual_regression %>% 
  filter(year == 2023,
         product %in% pco2_product_list) %>%
  group_by(name, lon, lat) %>%
  summarize(
    fit_mean = mean(fit),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n == length(pco2_product_list)) %>% 
  select(-n)

pco2_product_coarse_annual_regression_ensemble_baseline <-
left_join(
    pco2_product_coarse_annual_regression_ensemble_baseline,
    pco2_product_coarse_annual_regression %>% 
      filter(year == 2023,
             product %in% pco2_product_list)
  ) %>%
  mutate(`Baseline offset` = fit - fit_mean) %>% 
  select(name, lon, lat, product, `Baseline offset`)

full_join(
  pco2_product_coarse_annual_regression_ensemble,
  pco2_product_coarse_annual_regression_ensemble_baseline
) %>%
  pivot_longer(contains("offset"),
               names_to = "offset") %>% 
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = value)) +
      labs(title =  paste(2023, "offset from ensemble mean")) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$value, .01), quantile(.x$value, .99)),
        oob = squish) +
      theme(legend.title = element_markdown()) +
      facet_grid(product ~ offset)
  )

pco2_product_coarse_monthly_regression %>%
  filter(name %in% name_core) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid)) +
      labs(title = paste(2023, "anomaly")) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$resid, .01), quantile(.x$resid, .99)),
        oob = squish
      )+
      theme(legend.title = element_markdown()) +
      facet_grid(month ~ product) +
      theme(legend.position = "bottom")
  )

pco2_product_coarse_monthly_fCO2_decomposition <-
  fco2_decomposition(pco2_product_coarse_monthly_regression,
                     lon, lat, month)


pco2_product_coarse_monthly_fCO2_decomposition %>%
  group_split(product) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid)) +
      labs(title = paste(2023, "thermal vs non-thermal fCO2 anomaly"),
           subtitle = .x$product) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks("sfco2"),
        limits = c(quantile(.x$resid, .01), quantile(.x$resid, .99)),
        oob = squish
      )+
      theme(legend.title = element_markdown(),
            legend.position = "bottom") +
      facet_grid(
        month ~ name,
        labeller = labeller(name = x_axis_labels)
      )
  )

pco2_product_coarse_monthly_fCO2_decomposition %>%
  filter(product %in% pco2_product_list) %>%
  group_by(name, lon, lat, month) %>%
  summarize(
    resid_sd = sd(resid),
    resid_mean = mean(resid),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n == length(pco2_product_list)) %>% 
  select(-n) %>% 
  mutate(product = "Ensemble mean") %>% 
  group_split(product) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid_mean)) +
      geom_point(data = .x %>% filter(abs(resid_mean) < resid_sd),
                 aes(lon, lat, shape = "Ensemble mean\n< StDev")) +
      labs(title = paste(2023,"thermal vs non-thermal anomaly"),
           subtitle = .x$product) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks("sfco2"),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      )+
      scale_shape_manual(values = 46, name = "") +
      theme(legend.title = element_markdown(),
            legend.position = "bottom") +
      facet_grid(
        month ~ name,
        labeller = labeller(name = x_axis_labels)
      )
  )

pco2_product_coarse_annual_fCO2_decomposition <-
  pco2_product_coarse_monthly_fCO2_decomposition %>% 
  group_by(product, lat, lon, name) %>% 
  summarise(resid = mean(resid, na.rm = TRUE)) %>% 
  ungroup()

map +
  geom_tile(data = pco2_product_coarse_annual_fCO2_decomposition,
            aes(lon, lat, fill = resid)) +
  scale_fill_continuous_divergingx(
    palette = "RdBu",
    rev = TRUE,
    name = labels_breaks("sfco2"),
    limits = c(quantile(pco2_product_coarse_annual_fCO2_decomposition$resid, .01), 
               quantile(pco2_product_coarse_annual_fCO2_decomposition$resid, .99)),
    oob = squish
  )+
  labs(title = paste(2023, "annual mean fCO2 anomaly decomposition")) +
  theme(legend.title = element_markdown(),
        legend.position = "bottom") +
  facet_grid(product ~ name,
             labeller = labeller(name = x_axis_labels)) +
  theme(strip.text.x.top = element_markdown())

pco2_product_coarse_annual_fCO2_decomposition %>%
  filter(product %in% pco2_product_list) %>%
  group_by(name, lon, lat) %>%
  summarize(
    resid_sd = sd(resid),
    resid_mean = mean(resid),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n == length(pco2_product_list)) %>% 
  select(-n) %>% 
  mutate(product = "Ensemble mean") %>% 
  group_split(product) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid_mean)) +
      geom_point(data = .x %>% filter(abs(resid_mean) < resid_sd),
                 aes(lon, lat, shape = "Ensemble mean\n< StDev"))+
      labs(title = paste(2023, "fCO2 anomaly decomposition"),
           subtitle = .x$product) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks("sfco2"),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      ) +
      scale_shape_manual(values = 46, name = "")+
      theme(legend.title = element_markdown(),
            legend.position = "bottom") +
      facet_wrap(
         ~ name,
        labeller = labeller(name = x_axis_labels),
        ncol = 1
      ) +
      theme(strip.text.x.top = element_markdown())
  )

pco2_product_coarse_monthly_flux_attribution <-
  flux_attribution(pco2_product_coarse_monthly_regression,
                   lon, lat, month)

pco2_product_coarse_monthly_flux_attribution %>%
  drop_na() %>% 
  group_split(product) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid)) +
      labs(title = paste(2023, "monthly flux anomaly decomposition"),
           subtitle = .x$product) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks("fgco2"),
        limits = c(quantile(.x$resid, .01), quantile(.x$resid, .99)),
        oob = squish
      ) +
      theme(legend.title = element_markdown(), 
            legend.position = "bottom") +
      facet_grid(
        month ~ name,
        labeller = labeller(name = x_axis_labels)
      ) +
      theme(
        strip.text.x.top = element_markdown()
      )
  )

pco2_product_coarse_monthly_flux_attribution %>%
  drop_na() %>% 
  filter(product %in% pco2_product_list) %>%
  group_by(name, lon, lat, month) %>%
  summarize(
    resid_sd = sd(resid),
    resid_mean = mean(resid),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n == length(pco2_product_list)) %>% 
  select(-n) %>% 
  mutate(product = "Ensemble mean") %>% 
  drop_na() %>% 
  group_split(product) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid_mean)) +
      geom_point(data = .x %>% filter(abs(resid_mean) < resid_sd),
                 aes(lon, lat, shape = "Ensemble mean\n< StDev"))+
      labs(title = paste(2023, "flux anomaly decomposition"),
           subtitle = .x$product) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks("fgco2"),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      )+
      scale_shape_manual(values = 46, name = "") +
      theme(legend.title = element_markdown(),
            legend.position = "bottom") +
      facet_grid(
        month ~ name,
        labeller = labeller(name = x_axis_labels)
      ) +
      theme(strip.text.x.top = element_markdown())
  )

pco2_product_coarse_annual_flux_attribution <-
  pco2_product_coarse_monthly_flux_attribution %>% 
  group_by(product, lat, lon, name) %>% 
  summarise(resid = mean(resid, na.rm = TRUE)) %>% 
  ungroup()

map +
  geom_tile(data = pco2_product_coarse_annual_flux_attribution,
            aes(lon, lat, fill = resid)) +
  scale_fill_continuous_divergingx(
    palette = "RdBu",
    rev = TRUE,
    name = labels_breaks("fgco2"),
    limits = c(quantile(pco2_product_coarse_annual_flux_attribution$resid, .01, na.rm = TRUE), 
               quantile(pco2_product_coarse_annual_flux_attribution$resid, .99, na.rm = TRUE)),
    oob = squish
  )+
  labs(title = paste(2023, "annual mean flux anomaly decomposition")) +
  theme(legend.title = element_markdown(),
        legend.position = "bottom") +
  facet_grid(product ~ name,
             labeller = labeller(name = x_axis_labels)) +
  theme(strip.text.x.top = element_markdown())

pco2_product_coarse_annual_flux_attribution %>%
  filter(product %in% pco2_product_list) %>%
  group_by(name, lon, lat) %>%
  summarize(
    resid_sd = sd(resid),
    resid_mean = mean(resid),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n == length(pco2_product_list)) %>% 
  select(-n) %>% 
  mutate(product = "Ensemble mean") %>% 
  drop_na() %>% 
  group_split(product) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = resid_mean)) +
      geom_point(data = .x %>% filter(abs(resid_mean) < resid_sd),
                 aes(lon, lat, shape = "Ensemble mean\n< StDev"))+
      labs(title = paste(2023, "flux anomaly decomposition"),
           subtitle = .x$product) +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks("fgco2"),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      )+
      scale_shape_manual(values = 46, name = "") +
      theme(legend.title = element_markdown(),
            legend.position = "bottom") +
      facet_wrap(
         ~ name,
        labeller = labeller(name = x_axis_labels),
        ncol = 1
      ) +
      theme(strip.text.x.top = element_markdown())
  )

pco2_product_hovmoeller_monthly_regression %>%
  filter(name %in% name_core) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(decimal, lat, fill = resid)) +
      geom_raster() +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$resid,.01),quantile(.x$resid,.99)),
        oob = squish
      ) +
      theme(legend.title = element_markdown()) +
      coord_cartesian(expand = 0) +
      labs(title = "Monthly mean anomalies",
           y = "Latitude") +
      theme(axis.title.x = element_blank()) +
      facet_wrap(~ product, ncol = 1)
  )

pco2_product_hovmoeller_monthly_regression_ensemble <-
  pco2_product_hovmoeller_monthly_regression %>% 
  group_by(name, decimal, lat) %>%
  summarize(
    resid_range = max(resid) - min(resid),
    resid_mean = mean(resid),
    n = n()
  ) %>%
  ungroup() %>%
  filter(n > 1)
  

pco2_product_hovmoeller_monthly_regression_ensemble %>%
  mutate(product = "Ensemble mean") %>% 
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(decimal, lat, fill = resid_mean)) +
      geom_raster() +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      ) +
      theme(legend.title = element_markdown()) +
      coord_cartesian(expand = 0) +
      labs(title = "Monthly mean anomalies",
           y = "Latitude") +
      theme(axis.title.x = element_blank()) +
      facet_wrap( ~ product, ncol = 1)
  )

left_join(
    pco2_product_hovmoeller_monthly_regression_ensemble,
    pco2_product_hovmoeller_monthly_regression
  ) %>%
  mutate(resid_offset = resid - resid_mean) %>% 
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(decimal, lat, fill = resid_offset)) +
      geom_raster() +
      scale_fill_continuous_divergingx(
        palette = "RdBu",
        rev = TRUE,
        name = labels_breaks(.x %>% distinct(name)),
        limits = c(quantile(.x$resid_mean, .01), quantile(.x$resid_mean, .99)),
        oob = squish
      ) +
      theme(legend.title = element_markdown()) +
      coord_cartesian(expand = 0)+
      labs(title = "Monthly offset from ensemble mean anomalies",
           y = "Latitude") +
      theme(axis.title.x = element_blank()) +
      facet_wrap( ~ product, ncol = 1)
  )

pco2_product_annual_detrended %>% 
  filter(biome %in% "Global",
         name %in% c("fgco2", "temperature"),
         product %in% "OceanSODA") %>%
  p_season(title = "Anomalies from predicted annual mean | Global")

pco2_product_annual_detrended %>% 
  filter(biome %in% "Global",
         name %in% c("fgco2", "temperature"),
         product %in% pco2_product_list) %>%
  p_season(title = "Anomalies from predicted annual mean | Global")

pco2_product_annual_detrended %>% 
  filter(biome %in% "Global",
         name %in% name_core) %>%
  p_season(title = "Anomalies from predicted annual mean | Global")

pco2_product_annual_detrended %>%
  filter(biome %in% super_biomes,
         name %in% name_core) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(~ p_season(
    df = .x,
    title = paste("Anomalies from predicted annual mean |", .x$biome)
  ))

pco2_product_annual_detrended %>%
  filter(biome %in% key_biomes,
         name %in% name_core) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(~ p_season(
    df = .x,
    title = paste("Anomalies from predicted annual mean |", .x$biome)
  ))

pco2_product_monthly_detrended %>% 
  filter(biome %in% "Global",
         name %in% name_core) %>%
  p_season(title = "Anomalies from predicted monthly mean | Global")

pco2_product_monthly_detrended %>%
  filter(biome %in% super_biomes,
         name %in% name_core) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ p_season(df = .x,
               title = paste("Anomalies from predicted monthly mean |", .x$biome))
  )

pco2_product_monthly_detrended %>%
  filter(biome %in% key_biomes,
         name %in% name_core) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(~ p_season(
    df = .x,
    title = paste("Anomalies from predicted monthly mean |", .x$biome)
  ))

pco2_product_annual_regression %>%
  filter(biome %in% "Global",
         name %in% name_core) %>%
  select(-c(value, fit)) %>%
  pivot_wider(values_from = resid) %>%
  pivot_longer(-c(product, year, biome, fgco2_int))  %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(value, fgco2_int)) +
      geom_hline(yintercept = 0) +
      geom_point(
        data = . %>% filter(year != 2023),
        aes(fill = year),
        shape = 21
      ) +
      geom_smooth(
        data = . %>% filter(year != 2023),
        method = "lm",
        se = FALSE,
        fullrange = TRUE,
        aes(col = paste("Regression fit\nexcl.", 2023))
      ) +
      scale_color_grey() +
      scale_fill_grayC() +
      new_scale_fill() +
      geom_point(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(fill = as.factor(year)),
        shape = 21,
        size = 2
      )  +
      scale_fill_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      labs(y = labels_breaks("fgco2_int")$i_legend_title,
           x = labels_breaks(unique(.x$name))$i_legend_title) +
      facet_grid(
        biome ~ product,
        scales = "free_y"
      ) +
      theme(
        axis.title.x = element_markdown(),
        axis.title.y = element_markdown(),
        legend.title = element_blank()
      )
  )

pco2_product_annual_regression %>%
  filter(biome %in% super_biomes,
         name %in% name_core) %>%
  select(-c(value, fit)) %>%
  pivot_wider(values_from = resid) %>%
  pivot_longer(-c(product, year, biome, fgco2_int))  %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(value, fgco2_int)) +
      geom_hline(yintercept = 0) +
      geom_point(
        data = . %>% filter(year != 2023),
        aes(fill = year),
        shape = 21
      ) +
      geom_smooth(
        data = . %>% filter(year != 2023),
        method = "lm",
        se = FALSE,
        fullrange = TRUE,
        aes(col = paste("Regression fit\nexcl.", 2023))
      ) +
      scale_color_grey() +
      scale_fill_grayC() +
      new_scale_fill() +
      geom_point(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(fill = as.factor(year)),
        shape = 21,
        size = 2
      )  +
      scale_fill_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      labs(y = labels_breaks("fgco2_int")$i_legend_title,
           x = labels_breaks(unique(.x$name))$i_legend_title) +
      facet_grid(
        biome ~ product,
        scales = "free_y"
      ) +
      theme(
        axis.title.x = element_markdown(),
        axis.title.y = element_markdown(),
        legend.title = element_blank()
      )
  )