library(ggplot2)
library(dplyr)
library(sf)
library(terra)

# Define the file paths
raster_dir <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/laravel_app/storage/app/chemba/merged_final_tif"
polygon_path <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/r_app/experiments/pivot.geojson"

# Load the polygon
polygon <- st_read(polygon_path)

# Filter the polygon for the specific field
polygon_filtered <- subset(polygon, field == "1.1")

# List all raster files in the directory
raster_files <- list.files(raster_dir, pattern = "\\.tif$", full.names = FALSE)
raster_file <- raster_files[1]

# Initialize an empty list to store CI values
ci_values <- list()

# Define a function to process each raster file
process_raster_file <- function(raster_file) {
    # Extract the date from the file name (assuming the date is in YYYYMMDD format)
    date <- as.Date(substr(raster_file, 1, 10), format = "%Y-%m-%d")
    
    # Load the raster
    raster_data <- terra::rast(file.path(raster_dir, raster_file))
    
    # Crop the raster using the filtered polygon
    cropped_raster <- terra::crop(raster_data, polygon_filtered)
    
    # Extract the CI band (assuming the CI band is the first band)
    ci_band <- cropped_raster$CI
    
    # Extract the CI values from the CI band
    return(data.frame(Date = date, CI = terra::values(ci_band)) %>% 
        filter(!is.na(CI)) %>% 
        mutate(Polygon = unique(polygon_filtered$field)))
}

# Use walk to apply the function to each raster file
ci_values <- purrr::map(raster_files, process_raster_file)

# Combine all CI values into a single data frame
ci_values_df <- do.call(rbind, ci_values)

head(ci_values_df)

# Plot the first raster
plot(terra::rast(file.path(raster_dir, raster_files[1]))[[1]])

# Calculate the mean CI value
ci_stats <- ci_values_df %>%
    group_by(Date, Polygon) %>%
    summarize(
        mean_ci = mean(CI, na.rm = TRUE),
        q5 = quantile(CI, 0.05, na.rm = TRUE),
        q25 = quantile(CI, 0.25, na.rm = TRUE),
        q50 = quantile(CI, 0.50, na.rm = TRUE),
        q75 = quantile(CI, 0.75, na.rm = TRUE),
        q95 = quantile(CI, 0.95, na.rm = TRUE),
        cv = sd(CI, na.rm = TRUE) / mean(CI, na.rm = TRUE) * 100,
        .groups = "drop"
    )

# Plot the mean CI value over time
ggplot(ci_stats, aes(x = Date, y = mean_ci)) +
    geom_line() +
    geom_ribbon(aes(ymin = q25, ymax = q75), alpha = 0.2) +
    labs(title = "Mean CI value over time",
         x = "Date",
         y = "Mean CI") +
    theme_minimal()

# Plot the coefficient of variation over time
ggplot(ci_stats, aes(x = Date, y = cv)) +
    geom_line() +
    labs(title = "Coefficient of variation over time",
         x = "Date",
         y = "CV (%)") +
    theme_minimal()