pco2_product_hovmoeller_monthly_annual <-
  pco2_product %>%
  mutate(fgco2_int = fgco2) %>% 
  select(-c(lon, time, month, biome)) %>%
  group_by(year, lat) %>%
  summarise(across(-c(fgco2_int, area),
                   ~ weighted.mean(., area, na.rm = TRUE)),
            across(fgco2_int,
                   ~ sum(. * area, na.rm = TRUE) * 12.01 * 1e-15)) %>%
  ungroup() %>%
  rename(fgco2_hov = fgco2_int) %>% 
  filter(fgco2_hov != 0)

pco2_product_hovmoeller_monthly_annual <-
  pco2_product_hovmoeller_monthly_annual %>%
  pivot_longer(-c(year, lat)) %>% 
  drop_na()

pco2_product_hovmoeller_monthly_annual %>%
  filter(!(name %in% name_divergent)) %>% 
  group_split(name) %>%
  # tail(5) %>%
  map(
    ~ ggplot(data = .x,
             aes(year, lat, fill = value)) +
      geom_raster() +
      scale_fill_viridis_c(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      coord_cartesian(expand = 0) +
      labs(title = "Annual means",
           y = "Latitude") +
      theme(axis.title.x = element_blank())
  )

pco2_product_hovmoeller_monthly_annual %>%
  filter(name %in% name_divergent) %>% 
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(year, lat, fill = value)) +
      geom_raster() +
      scale_fill_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      coord_cartesian(expand = 0) +
      labs(title = "Annual means",
           y = "Latitude") +
      theme(axis.title.x = element_blank())
  )

pco2_product_hovmoeller_monthly <-
  pco2_product %>%
  mutate(fgco2_int = fgco2) %>% 
  select(-c(lon, time, biome)) %>%
  group_by(year, month, lat) %>%
  summarise(across(-c(fgco2_int, area),
                   ~ weighted.mean(., area, na.rm = TRUE)),
            across(fgco2_int,
                   ~ sum(. * area, na.rm = TRUE) * 12.01 * 1e-15)) %>%
  ungroup() %>%
  rename(fgco2_hov = fgco2_int) %>% 
  filter(fgco2_hov != 0)


pco2_product_hovmoeller_monthly <-
  pco2_product_hovmoeller_monthly %>%
  pivot_longer(-c(year, month, lat)) %>% 
  drop_na()

pco2_product_hovmoeller_monthly <-
  pco2_product_hovmoeller_monthly %>% 
  mutate(decimal = year + (month-1) / 12)

pco2_product_hovmoeller_monthly %>%
  filter(!(name %in% name_divergent)) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(decimal, lat, fill = value)) +
      geom_raster() +
      scale_fill_viridis_c(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      labs(title = "Monthly means",
           y = "Latitude") +
      coord_cartesian(expand = 0) +
      theme(axis.title.x = element_blank())
  )

pco2_product_hovmoeller_monthly %>%
  filter(name %in% name_divergent) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(decimal, lat, fill = value)) +
      geom_raster() +
      scale_fill_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      labs(title = "Monthly means",
           y = "Latitude") +
      coord_cartesian(expand = 0) +
      theme(axis.title.x = element_blank())
  )

pCO2productanalysis_2023 <-
  knitr::knit_expand(
    file = here::here("analysis/child/pCO2_product_analysis.Rmd"),
    product_name = "CMEMS",
    year_anom = 2023
  )

pco2_product_coarse_annual_regression <-
  pco2_product_coarse_annual %>%
  drop_na() %>% 
  anomaly_determination(lon, lat)

pco2_product_coarse_annual_regression <-
  pco2_product_coarse_annual_regression %>%
  drop_na()

pco2_product_coarse_annual_regression %>%
  filter(year == 2023,
         !(name %in% name_divergent)) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = value)) +
      labs(title = paste("Annual mean", 2023)) +
      scale_fill_viridis_c(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown())
  )

pco2_product_coarse_annual_regression %>%
  filter(year == 2023,
         name %in% name_divergent) %>%
  group_split(name) %>% 
  # head(1) %>%
  map( ~ map +
         geom_tile(data = .x,
                   aes(lon, lat, fill = value)) +
         labs(title = paste("Annual mean", 2023)) +
         scale_fill_divergent(
           name = labels_breaks(.x %>% distinct(name))) +
         theme(legend.title = element_markdown())
  )

pco2_product_coarse_annual_regression <-
  pco2_product_coarse_annual_regression %>%
  group_by(name) %>%
  filter(year %in% c(min(year), max(year), 2023)) %>%
  ungroup()

pco2_product_coarse_annual_regression %>%
  select(-c(value, resid)) %>%
  filter(year %in% c(min(year), max(year))) %>%
  arrange(year) %>%
  group_by(lon, lat, name) %>%
  mutate(change = fit - lag(fit),
         period = paste(lag(year), year, sep = "-")) %>%
  ungroup() %>%
  filter(!is.na(change)) %>% 
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = change)) +
      labs(title =  paste("Change: ",.x$period)) +
      scale_fill_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown())
  )

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_divergent(
           name = labels_breaks(.x %>% distinct(name))) +
         theme(legend.title = element_markdown())
  )

pco2_product_coarse_annual_regression %>%
  write_csv(paste0("../data/","CMEMS","_","2023","_anomaly_map_annual.csv"))

pco2_product_coarse_monthly_regression <-
  pco2_product_coarse_monthly %>%
  drop_na() %>% 
  anomaly_determination(lon, lat, month)

pco2_product_coarse_monthly_regression <-
  pco2_product_coarse_monthly_regression %>% 
  drop_na()


pco2_product_coarse_monthly_regression %>%
  filter(year == 2023,
         !(name %in% name_divergent)) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = value)) +
      labs(title = paste("Monthly means", 2023)) +
      scale_fill_viridis_c(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      facet_wrap( ~ month, ncol = 2)
  )

pco2_product_coarse_monthly_regression %>%
  filter(year == 2023,
         name %in% name_divergent) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = value)) +
      labs(title = paste("Monthly means", 2023)) +
      scale_fill_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      facet_wrap( ~ month, ncol = 2)
  )

pco2_product_coarse_monthly_regression %>%
  group_by(name) %>%
  filter(year %in% c(min(year), max(year))) %>%
  ungroup() %>%
  select(-c(value, resid)) %>%
  arrange(year) %>%
  group_by(lon, lat, name, month) %>%
  mutate(change = fit - lag(fit),
         period = paste(lag(year), year, sep = "-")) %>%
  ungroup() %>%
  filter(!is.na(change)) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ map +
      geom_tile(data = .x,
                aes(lon, lat, fill = change)) +
      labs(title =  paste("Change: ", .x$period)) +
      scale_fill_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      facet_wrap(~ month, ncol = 2)
  )

pco2_product_coarse_monthly_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_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      facet_wrap( ~ month, ncol = 2)
  )

pco2_product_coarse_monthly_regression %>%
  filter(year == 2023) %>%
  write_csv(paste0("../data/","CMEMS","_","2023","_anomaly_map_monthly.csv"))

pco2_product_hovmoeller_monthly_regression %>%
  filter(between(year, 2023-2, 2023)) %>%
  group_split(name) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(decimal, lat, fill = resid)) +
      geom_raster() +
      scale_fill_divergent(name = labels_breaks(.x %>% distinct(name))) +
      theme(legend.title = element_markdown()) +
      coord_cartesian(expand = 0) +
      labs(title = "Monthly mean anomalies",
           y = "Latitude") +
      theme(axis.title.x = element_blank())
  )

fig.height <- pco2_product_monthly %>% 
  distinct(name) %>% 
  nrow()

fig.height <- (fig.height + 2) * 0.1

pco2_product_monthly %>%
  filter(biome %in% "Global") %>%
  ggplot(aes(month, value, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Absolute values | Global") +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    legend.title = element_blank(),
    axis.title.y = element_blank()
  )

pco2_product_monthly %>%
  filter(biome %in% key_biomes) %>%
  ggplot(aes(month, value, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Absolute values | Selected biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    legend.title = element_blank(),
    axis.title.y = element_blank()
  )

pco2_product_monthly %>%
  filter(biome %in% key_biomes) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(month, value, group = as.factor(year))) +
      geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
                aes(col = year)) +
      scale_color_grayC() +
      new_scale_color() +
      geom_path(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(col = as.factor(year)),
        linewidth = 1
      ) +
      scale_color_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = paste("Absolute values |", .x$biome)) +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    legend.title = element_blank(),
    axis.title.y = element_blank()
  )
  )

pco2_product_monthly %>%
  filter(biome %in% super_biomes) %>%
  ggplot(aes(month, value, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Absolute values | Selected super biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    legend.title = element_blank(),
    axis.title.y = element_blank()
  )

pco2_product_monthly %>%
  filter(biome %in% super_biomes) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(month, value, group = as.factor(year))) +
      geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
                aes(col = year)) +
      scale_color_grayC() +
      new_scale_color() +
      geom_path(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(col = as.factor(year)),
        linewidth = 1
      ) +
      scale_color_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = paste("Absolute values |", .x$biome)) +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    legend.title = element_blank(),
    axis.title.y = element_blank()
  )
  )

pco2_product_annual <-
  pco2_product_monthly %>%
  group_by(year, biome, name) %>%
  summarise(value = mean(value)) %>%
  ungroup()
  
pco2_product_annual %>% 
  filter(biome %in% "Global") %>%
  ggplot(aes(year, value)) +
  geom_path() +
  geom_point(data = . %>% filter(year != 2023)) +
  geom_point(data = . %>% filter(year == 2023),
             shape = 1) +
  geom_smooth(data = . %>% filter(year != 2023,
                                  !(name %in% name_quadratic_fit)),
              method = "lm",
              fullrange = TRUE,
              aes(col = "linear"),
              se = FALSE) +
  geom_smooth(data = . %>% filter(year != 2023,
                                  name %in% name_quadratic_fit),
              method = "lm",
              fullrange = TRUE,
              formula = y ~ x + I(x^2),
              aes(col = "quadratic"),
              se = FALSE) +
  scale_color_brewer(
    palette = "Set1",
    name = paste("Regression fit\nexcl.", 2023)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Annual mean trends | Global") +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank()
  )

pco2_product_annual %>% 
  filter(biome %in% key_biomes) %>%
  ggplot(aes(year, value)) +
  geom_path() +
  geom_point(data = . %>% filter(year != 2023),
             size = 0.2) +
  geom_point(data = . %>% filter(year == 2023),
             shape = 1) +
  geom_smooth(data = . %>% filter(year != 2023,
                                  !(name %in% name_quadratic_fit)),
              method = "lm",
              fullrange = TRUE,
              aes(col = "linear"),
              se = FALSE) +
  geom_smooth(data = . %>% filter(year != 2023,
                                  name %in% name_quadratic_fit),
              method = "lm",
              fullrange = TRUE,
              formula = y ~ x + I(x^2),
              aes(col = "quadratic"),
              se = FALSE) +
  scale_color_brewer(
    palette = "Set1",
    name = paste("Regression fit\nexcl.", 2023)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Annual mean trends | Selected biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank()
  )

pco2_product_annual %>% 
  filter(biome %in% super_biomes) %>%
  ggplot(aes(year, value)) +
  geom_path() +
  geom_point(data = . %>% filter(year != 2023),
             size = 0.2) +
  geom_point(data = . %>% filter(year == 2023),
             shape = 1) +
  geom_smooth(data = . %>% filter(year != 2023,
                                  !(name %in% name_quadratic_fit)),
              method = "lm",
              fullrange = TRUE,
              aes(col = "linear"),
              se = FALSE) +
  geom_smooth(data = . %>% filter(year != 2023,
                                  name %in% name_quadratic_fit),
              method = "lm",
              fullrange = TRUE,
              formula = y ~ x + I(x^2),
              aes(col = "quadratic"),
              se = FALSE) +
  scale_color_brewer(
    palette = "Set1",
    name = paste("Regression fit\nexcl.", 2023)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Annual mean trends | Selected super biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank()
  )

pco2_product_annual_regression <-
  pco2_product_annual %>%
  anomaly_determination(biome)
  

pco2_product_annual_regression %>%
  filter(biome %in% "Global") %>%
  ggplot(aes(year, resid)) +
  geom_hline(yintercept = 0) +
  geom_path() +
  geom_point(data = . %>% filter(!between(year, 2023-1, 2023)),
             aes(fill = year),
             shape = 21) +
  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)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Residual from fit to annual means (actual - predicted) | Global") +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank(),
    legend.title = element_blank()
  )

pco2_product_annual_regression %>%
  filter(biome %in% key_biomes) %>%
  ggplot(aes(year, resid)) +
  geom_hline(yintercept = 0) +
  geom_path() +
  geom_point(data = . %>% filter(!between(year, 2023-1, 2023)),
             aes(fill = year),
             shape = 21) +
  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)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Residual from fit to annual means (actual - predicted) | Selected biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank(),
    legend.title = element_blank()
  )

pco2_product_annual_regression %>%
  filter(biome %in% super_biomes) %>%
  ggplot(aes(year, resid)) +
  geom_hline(yintercept = 0) +
  geom_path() +
  geom_point(data = . %>% filter(!between(year, 2023-1, 2023)),
             aes(fill = year),
             shape = 21) +
  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)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Residual from fit to annual means (actual - predicted) | Selected super biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank(),
    legend.title = element_blank()
  )

pco2_product_annual_regression %>%
  filter(biome != "global") %>%
  ggplot(aes(biome, resid)) +
  geom_hline(yintercept = 0) +
  geom_jitter(data = . %>% filter(!between(year, 2023-1, 2023)),
             aes(fill = year),
             shape = 21, width = 0.2) +
  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(title = "Residual from fit to annual means (actual - predicted) | All biomes") +
  facet_grid(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank(),
    legend.title = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
  )

pco2_product_annual_regression %>%
  write_csv(paste0("../data/","CMEMS","_","2023","_biome_annual_regression.csv"))

pco2_product_annual_detrended <-
  full_join(pco2_product_monthly,
            pco2_product_annual_regression %>% select(-c(value, resid))) %>%
  mutate(resid = value - fit)

pco2_product_annual_detrended %>% 
  filter(biome %in% "Global") %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted annual mean | Global") +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank()
  )

pco2_product_annual_detrended %>% 
  filter(biome %in% key_biomes) %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted annual mean | Selected biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
  )

pco2_product_annual_detrended %>% 
  filter(biome %in% super_biomes) %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted annual mean | Selected super biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
  )

pco2_product_annual_detrended %>%
  filter(biome %in% key_biomes) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(month, resid, group = as.factor(year))) +
      geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
                aes(col = year)) +
      scale_color_grayC() +
      new_scale_color() +
      geom_path(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(col = as.factor(year)),
        linewidth = 1
      ) +
      scale_color_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
      labs(title = paste("Anomalies from predicted annual mean |", .x$biome)) +
      facet_wrap(
        name ~ .,
        scales = "free_y",
        labeller = labeller(name = x_axis_labels),
        strip.position = "left",
        ncol = 2
      ) +
      theme(
        strip.text.y.left = element_markdown(),
        strip.placement = "outside",
        strip.background.y = element_blank(),
        axis.title = element_blank(),
        legend.title = element_blank()
      )
  )

pco2_product_annual_detrended %>%
  filter(biome %in% super_biomes) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(month, resid, group = as.factor(year))) +
      geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
                aes(col = year)) +
      scale_color_grayC() +
      new_scale_color() +
      geom_path(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(col = as.factor(year)),
        linewidth = 1
      ) +
      scale_color_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
      labs(title = paste("Anomalies from predicted annual mean |", .x$biome)) +
      facet_wrap(
        name ~ .,
        scales = "free_y",
        labeller = labeller(name = x_axis_labels),
        strip.position = "left",
        ncol = 2
      ) +
      theme(
        strip.text.y.left = element_markdown(),
        strip.placement = "outside",
        strip.background.y = element_blank(),
        axis.title = element_blank(),
        legend.title = element_blank()
      )
  )

pco2_product_annual_detrended %>%
  write_csv(paste0("../data/","CMEMS","_","2023","_biome_annual_detrended.csv"))

pco2_product_monthly %>% 
  filter(biome %in% "Global") %>%
  mutate(month = as.factor(month)) %>% 
  ggplot(aes(year, value, col = month)) +
  # geom_point() +
  geom_smooth(data = . %>% filter(year != 2023,
                                  !(name %in% name_quadratic_fit)),
              method = "lm",
              se = FALSE,
              fullrange = TRUE) +
  geom_smooth(
    data = . %>% filter(year != 2023,
                        name %in% name_quadratic_fit),
    method = "lm",
    fullrange = TRUE,
    formula = y ~ x + I(x ^ 2),
    se = FALSE
  ) +
  scale_color_scico_d(palette = "romaO", 
    name = paste("Regression fit\nexcl.", 2023)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Monthly mean trends | Global") +
  facet_wrap(name ~ .,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             strip.position = "left",
             ncol = 2) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank()
  )

pco2_product_monthly %>% 
  filter(biome %in% key_biomes) %>%
  mutate(month = as.factor(month)) %>% 
  ggplot(aes(year, value, col = month)) +
  # geom_point() +
  geom_smooth(data = . %>% filter(year != 2023,
                                  !(name %in% name_quadratic_fit)),
              method = "lm",
              fullrange = TRUE,
              se = FALSE) +
  geom_smooth(
    data = . %>% filter(year != 2023,
                        name %in% name_quadratic_fit),
    method = "lm",
    fullrange = TRUE,
    formula = y ~ x + I(x ^ 2),
    se = FALSE
  ) +
  scale_color_scico_d(palette = "romaO",
    name = paste("Regression fit\nexcl.", 2023)) +
  scale_x_continuous(breaks = seq(1980, 2020, 20)) +
  labs(title = "Monthly mean trends | Selected biomes") +
  facet_grid(name ~ biome,
             scales = "free_y",
             labeller = labeller(name = x_axis_labels),
             switch = "y") +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title = element_blank()
  )

pco2_product_monthly_regression <-
  pco2_product_monthly %>%
  anomaly_determination(biome, month)


pco2_product_monthly_regression %>%
  filter(biome %in% "Global") %>%
  group_split(month) %>%
  head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(year, resid)) +
      geom_hline(yintercept = 0) +
      geom_path() +
      geom_point(
        data = . %>% filter(!between(year, 2023-1, 2023)),
        aes(fill = year),
        shape = 21
      ) +
      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)
      ) +
      scale_x_continuous(breaks = seq(1980, 2020, 20)) +
      labs(title = "Residual from fit to monthly means (actual - predicted) | Global",
           subtitle = paste("Month:", .x$month)) +
      facet_wrap(
        name ~ .,
        scales = "free_y",
        labeller = labeller(name = x_axis_labels),
        strip.position = "left",
        ncol = 2
      ) +
      theme(
        strip.text.y.left = element_markdown(),
        strip.placement = "outside",
        strip.background.y = element_blank(),
        axis.title = element_blank(),
        legend.title = element_blank()
      )
  )

pco2_product_monthly_regression %>%
  filter(biome %in% key_biomes) %>%
  group_split(month) %>%
  head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(year, resid)) +
      geom_hline(yintercept = 0) +
      geom_path() +
      geom_point(
        data = . %>% filter(!between(year, 2023-1, 2023)),
        aes(fill = year),
        shape = 21
      ) +
      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)
      ) +
      scale_x_continuous(breaks = seq(1980, 2020, 20)) +
      labs(title = "Residual from fit to monthly means (actual - predicted) | Selected biomes",
           subtitle = paste("Month:", .x$month)) +
      facet_grid(
        name ~ biome,
        scales = "free_y",
        labeller = labeller(name = x_axis_labels),
        switch = "y"
      ) +
      theme(
        strip.text.y.left = element_markdown(),
        strip.placement = "outside",
        strip.background.y = element_blank(),
        axis.title = element_blank(),
        legend.title = element_blank()
      )
  )

pco2_product_monthly_regression %>%
  write_csv(paste0("../data/","CMEMS","_","2023","_biome_monthly_regression.csv"))

pco2_product_monthly_detrended <-
  full_join(pco2_product_monthly,
            pco2_product_monthly_regression %>% select(-c(value, resid, time))) %>%
  mutate(resid = value - fit)

pco2_product_monthly_detrended %>% 
  filter(biome %in% "Global") %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted monthly mean | Global") +
  facet_wrap(
    name ~ .,
    scales = "free_y",
    labeller = labeller(name = x_axis_labels),
    strip.position = "left",
    ncol = 2
  ) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank()
  )

pco2_product_monthly_detrended %>% 
  filter(biome %in% key_biomes) %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted monthly mean | Selected biomes") +
  facet_grid(
    name ~ biome,
    scales = "free_y",
    labeller = labeller(name = x_axis_labels),
    switch = "y"
  ) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank()
  )

pco2_product_monthly_detrended %>% 
  filter(biome %in% super_biomes) %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted monthly mean | Selected super biomes") +
  facet_grid(
    name ~ biome,
    scales = "free_y",
    labeller = labeller(name = x_axis_labels),
    switch = "y"
  ) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank()
  )

pco2_product_monthly_detrended %>% 
  filter(biome %in% key_biomes) %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted monthly mean | Selected biomes") +
  facet_grid(
    name ~ biome,
    scales = "free_y",
    labeller = labeller(name = x_axis_labels),
    switch = "y"
  ) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank()
  )

pco2_product_monthly_detrended %>%
  filter(biome %in% key_biomes) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(month, resid, group = as.factor(year))) +
      geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
                aes(col = year)) +
      scale_color_grayC() +
      new_scale_color() +
      geom_path(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(col = as.factor(year)),
        linewidth = 1
      ) +
      scale_color_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
      labs(title = paste("Anomalies from predicted monthly mean |", .x$biome)) +
      facet_wrap(
        name ~ .,
        scales = "free_y",
        labeller = labeller(name = x_axis_labels),
        strip.position = "left",
        ncol = 2
      ) +
      theme(
        strip.text.y.left = element_markdown(),
        strip.placement = "outside",
        strip.background.y = element_blank(),
        axis.title.y = element_blank(),
        legend.title = element_blank()
      )
  )

pco2_product_monthly_detrended %>% 
  filter(biome %in% super_biomes) %>%
  ggplot(aes(month, resid, group = as.factor(year))) +
  geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
            aes(col = year)) +
  scale_color_grayC() +
  new_scale_color() +
  geom_path(data = . %>% filter(between(year, 2023-1, 2023)),
            aes(col = as.factor(year)),
            linewidth = 1) +
  scale_color_manual(values = c("orange", "red"),
                     guide = guide_legend(reverse = TRUE,
                                          order = 1)) +
  scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
  labs(title = "Anomalies from predicted monthly mean | Selected super biomes") +
  facet_grid(
    name ~ biome,
    scales = "free_y",
    labeller = labeller(name = x_axis_labels),
    switch = "y"
  ) +
  theme(
    strip.text.y.left = element_markdown(),
    strip.placement = "outside",
    strip.background.y = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank()
  )

pco2_product_monthly_detrended %>%
  filter(biome %in% super_biomes) %>%
  group_split(biome) %>%
  # head(1) %>%
  map(
    ~ ggplot(data = .x,
             aes(month, resid, group = as.factor(year))) +
      geom_path(data = . %>% filter(!between(year, 2023-1, 2023)),
                aes(col = year)) +
      scale_color_grayC() +
      new_scale_color() +
      geom_path(
        data = . %>% filter(between(year, 2023-1, 2023)),
        aes(col = as.factor(year)),
        linewidth = 1
      ) +
      scale_color_manual(
        values = c("orange", "red"),
        guide = guide_legend(reverse = TRUE,
                             order = 1)
      ) +
      scale_x_continuous(breaks = seq(1, 12, 3), expand = c(0, 0)) +
      labs(title = paste("Anomalies from predicted monthly mean |", .x$biome)) +
      facet_wrap(
        name ~ .,
        scales = "free_y",
        labeller = labeller(name = x_axis_labels),
        strip.position = "left",
        ncol = 2
      ) +
      theme(
        strip.text.y.left = element_markdown(),
        strip.placement = "outside",
        strip.background.y = element_blank(),
        axis.title.y = element_blank(),
        legend.title = element_blank()
      )
  )

pco2_product_monthly_detrended %>%
  write_csv(paste0("../data/","CMEMS","_","2023","_biome_monthly_detrended.csv"))

pco2_product_annual_regression %>%
  filter(biome == "Global") %>%
  select(-c(value, fit)) %>% 
  pivot_wider(values_from = resid) %>% 
  pivot_longer(-c(year, biome, fgco2_int))  %>%
  ggplot(aes(value, fgco2_int)) +
  geom_hline(yintercept = 0) +
  geom_point(data = . %>% filter(!between(year, 2023-1, 2023)),
             aes(fill = year),
             shape = 21) +
  geom_smooth(
    data = . %>% filter(!between(year, 2023-1, 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(title = "Globally integrated fluxes",
       y = labels_breaks("fgco2_int")$i_legend_title) +
  facet_wrap(
    ~ name,
    scales = "free_x",
    labeller = labeller(name = x_axis_labels),
    strip.position = "bottom",
    ncol = 2
  ) +
  theme(
    strip.text.x.bottom = element_markdown(),
    strip.placement = "outside",
    strip.background.x = element_blank(),
    axis.title.y = element_markdown(),
    axis.title.x = element_blank(),
    legend.title = element_blank()
  )

pco2_product_monthly_detrended_anomaly <-
  pco2_product_monthly_detrended %>%
  select(year, month, biome, name, resid) %>%
  pivot_wider(names_from = name,
              values_from = resid)


pco2_product_monthly_detrended_anomaly %>%
  filter(biome == "Global") %>%
  pivot_longer(-c(year, month, biome, fgco2_int))  %>%
  ggplot(aes(value, fgco2_int)) +
  geom_hline(yintercept = 0) +
  geom_point(data = . %>% filter(year != 2023),
             aes(col = paste(min(year), max(year), sep = "-")),
             alpha = 0.2) +
  geom_smooth(
    data = . %>% filter(year != 2023),
    aes(col = paste(min(year), max(year), sep = "-")),
    method = "lm",
    se = FALSE,
    fullrange = TRUE
  )  +
  scale_color_grey(name = "") +
  new_scale_color() +
  geom_path(data = . %>% filter(year == 2023),
            aes(col = as.factor(month), group = 1))  +
  geom_point(data = . %>% filter(year == 2023),
             aes(fill =  as.factor(month)),
             shape = 21,
             size = 3)  +
  scale_color_scico_d(palette = "buda",
                     guide = guide_legend(reverse = TRUE,
                                          order = 1),
                     name = paste("Month\nof", 2023)) +
  scale_fill_scico_d(palette = "buda",
                     guide = guide_legend(reverse = TRUE,
                                          order = 1),
                     name = paste("Month\nof", 2023)) +
  labs(title = "Globally integrated fluxes",
       y = labels_breaks("fgco2_int")$i_legend_title) +
  facet_wrap(
    ~ name,
    scales = "free_x",
    labeller = labeller(name = x_axis_labels),
    strip.position = "bottom",
    ncol = 2
  ) +
  theme(
    strip.text.x.bottom = element_markdown(),
    strip.placement = "outside",
    strip.background.x = element_blank(),
    axis.title.y = element_markdown(),
    axis.title.x = element_blank()
  )

pco2_product_monthly_detrended_anomaly %>%
  filter(!(biome %in% c(super_biomes, "Global"))) %>%
  pivot_longer(-c(year, month, 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(col = paste(min(year), max(year), sep = "-")),
        alpha = 0.2
      ) +
      geom_smooth(
        data = . %>% filter(year != 2023),
        aes(col = paste(min(year), max(year), sep = "-")),
        method = "lm",
        se = FALSE,
        fullrange = TRUE
      )  +
      scale_color_grey(name = "") +
      new_scale_color() +
      geom_path(data = . %>% filter(year == 2023),
                aes(col = as.factor(month), group = 1))  +
      geom_point(
        data = . %>% filter(year == 2023),
        aes(fill = as.factor(month)),
        shape = 21,
        size = 3
      )  +
      scale_color_scico_d(
        palette = "buda",
        guide = guide_legend(reverse = TRUE,
                             order = 1),
        name = paste("Month\nof", 2023)
      ) +
      scale_fill_scico_d(
        palette = "buda",
        guide = guide_legend(reverse = TRUE,
                             order = 1),
        name = paste("Month\nof", 2023)
      ) +
      facet_wrap( ~ biome, ncol = 3, scales = "free") +
      labs(
        title = "Biome integrated fluxes",
        y = labels_breaks("fgco2_int")$i_legend_title,
        x = labels_breaks(.x %>% distinct(name))$i_legend_title
      ) +
      theme(axis.title.x = element_markdown(),
            axis.title.y = element_markdown())
  )