library(sf)
library(terra)
library(tidyverse)
library(lubridate)
library(exactextractr)

# Vang alle command line argumenten op
args <- commandArgs(trailingOnly = TRUE)

# Converteer het tweede argument naar een string waarde
project_dir <- as.character(args[1])

# Controleer of data_dir een geldige waarde is
if (!is.character(project_dir)) {
  project_dir <- "chemba"
}

source("parameters_project.R")
source("ci_extraction_utils.R")

pivot_stats2 <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) %>%
    ungroup() %>%
    group_by(field, sub_field) %>%
    summarise(across(everything(), ~ first(na.omit(.))), .groups = "drop")

   #%>% drop_na(pivot_quadrant)

# gather data into long format for easier calculation and visualisation
pivot_stats_long <- pivot_stats2 %>%
  tidyr::gather("Date", value, -field, -sub_field ) %>%
  mutate(#Date = right(Date, 8),
         Date = lubridate::ymd(Date)
  ) %>%
  drop_na(c("value","Date")) %>%
  mutate(value = as.numeric(value))%>%
  filter_all(all_vars(!is.infinite(.))) %>%
  # rename(field = pivot_quadrant,
  #        sub_field = field) %>%
  unique()

years <- harvesting_data %>% 
  filter(!is.na(season_start)) %>% 
  distinct(year) %>% 
  pull(year)

extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
  # Filter harvesting data for the given season and field names
  filtered_harvesting_data <- harvesting_data %>%
    filter(year == season, sub_field %in% field_names)
  
  # Filter field CI data for the given field names
  filtered_field_CI_data <- field_CI_data %>%
    filter(sub_field %in% field_names)
  
  # Return an empty data frame if no CI data is found
  if (nrow(filtered_field_CI_data) == 0) {
    return(data.frame())
  }
  
  # Create a linear interpolation function for the CI data
  ApproxFun <- approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value)
  Dates <- seq.Date(ymd(min(filtered_field_CI_data$Date)), ymd(max(filtered_field_CI_data$Date)), by = 1)
  LinearFit <- ApproxFun(Dates)
  
  # Combine interpolated data with the original CI data
  CI <- data.frame(Date = Dates, FitData = LinearFit) %>%
    left_join(., filtered_field_CI_data, by = "Date") %>%
    filter(Date > filtered_harvesting_data$season_start & Date < filtered_harvesting_data$season_end)
  
  # If CI is empty after filtering, return an empty dataframe
  if (nrow(CI) == 0) {
    return(data.frame())
  }
  
  # Add additional columns if data exists
  CI <- CI %>%
    mutate(DOY = seq(1, n(), 1),
           model = paste0("Data", season, " : ", field_names),
           season = season,
           subField = field_names)
  
  return(CI)
}



CI_all <- map_df(years, function(yr) {
  # yr = 2021
  message(yr)
  
  # Get the fields harvested in this year
  sub_fields <- harvesting_data %>%
    filter(year == yr) %>%
    filter(!is.na(season_start)) %>%
    pull(sub_field)
  
  # Filter sub_fields to only include those with value data in pivot_stats_long
  valid_sub_fields <- sub_fields %>%
    keep(~ any(pivot_stats_long$sub_field == .x))
  
  # Extract data for each valid field
  map(valid_sub_fields, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = yr)) %>%
    list_rbind()
})

CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitData),
                                                cumulative_CI = cumsum(FitData))


message('CI_all cumulative')

head(CI_all)
message('show head')

saveRDS(CI_all, here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))
message('rds saved')