diff --git a/laravel_app/app/Jobs/ProjectDownloadRDSJob.php b/laravel_app/app/Jobs/ProjectDownloadRDSJob.php
index 688f77c..7c27b8f 100644
--- a/laravel_app/app/Jobs/ProjectDownloadRDSJob.php
+++ b/laravel_app/app/Jobs/ProjectDownloadRDSJob.php
@@ -34,7 +34,7 @@ public function handle(): void
     {
 
         $command = [
-            sprintf('%sbuild_RDS.sh', base_path('../')),
+            sprintf('%supdate_RDS.sh', base_path('../')),
             sprintf('--end_date=%s', $this->date->format('Y-m-d')),
             sprintf('--offset=%d', $this->offset),
             sprintf('--project_dir=%s', $this->project->download_path),
diff --git a/r_app/2_CI_data_prep.R b/r_app/2_CI_data_prep.R
deleted file mode 100644
index 30b436a..0000000
--- a/r_app/2_CI_data_prep.R
+++ /dev/null
@@ -1,855 +0,0 @@
-# activeer de renv omgeving;
-
-# renv::activate('~/smartCane/r_app')
-# renv::restore()
-
-
-
-library(here)
-library(sf)
-library(terra)
-library(tidyverse)
-library(lubridate)
-library(exactextractr)
-library(readxl)
-#funcion CI_prep package
-
-date_list <- function(end_date, offset){
-  offset <- as.numeric(offset) - 1
-  end_date <- as.Date(end_date)
-  start_date <- end_date - lubridate::days(offset)
-
-  week <- week(start_date)
-  year <- year(start_date)
-  days_filter <- seq(from = start_date, to = end_date, by = "day")
-
-  return(list("week" = week, "year" = year, "days_filter" = days_filter))
-}
-
-# date_list <- function(weeks_in_the_paste){
-#   week <- week(Sys.Date()- weeks(weeks_in_the_paste) )
-#   year <- year(Sys.Date()- weeks(weeks_in_the_paste) )
-#   days_filter <- Sys.Date() - weeks(weeks_in_the_paste) - days(0:6)
-#
-#   return(c("week" = week,
-#            "year" = year,
-#            "days_filter" = list(days_filter)))
-#
-# }
-
-CI_func <- function(x, drop_layers = FALSE){
-  CI <- x[[4]]/x[[2]]-1
-  add(x) <- CI
-  names(x) <- c("red", "green", "blue","nir", "cloud" ,"CI")
-  if(drop_layers == FALSE){
-    return(x)
-  }else{
-    return(x$CI)
-  }
-}
-
-mask_raster <- function(raster, fields){
-  # x <- rast(filtered_files[1])
-  x <- rast(raster) 
-  emtpy_or_full <- global(x, sum)
-  
-  if(emtpy_or_full[1,] >=  2000000){
-    names(x) <- c("red", "green", "blue","nir", "cloud")
-    cloud <- x$cloud
-    
-    cloud[cloud == 0 ] <- NA
-    x_masked <- mask(x, cloud, inverse = T) %>% crop(.,fields, mask = TRUE )
-    x_masked <-  x_masked %>% CI_func()
-    
-    message(raster, " masked")
-    return(x_masked)
-  }
-}
-
-date_extract <- function(file_path) {
-  str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}")
-}
-
-extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) {
-  # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE)
-  # date <- date_extract(filtered_files[1])
-  # field_geojson <- sf::st_as_sf(pivot_sf_q)
-  
-  field_geojson <- sf::st_as_sf(field_geojson)
-  x <- rast(file[1]) %>% CI_func(drop_layers = TRUE)
-  date <- date_extract(file)
-  
-  pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x, field_geojson, fun = "mean"), 2)) %>%
-    st_drop_geometry() %>% rename("{date}" := mean_CI) 
-  
-  save_suffix <- if (quadrants){"quadrant"} else {"whole_field"}
-  
-  save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
-  
-  saveRDS(pivot_stats, save_path)
-}
-
-right = function(text, num_char) {
-  substr(text, nchar(text) - (num_char-1), nchar(text))
-}
-
-extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
-  #  field_names =  "1.2A"
-  #  harvesting_data = harvesting_data
-  # field_CI_data = pivot_stats_long
-  #  season= 2021
-  
-  filtered_harvesting_data <- harvesting_data %>%
-    filter(year == season, field %in% field_names)
-  
-  filtered_field_CI_data <- field_CI_data %>%
-    filter(field %in% field_names)
-  
-  # CI <- map_df(field_names, ~ {
-  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)
-  
-  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) %>%
-    mutate(DOY = seq(1, n(), 1),
-           model = paste0("Data", season, " : ", field_names),
-           season = season,
-           field = field_names)
-  # }) #%>%
-  #{if (length(field_names) > 0) message("Done!")}
-  
-  return(CI)
-}
-
-load_fields <- function(geojson_path) {
-  field_geojson <- st_read(geojson_path) %>%
-    select(pivot, pivot_quadrant) %>%
-    vect()
-  return(field_geojson)
-}
-
-load_harvest_data <- function(havest_data_path){
-  harvest_data <- readRDS(havest_data_path)
-  return(harvest_data)
-}
-
-load_rasters <- function(raster_path, dates) {
-  raster_files <- list.files(raster_path, full.names = TRUE, pattern = ".tif")
-  filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
-    compact() %>%
-    flatten_chr()
-  
-  return(filtered_files)
-}
-
-mask_and_set_names <- function(filtered_files, fields) {
-  rasters_masked <- map(filtered_files, mask_raster, fields = fields) %>% set_names(filtered_files)
-  rasters_masked[sapply(rasters_masked, is.null)] <- NULL
-  rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract))
-  
-  return(rasters_masked)
-}
-
-calculate_total_pix_area <- function(filtered_files, fields_geojson) {
-  # total_pix_area <- rast(filtered_files[1]) %>%
-  #   subset(1) %>%
-  #   crop(fields_geojson, mask = TRUE)%>%
-  #   global(.data, fun = "notNA")
-  total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(fields_geojson, mask = TRUE) %>% freq(., usenames = TRUE)
-  
-  return(total_pix_area)
-}
-
-cloud_layer_extract <- function(rasters_masked){
-  cloud_layer_rast <- map(rasters_masked, function(spatraster) {
-    spatraster[[5]]
-  }) %>% rast()
-  
-  return(cloud_layer_rast)
-}
-
-calculate_cloud_coverage <- function(cloud_layer_rast, total_pix_area) {
-  cloud_perc_list <- freq(cloud_layer_rast, usenames = TRUE) %>%
-    mutate(cloud_perc = (100 -((count/total_pix_area$count)*100)),
-           cloud_thres_5perc = as.integer(cloud_perc < 5),
-           cloud_thres_40perc = as.integer(cloud_perc < 40)) %>%
-    rename(Date = layer) %>% select(-value, -count)
-  
-  cloud_index_5perc <- which(cloud_perc_list$cloud_thres_5perc == max(cloud_perc_list$cloud_thres_5perc))
-  cloud_index_40perc <- which(cloud_perc_list$cloud_thres_40perc == max(cloud_perc_list$cloud_thres_40perc))
-  
-  return(list(cloud_perc_list = cloud_perc_list, cloud_index_5perc = cloud_index_5perc, cloud_index_40perc = cloud_index_40perc))
-}
-
-process_cloud_coverage <- function(cloud_coverage, rasters_masked) {
-  if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) > 1) {
-    message("More than 1 raster without clouds (<5%), max mosaic made ")
-    
-    cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_5perc]
-    rsrc <- sprc(cloudy_rasters_list)
-    x <- mosaic(rsrc)
-    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
-    
-  } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) == 1) {
-    message("Only 1 raster without clouds (<5%)")
-    
-    x <- rast(rasters_masked[cloud_coverage$cloud_index_5perc[1]])
-    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
-    
-  } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) > 1) {
-    message("More than 1 image contains clouds, composite made of <40% cloud cover images")
-    
-    cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_40perc]
-    rsrc <- sprc(cloudy_rasters_list)
-    x <- mosaic(rsrc)
-    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
-    
-  } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) == 0) {
-    message("No cloud free images available")
-    x <-  rast(rasters_masked[1])
-    x[x] <- NA
-    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
-  }
-  
-  return(x)
-}
-
-extract_rasters_daily_func <- function(daily_vals_dir, filtered_files, fields_geojson) {
-  extracted_files <- walk(filtered_files, extract_rasters_daily, field_geojson = fields_geojson, quadrants = TRUE, daily_vals_dir)
-  return(extracted_files)
-}
-
-CI_load <- function(daily_vals_dir, grouping_variable){
-  extracted_values <-  list.files(here(daily_vals_dir), full.names = TRUE)
-  
-  field_CI_values <- extracted_values %>%
-    map_dfr(readRDS) %>%
-    group_by(.data[[grouping_variable]]) %>%
-    summarise(across(everything(), ~ first(na.omit(.))))
-  return(field_CI_values)
-}
-
-CI_long <- function(field_CI_values, pivot_long_cols){
-  field_CI_long <- field_CI_values %>%
-    gather("Date", value, -pivot_long_cols) %>%
-    mutate(Date = right(Date, 8),
-           Date = ymd(Date)
-    ) %>%
-    drop_na(c("value","Date")) %>%
-    mutate(value = as.numeric(value))%>%
-    filter_all(all_vars(!is.infinite(.)))%>%
-    rename(field = pivot_quadrant)
-  
-  return(field_CI_long)
-}
-
-process_year_data <- function(year, harvest_data, field_CI_long) {
-  pivots_dates_year <- harvest_data %>% na.omit() %>% filter(year == year)
-  pivot_select_model_year <- unique(pivots_dates_year$field)
-  
-  data <- map_dfr(pivot_select_model_year, ~ extract_CI_data(.x, harvest_data, field_CI_long, season = year))
-  
-  return(data)
-}
-
-
-
-#functions for CI_data_prep
-create_mask_and_crop <- function(file, pivot_sf_q) {
-  # file <- filtered_files[5]
-  message("starting ", file)
-  loaded_raster <- rast(file)
-  names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
-  #   names(CI) <- c("green","nir")
-  message("raster loaded")
-  
-  # CI <- CI[[2]] / CI[[1]] - 1
-  CI <- loaded_raster$NIR / loaded_raster$Green - 1
-  
-  loaded_raster$CI <- CI
-  #   CI <- CI$nir/CI$green-1
-  message("CI calculated")
-  loaded_raster <- terra::crop(loaded_raster, pivot_sf_q, mask = TRUE) #%>% CI_func()
-  
-  loaded_raster[loaded_raster == 0] <- NA
-  # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI")
-  # v_crop$CI <- v_crop$CI - 1
-  new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
-  writeRaster(loaded_raster, new_file, overwrite = TRUE)
-  
-  vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
-  terra::vrt(new_file, vrt_file, overwrite = TRUE)
-  
-  # v_crop <- mask_raster(v, pivot_sf_q)
-  return(loaded_raster)
-}
-
-
-extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) {
-  # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE)
-  # date <- date_extract(filtered_files[1])
-  # field_geojson <- sf::st_as_sf(pivot_sf_q)
-  field_geojson <- sf::st_as_sf(field_geojson)
-  x <- rast(file[1])
-  date <- date_extract(file)
-  pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>%
-    st_drop_geometry() %>% rename("{date}" := mean_CI)
-  save_suffix <- if (quadrants){"quadrant"} else {"whole_field"}
-  save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
-  saveRDS(pivot_stats, save_path)
-}
-
-#functions for rmarkdown file
-
-
-create_RGB_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = TRUE, legend_is_portrait = FALSE, week, age, red = TRUE) {
-  r <- if (red) 1 else 4  # Set r based on the value of red
-  title <- if (red) paste0("RGB image of the fields") else paste0("False colour image of the fields")
-  
-  tm_shape(pivot_raster, unit = "m") + tm_rgb(r = r, g = 2, b = 3, max.value = 255) +
-    tm_layout(main.title = title, 
-              main.title.size = 1) +
-    tm_scale_bar(position = c("right", "top"), text.color = "black") +
-    tm_compass(position = c("right", "top"), text.color = "black") +
-    tm_shape(pivot_shape) + tm_borders(col = "gray") +
-    tm_text("sub_field", size = 0.6, col = "gray") +
-    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha = 0.5)
-}
-
-create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, legend_only = F){
-  tm_shape(pivot_raster, unit = "m")+
-    tm_raster(breaks =  CI_breaks, palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) +
-    tm_layout(main.title = paste0("Max CI week ", week,"\n",  age, " weeks old"),
-              main.title.size = 1, legend.show = show_legend, legend.only = legend_only) +
-    tm_shape(pivot_shape) +
-    tm_borders(lwd = 3) +  tm_text("sub_field", size = 1/2) +
-    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5) +tmap_options(check.and.fix = TRUE)
-}
-
-create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age){
-  tm_shape(pivot_raster, unit = "m")+
-    tm_raster(breaks = CI_diff_breaks, palette = "PRGn",legend.is.portrait = legend_is_portrait ,midpoint = 0, title = "CI difference") +
-    tm_layout(main.title =  paste0("CI change week ", week_1, "- week ",week_2, "\n", age," weeks old"),
-              main.title.size = 1, legend.show = show_legend) +
-    tm_shape(pivot_shape) +
-    tm_borders(lwd = 3) +  tm_text("sub_field", size = 1/2) +
-    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
-}
-
-ci_plot <- function(pivotName){
-  # pivotName = "MV2B09"
-  # pivotName = "1.1"
-  
-  pivotShape <- AllPivots_merged %>% terra::subset(field %in% pivotName)  %>% st_transform(crs(CI))
-  # age <- AllPivots %>% dplyr::filter(field %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique() %>%
-  #   mutate(Age = round(Age))
-  
-  age <- AllPivots %>%
-    group_by(field) %>%
-    filter(Season == max(Season, na.rm = TRUE), field %in% pivotName) %>%
-    dplyr::select(Age)%>% st_drop_geometry()  %>% unique() 
-  
-  AllPivots2 <- AllPivots0 %>% dplyr::filter(field %in% pivotName)
-  
-  singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape)
-  singlePivot_m1 <- CI_m1 %>% crop(., pivotShape) %>% mask(., pivotShape)
-  singlePivot_m2 <- CI_m2 %>% crop(., pivotShape) %>% mask(., pivotShape)
-  # singlePivot_m3 <- CI_m3 %>% crop(., pivotShape) %>% mask(., pivotShape)
-  
-  singlePivot_RGB <- RGB_raster %>% crop(., pivotShape) %>% mask(., pivotShape)
-  singlePivot_false <- RGB_raster_stretch %>% crop(., pivotShape) %>% mask(., pivotShape)
-  
-  abs_CI_last_week <- last_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape)
-  abs_CI_three_week <- three_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape)
-  
-  # planting_date <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% ungroup()  %>% dplyr::select(planting_date)  %>% unique()
-  
-  joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape))  %>%  dplyr::filter(field %in% pivotName) %>% st_crop(., pivotShape)
-  
-  # CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait = T, week = week_minus_2, age = age  -2)
-  Legend_map <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age  -1, legend_only = T)
-  CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age  -1)
-  CImap <- create_CI_map(singlePivot, AllPivots2,  joined_spans2, show_legend= F, legend_is_portrait = T, week = week, age = age )
-  RGBmap <- create_RGB_map(singlePivot_false, AllPivots2,  joined_spans2, show_legend= F,  week = week, age = age, red =T )
-  Falsemap <- create_RGB_map(singlePivot_false, AllPivots2,  joined_spans2, show_legend= F,  week = week, age = age, red =F )
-  
-  
-  CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week,AllPivots2,  joined_spans2,  show_legend = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_1, age = age)
-  CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2,  joined_spans2, show_legend = F, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age)
-  
-  # tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1)
-  tst <- tmap_arrange(RGBmap,Falsemap,
-                      CImap_m1, CImap, 
-                      CI_max_abs_last_week, CI_max_abs_three_week, 
-                      ncol = 2)
-  
-  cat(paste("## field", pivotName, "-", age$Age[1], "weeks after planting/harvest", "\n"))
-  # cat("\n")
-  # cat('<h2> Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest <h2>')
-  # cat(paste("# Pivot",pivots$pivot[i],"\n"))
-  print(tst)
-  
-}
-
-
-subchunkify <- function(g, fig_height=7, fig_width=5) {
-  g_deparsed <- paste0(deparse(
-    function() {g}
-  ), collapse = '')
-  
-  sub_chunk <- paste0("```{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", echo=FALSE}",
-                      "\n(",
-                      g_deparsed
-                      , ")()",
-                      "\n```
-  ")
-  
-  cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE))
-}
-
-cum_ci_plot <- function(pivotName){
-  
-  # pivotName = "2.1"
-  
-  # Check if pivotName exists in the data
-  if (!pivotName %in% unique(CI_quadrant$field)) {
-    # message("PivotName '", pivotName, "' not found. Plotting empty graph.")
-    g <- ggplot() +   labs(title = "Empty Graph - Yield dates missing")
-    return(
-      subchunkify(g, fig_height=2, fig_width = 10)
-    )
-  } else {
-    # message("PivotName '", pivotName, "'  found. Plotting normal graph.")
-    data_ci <- CI_quadrant %>% filter(field %in% pivotName)
-    
-    
-    
-    data_ci2 <- data_ci  %>% mutate(CI_rate = cumulative_CI/DOY,
-                                    week = week(Date))%>% group_by(sub_field) %>%
-      mutate(mean_rolling10 = rollapplyr(CI_rate  , width = 10, FUN = mean, partial = TRUE))
-    
-    # date_preperation_perfect_pivot <- data_ci2 %>% group_by(season)  %>% summarise(min_date = min(Date),
-    #                                                                                max_date = max(Date),
-    #                                                                                days = max_date - min_date)
-    
-    # Identify unique seasons
-    filtered_data <- data_ci2 %>%
-      group_by(season) %>%
-      mutate(rank = dense_rank(desc(season))) %>%
-      filter(rank <= 2) %>%
-      ungroup() %>%
-      dplyr::select(-rank)
-    
-    
-    # g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) +
-    g <- ggplot(data= filtered_data ) +
-      # geom_line(aes(Date, mean_rolling10, col = sub_field)) +
-      geom_line(aes(Date, CI_rate, col = sub_field)) +
-      facet_wrap(~season, scales = "free_x") +
-      # geom_line(data= perfect_pivot, aes(Date , mean_rolling10, col = "Model CI (p5.1 Data 2022, \n date x axis is fictive)"), lty="11",size=1) +
-      labs(title = paste("CI rate - field", pivotName),
-           y = "CI rate (cumulative CI / Age)")+
-      # scale_y_continuous(limits=c(0.5,3), breaks = seq(0.5, 3, 0.5))+
-      scale_x_date(date_breaks = "1 month", date_labels =  "%m-%Y") +
-      theme(axis.text.x = element_text(angle = 60, hjust = 1),
-            legend.justification=c(1,0), legend.position = c(1, 0),
-            legend.title = element_text(size = 8),
-            legend.text = element_text(size = 8)) +
-      guides(color = guide_legend(nrow = 2, byrow = TRUE))
-    subchunkify(g, fig_height=6, fig_width = 10)
-  }
-}
-
-
-
-# Vang alle command line argumenten op
-args <- commandArgs(trailingOnly = TRUE)
-
-# Controleer of er ten minste één argument is doorgegeven
-if (length(args) == 0) {
-  stop("Geen argumenten doorgegeven aan het script")
-}
-
-# Converteer het eerste argument naar een numerieke waarde
-end_date <- as.Date(args[1])
-
-
-offset <- as.numeric(args[2])
-# Controleer of weeks_ago een geldig getal is
-if (is.na(offset)) {
-  # stop("Het argument is geen geldig getal")
-   offset <- 7
-}
-
-# Converteer het tweede argument naar een string waarde
-project_dir <- as.character(args[3])
-
-# Controleer of data_dir een geldige waarde is
-if (!is.character(project_dir)) {
-  project_dir <- "chemba"
-}
-
-
-laravel_storage_dir <- here("laravel_app/storage/app", project_dir)
-#preparing directories
-planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
-merged_final <- here(laravel_storage_dir, "merged_final_tif")
-
-new_project_question = FALSE
-planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
-merged_final <- here(laravel_storage_dir, "merged_final_tif")
-
-data_dir <- here(laravel_storage_dir, "Data")
-extracted_CI_dir <- here(data_dir, "extracted_ci")
-daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals")
-cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals")
-
-weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
-daily_vrt <- here(data_dir, "vrt")
-harvest_dir <-  here(data_dir, "HarvestData")
-
-source("parameters_project.R")
-
-dir.create(here(laravel_storage_dir))
-dir.create(here(data_dir))
-dir.create(here(extracted_CI_dir))
-dir.create(here(daily_CI_vals_dir))
-dir.create(here(cumulative_CI_vals_dir))
-dir.create(here(weekly_CI_mosaic))
-dir.create(here(daily_vrt))
-dir.create(merged_final)
-dir.create(harvest_dir)
-
-week <- week(end_date)
-
-
-#weeks_ago = 0
-# Creating weekly mosaic
-#dates <- date_list(weeks_ago)
-dates <- date_list(end_date, offset)
-print(dates)
-#load pivot geojson
-# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
-
-raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
-
-filtered_files  <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
-  compact() %>%
-  flatten_chr()
-head(filtered_files)
-
-# filtered_files <- raster_files #for first CI extraction
-
-# create_mask_and_crop <- function(file, field_boundaries) {
-# message("starting ", file)
-#   CI <- rast(file)
-# #   names(CI) <- c("green","nir")
-# message("raster loaded")
-# 
-#   CI <- CI[[2]]/CI[[1]]-1
-# #   CI <- CI$nir/CI$green-1
-# message("CI calculated")
-#   CI <- terra::crop(CI, field_boundaries, mask = TRUE) #%>% CI_func()
-# 
-#   new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
-#   writeRaster(CI, new_file, overwrite = TRUE)
-# 
-#   vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
-#   terra::vrt(new_file, vrt_file, overwrite = TRUE)
-# 
-#   return(CI)
-# }
-create_mask_and_crop <- function(file, field_boundaries) {
-  # file <- filtered_files[5]
-  message("starting ", file)
-  loaded_raster <- rast(file)
-  names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
-  #   names(CI) <- c("green","nir")
-  message("raster loaded")
-  
-  # CI <- CI[[2]] / CI[[1]] - 1
-  CI <- loaded_raster$NIR / loaded_raster$Green - 1
-  
-  loaded_raster$CI <- CI
-  #   CI <- CI$nir/CI$green-1
-  message("CI calculated")
-  loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE) #%>% CI_func()
-  
-  loaded_raster[loaded_raster == 0] <- NA
-  # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI")
-  # v_crop$CI <- v_crop$CI - 1
-  new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
-  writeRaster(loaded_raster, new_file, overwrite = TRUE)
-  
-  vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
-  terra::vrt(new_file, vrt_file, overwrite = TRUE)
-  
-  # v_crop <- mask_raster(v, pivot_sf_q)
-  return(loaded_raster)
-}
-
-vrt_list <- list()
-
-# for (file in raster_files) {
-#   v_crop <- create_mask_and_crop(file, field_boundaries)
-#   message(file, " processed")
-#   gc()
-# }
-
-for (file in filtered_files) {
-  v_crop <- create_mask_and_crop(file, field_boundaries)
-  emtpy_or_full <- global(v_crop, "notNA")
-
-  vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
-
-  if(emtpy_or_full[1,] >  10000){
-    vrt_list[vrt_file] <- vrt_file
-
-  }else{
-    file.remove(vrt_file)
-
-  }
-
-  message(file, " processed")
-  gc()
-}
-
-# list_global <- list_global %>% flatten_chr()
-vrt_list <- vrt_list %>% flatten_chr()
-
-total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>%
-  crop(field_boundaries, mask = TRUE) %>%
-  global(., fun="notNA") #%>%
-
-layer_5_list <- purrr::map(vrt_list, function(vrt_list) {
-  rast(vrt_list[1]) %>% terra::subset(1)
-}) %>% rast()
-
-missing_pixels_count <- layer_5_list %>% global(., fun="notNA") %>%
-  mutate(
-    total_pixels = total_pix_area$notNA,
-    missing_pixels_percentage = round(100 -((notNA/total_pix_area$notNA)*100)),
-    thres_5perc = as.integer(missing_pixels_percentage < 5),
-    thres_40perc = as.integer(missing_pixels_percentage < 45)
-    )
-
-index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc)  )
-index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc))
-
-## Create mosaic
-
-if(sum(missing_pixels_count$thres_5perc)>1){
-  message("More than 1 raster without clouds (<5%), max composite made")
-  
-  cloudy_rasters_list <- vrt_list[index_5perc]
-  
-  rsrc <- sprc(cloudy_rasters_list)
-  x <- mosaic(rsrc, fun = "max")
-  
-  # names(x) <-  "CI"
-  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
-}else if(sum(missing_pixels_count$thres_5perc)==1){
-  message("Only 1 raster without clouds (<5%)")
-  
-  x <- rast(vrt_list[index_5perc[1]])
-  # names(x) <-  c("CI")
-  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
-}else if(sum(missing_pixels_count$thres_40perc)>1){
-  message("More than 1 image contains clouds, composite made of <40% cloud cover images")
-  
-  cloudy_rasters_list <- vrt_list[index_40perc]
-  
-  rsrc <- sprc(cloudy_rasters_list)
-  x <- mosaic(rsrc, fun = "max")
-  # names(x) <-  "CI"
-  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
-}else if(sum(missing_pixels_count$thres_40perc)==1){
-  message("Only 1 image available but contains clouds")
-  
-  x <- rast(vrt_list[index_40perc[1]])
-  # names(x) <-  c("CI")
-  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
-}else{
-  message("No cloud free images available, all images combined")
-
-  rsrc <- sprc(vrt_list)
-  x <- mosaic(rsrc, fun = "max")
-  # x <-  rast(vrt_list[1]) %>% setValues(NA)
-  # names(x) <-  c("CI")
-  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
-}
-
-plot(x$CI, main = paste("CI map ", dates$week))
-plotRGB(x, main = paste("RGB map ", dates$week))
-
-file_path_tif <- here(weekly_CI_mosaic ,paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif"))
-writeRaster(x, file_path_tif, overwrite=TRUE)
-message("Raster written/made at: ", file_path_tif)
-
-
-# Extracting CI
-extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) {
-  # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE)
-  # date <- date_extract(filtered_files[1])
-  # field_geojson <- sf::st_as_sf(pivot_sf_q)
-  field_geojson <- sf::st_as_sf(field_geojson)
-  x <- rast(file[1])
-  date <- date_extract(file)
-  pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>%
-    st_drop_geometry() %>% rename("{date}" := mean_CI)
-  save_suffix <- if (quadrants){"quadrant"} else {"whole_field"}
-  save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
-  saveRDS(pivot_stats, save_path)
-}
-
-# pivot_sf_q <- st_read(here("..", "Data", "pivot_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
-# pivot_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant)   %>% group_by(pivot) %>% summarise() %>% vect()
-# message("pivot loaded")
-
-raster_files_NEW <- list.files(merged_final,full.names = T, pattern = ".tif")
-
-# pivots_dates0 <- readRDS(here(harvest_dir, "harvest_data_new")) %>% filter(
-#   pivot %in% c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13",
-#                "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3",
-#                "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3")
-# )
-
-# harvesting_data <- pivots_dates0 %>%
-#   select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023", "season_start_2024", "season_end_2024")) %>%
-#   pivot_longer(cols = starts_with("season"), names_to = "Year", values_to = "value") %>%
-#   separate(Year, into = c("name", "Year"), sep = "(?<=season_start|season_end)\\_", remove = FALSE) %>%
-#   mutate(name = str_to_title(name)) %>%
-#   pivot_wider(names_from = name, values_from = value) %>%
-#   rename(field = pivot_quadrant)
-
-#If run for the firsttime, it will extract all data since the start and put into a table.rds. otherwise it will only add on to the existing table.
-
-
-    # Define the path to the file
-    file_path <- here(cumulative_CI_vals_dir,"combined_CI_data.rds")
-
-    # Check if the file exists
-    if (!file.exists(file_path)) {
-      # Create the file with columns "column1" and "column2"
-      data <- data.frame(sub_field=NA, field=NA)
-      saveRDS(data, file_path)
-    }
-
-  print("combined_CI_data.rds exists, adding the latest image data to the table.")
-
-  filtered_files  <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>%
-    compact() %>%
-    flatten_chr()
-
-  walk(filtered_files, extract_rasters_daily, field_geojson= field_boundaries, quadrants = TRUE, daily_CI_vals_dir)
-
-  extracted_values <-  list.files(daily_CI_vals_dir, full.names = TRUE)
-
-  extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>%
-    compact() %>%
-    flatten_chr()
-
-  pivot_stats <- extracted_values %>%
-    map(readRDS) %>% list_rbind() %>%
-    group_by(sub_field) %>%
-    summarise(across(everything(), ~ first(na.omit(.))))
-
-  combined_CI_data <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #%>% drop_na(pivot_quadrant)
-  pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data)
-  # pivot_stats2 <- combined_CI_data
-  print("All CI values extracted from latest 7 images.")
-  saveRDS(combined_CI_data, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file
-
-
-# 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()
-
-
-# #2021
-# pivot_select_model_Data_2021 <- harvesting_data %>% filter(Year == 2021) %>% pull(field)
-# 
-# pivot_select_model_Data_2022 <- harvesting_data %>% filter(Year == 2022) %>% pull(field)
-
-# pivot_select_model_Data_2023 <- harvesting_data %>% filter(year == 2023) %>% filter(!is.na(season_start))  %>% pull(sub_field)
-
-pivot_select_model_Data_2024 <- harvesting_data %>% filter(year == 2024)%>% filter(!is.na(season_start)) %>% pull(sub_field)
-
-# pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022))
-# pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant )
-# 
-# pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023))
-# pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant)
-# 
-# pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024))
-# pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant)
-
-extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
-  #  field_names =  "AG1D06P"
-  #  field_names =  "1.13A"
-  #  harvesting_data = harvesting_data
-  # field_CI_data = pivot_stats_long
-  #  season= 2023
-  
-  filtered_harvesting_data <- harvesting_data %>%
-    na.omit() %>%
-    filter(year == season, sub_field %in% field_names)
-  
-  filtered_field_CI_data <- field_CI_data %>%
-    filter(sub_field %in% field_names)
-  
-  # CI <- map_df(field_names, ~ {
-  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)
-  
-  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) %>%
-    mutate(DOY = seq(1, n(), 1),
-           model = paste0("Data", season, " : ", field_names),
-           season = season,
-           sub_field = field_names)
-  # }) #%>%
-  #{if (length(field_names) > 0) message("Done!")}
-  
-  return(CI)
-}
-
-## Extracting the correct CI values
-# Data_2021 <- map(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021)) %>% list_rbind()
-# message('2021')
-# Data_2022 <- map(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022)) %>% list_rbind()
-# message('2022')
-# Data_2023 <- map(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023)) %>% list_rbind()
-# message('2023')
-Data_2024 <- map(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024)) %>% list_rbind()
-message('2024')
-
-
-#CI_all <- rbind(Data_2023, Data_2024)
-CI_all <-  Data_2024
-message('CI_all created')
-#CI_all <-  Data_2023
-
-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')
-
diff --git a/r_app/ci_extraction.R b/r_app/ci_extraction.R
new file mode 100644
index 0000000..0522ce7
--- /dev/null
+++ b/r_app/ci_extraction.R
@@ -0,0 +1,241 @@
+library(here)
+library(sf)
+library(terra)
+library(tidyverse)
+library(lubridate)
+library(exactextractr)
+library(readxl)
+
+
+# Vang alle command line argumenten op
+args <- commandArgs(trailingOnly = TRUE)
+
+# Controleer of er ten minste één argument is doorgegeven
+if (length(args) == 0) {
+  stop("Geen argumenten doorgegeven aan het script")
+}
+
+# Converteer het eerste argument naar een numerieke waarde
+end_date <- as.Date(args[1])
+
+
+offset <- as.numeric(args[2])
+# Controleer of weeks_ago een geldig getal is
+if (is.na(offset)) {
+  # stop("Het argument is geen geldig getal")
+  offset <- 7
+}
+
+# Converteer het tweede argument naar een string waarde
+project_dir <- as.character(args[3])
+
+# Controleer of data_dir een geldige waarde is
+if (!is.character(project_dir)) {
+  project_dir <- "chemba"
+}
+
+
+laravel_storage_dir <- here("laravel_app/storage/app", project_dir)
+#preparing directories
+planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
+merged_final <- here(laravel_storage_dir, "merged_final_tif")
+
+new_project_question = FALSE
+planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
+merged_final <- here(laravel_storage_dir, "merged_final_tif")
+
+data_dir <- here(laravel_storage_dir, "Data")
+extracted_CI_dir <- here(data_dir, "extracted_ci")
+daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals")
+cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals")
+
+weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
+daily_vrt <- here(data_dir, "vrt")
+harvest_dir <-  here(data_dir, "HarvestData")
+
+source("parameters_project.R")
+source("utils_2.R")
+
+dir.create(here(laravel_storage_dir))
+dir.create(here(data_dir))
+dir.create(here(extracted_CI_dir))
+dir.create(here(daily_CI_vals_dir))
+dir.create(here(cumulative_CI_vals_dir))
+dir.create(here(weekly_CI_mosaic))
+dir.create(here(daily_vrt))
+dir.create(merged_final)
+dir.create(harvest_dir)
+
+# end_date <- lubridate::dmy("20-6-2024")
+week <- week(end_date)
+
+
+#weeks_ago = 0
+# Creating weekly mosaic
+#dates <- date_list(weeks_ago)
+dates <- date_list(end_date, offset)
+print(dates)
+#load pivot geojson
+# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
+
+# raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
+
+# filtered_files  <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
+#   compact() %>%
+#   flatten_chr()
+# head(filtered_files)
+
+
+
+# Extracting CI
+
+
+# pivot_sf_q <- st_read(here("..", "Data", "pivot_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
+# pivot_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant)   %>% group_by(pivot) %>% summarise() %>% vect()
+# message("pivot loaded")
+
+raster_files_NEW <- list.files(merged_final,full.names = T, pattern = ".tif")
+
+# pivots_dates0 <- readRDS(here(harvest_dir, "harvest_data_new")) %>% filter(
+#   pivot %in% c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13",
+#                "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3",
+#                "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3")
+# )
+
+# harvesting_data <- pivots_dates0 %>%
+#   select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023", "season_start_2024", "season_end_2024")) %>%
+#   pivot_longer(cols = starts_with("season"), names_to = "Year", values_to = "value") %>%
+#   separate(Year, into = c("name", "Year"), sep = "(?<=season_start|season_end)\\_", remove = FALSE) %>%
+#   mutate(name = str_to_title(name)) %>%
+#   pivot_wider(names_from = name, values_from = value) %>%
+#   rename(field = pivot_quadrant)
+
+#If run for the firsttime, it will extract all data since the start and put into a table.rds. otherwise it will only add on to the existing table.
+
+
+# Define the path to the file
+file_path <- here(cumulative_CI_vals_dir,"combined_CI_data.rds")
+
+# Check if the file exists
+if (!file.exists(file_path)) {
+  # Create the file with columns "column1" and "column2"
+  data <- data.frame(sub_field=NA, field=NA)
+  saveRDS(data, file_path)
+}
+
+print("combined_CI_data.rds exists, adding the latest image data to the table.")
+
+filtered_files  <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>%
+  compact() %>%
+  flatten_chr()
+
+walk(filtered_files, extract_rasters_daily, field_geojson= field_boundaries, quadrants = TRUE, daily_CI_vals_dir)
+
+extracted_values <-  list.files(daily_CI_vals_dir, full.names = TRUE)
+
+extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>%
+  compact() %>%
+  flatten_chr()
+
+pivot_stats <- extracted_values %>%
+  map(readRDS) %>% list_rbind() %>%
+  group_by(sub_field) %>%
+  summarise(across(everything(), ~ first(na.omit(.))))
+
+combined_CI_data <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #%>% drop_na(pivot_quadrant)
+pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data)
+# pivot_stats2 <- combined_CI_data
+print("All CI values extracted from latest 7 images.")
+saveRDS(combined_CI_data, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file
+
+
+# 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()
+
+
+# #2021
+# pivot_select_model_Data_2021 <- harvesting_data %>% filter(Year == 2021) %>% pull(field)
+# 
+# pivot_select_model_Data_2022 <- harvesting_data %>% filter(Year == 2022) %>% pull(field)
+
+# pivot_select_model_Data_2023 <- harvesting_data %>% filter(year == 2023) %>% filter(!is.na(season_start))  %>% pull(sub_field)
+
+pivot_select_model_Data_2024 <- harvesting_data %>% filter(year == 2024)%>% filter(!is.na(season_start)) %>% pull(sub_field)
+
+# pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022))
+# pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant )
+# 
+# pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023))
+# pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant)
+# 
+# pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024))
+# pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant)
+
+extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
+  #  field_names =  "AG1D06P"
+  #  field_names =  "1.13A"
+  #  harvesting_data = harvesting_data
+  # field_CI_data = pivot_stats_long
+  #  season= 2023
+  
+  filtered_harvesting_data <- harvesting_data %>%
+    na.omit() %>%
+    filter(year == season, sub_field %in% field_names)
+  
+  filtered_field_CI_data <- field_CI_data %>%
+    filter(sub_field %in% field_names)
+  
+  # CI <- map_df(field_names, ~ {
+  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)
+  
+  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) %>%
+    mutate(DOY = seq(1, n(), 1),
+           model = paste0("Data", season, " : ", field_names),
+           season = season,
+           sub_field = field_names)
+  # }) #%>%
+  #{if (length(field_names) > 0) message("Done!")}
+  
+  return(CI)
+}
+
+## Extracting the correct CI values
+# Data_2021 <- map(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021)) %>% list_rbind()
+# message('2021')
+# Data_2022 <- map(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022)) %>% list_rbind()
+# message('2022')
+# Data_2023 <- map(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023)) %>% list_rbind()
+# message('2023')
+Data_2024 <- map(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024)) %>% list_rbind()
+message('2024')
+
+
+#CI_all <- rbind(Data_2023, Data_2024)
+CI_all <-  Data_2024
+message('CI_all created')
+#CI_all <-  Data_2023
+
+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')
+
diff --git a/r_app/ci_extraction_utils.R b/r_app/ci_extraction_utils.R
new file mode 100644
index 0000000..f092da4
--- /dev/null
+++ b/r_app/ci_extraction_utils.R
@@ -0,0 +1,284 @@
+# Utils for ci extraction
+
+date_list <- function(end_date, offset){
+  offset <- as.numeric(offset) - 1
+  end_date <- as.Date(end_date)
+  start_date <- end_date - lubridate::days(offset)
+  
+  week <- week(start_date)
+  year <- year(start_date)
+  days_filter <- seq(from = start_date, to = end_date, by = "day")
+  
+  return(list("week" = week, "year" = year, "days_filter" = days_filter))
+}
+
+# date_list <- function(weeks_in_the_paste){
+#   week <- week(Sys.Date()- weeks(weeks_in_the_paste) )
+#   year <- year(Sys.Date()- weeks(weeks_in_the_paste) )
+#   days_filter <- Sys.Date() - weeks(weeks_in_the_paste) - days(0:6)
+#
+#   return(c("week" = week,
+#            "year" = year,
+#            "days_filter" = list(days_filter)))
+#
+# }
+
+CI_func <- function(x, drop_layers = FALSE){
+  CI <- x[[4]]/x[[2]]-1
+  add(x) <- CI
+  names(x) <- c("red", "green", "blue","nir", "cloud" ,"CI")
+  if(drop_layers == FALSE){
+    return(x)
+  }else{
+    return(x$CI)
+  }
+}
+
+mask_raster <- function(raster, fields){
+  # x <- rast(filtered_files[1])
+  x <- rast(raster) 
+  emtpy_or_full <- global(x, sum)
+  
+  if(emtpy_or_full[1,] >=  2000000){
+    names(x) <- c("red", "green", "blue","nir", "cloud")
+    cloud <- x$cloud
+    
+    cloud[cloud == 0 ] <- NA
+    x_masked <- mask(x, cloud, inverse = T) %>% crop(.,fields, mask = TRUE )
+    x_masked <-  x_masked %>% CI_func()
+    
+    message(raster, " masked")
+    return(x_masked)
+  }
+}
+
+date_extract <- function(file_path) {
+  str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}")
+}
+
+extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) {
+  # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE)
+  # date <- date_extract(filtered_files[1])
+  # field_geojson <- sf::st_as_sf(pivot_sf_q)
+  field_geojson <- sf::st_as_sf(field_geojson)
+  x <- rast(file[1])
+  date <- date_extract(file)
+  pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>%
+    st_drop_geometry() %>% rename("{date}" := mean_CI)
+  save_suffix <- if (quadrants){"quadrant"} else {"whole_field"}
+  save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
+  saveRDS(pivot_stats, save_path)
+}
+
+right = function(text, num_char) {
+  substr(text, nchar(text) - (num_char-1), nchar(text))
+}
+
+extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
+  #  field_names =  "1.2A"
+  #  harvesting_data = harvesting_data
+  # field_CI_data = pivot_stats_long
+  #  season= 2021
+  
+  filtered_harvesting_data <- harvesting_data %>%
+    filter(year == season, field %in% field_names)
+  
+  filtered_field_CI_data <- field_CI_data %>%
+    filter(field %in% field_names)
+  
+  # CI <- map_df(field_names, ~ {
+  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)
+  
+  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) %>%
+    mutate(DOY = seq(1, n(), 1),
+           model = paste0("Data", season, " : ", field_names),
+           season = season,
+           field = field_names)
+  # }) #%>%
+  #{if (length(field_names) > 0) message("Done!")}
+  
+  return(CI)
+}
+# 
+# load_fields <- function(geojson_path) {
+#   field_geojson <- st_read(geojson_path) %>%
+#     select(pivot, pivot_quadrant) %>%
+#     vect()
+#   return(field_geojson)
+# }
+# 
+# load_harvest_data <- function(havest_data_path){
+#   harvest_data <- readRDS(havest_data_path)
+#   return(harvest_data)
+# }
+# 
+# load_rasters <- function(raster_path, dates) {
+#   raster_files <- list.files(raster_path, full.names = TRUE, pattern = ".tif")
+#   filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
+#     compact() %>%
+#     flatten_chr()
+#   
+#   return(filtered_files)
+# }
+# 
+# mask_and_set_names <- function(filtered_files, fields) {
+#   rasters_masked <- map(filtered_files, mask_raster, fields = fields) %>% set_names(filtered_files)
+#   rasters_masked[sapply(rasters_masked, is.null)] <- NULL
+#   rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract))
+#   
+#   return(rasters_masked)
+# }
+# 
+# calculate_total_pix_area <- function(filtered_files, fields_geojson) {
+#   # total_pix_area <- rast(filtered_files[1]) %>%
+#   #   subset(1) %>%
+#   #   crop(fields_geojson, mask = TRUE)%>%
+#   #   global(.data, fun = "notNA")
+#   total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(fields_geojson, mask = TRUE) %>% freq(., usenames = TRUE)
+#   
+#   return(total_pix_area)
+# }
+# 
+# cloud_layer_extract <- function(rasters_masked){
+#   cloud_layer_rast <- map(rasters_masked, function(spatraster) {
+#     spatraster[[5]]
+#   }) %>% rast()
+#   
+#   return(cloud_layer_rast)
+# }
+# 
+# calculate_cloud_coverage <- function(cloud_layer_rast, total_pix_area) {
+#   cloud_perc_list <- freq(cloud_layer_rast, usenames = TRUE) %>%
+#     mutate(cloud_perc = (100 -((count/total_pix_area$count)*100)),
+#            cloud_thres_5perc = as.integer(cloud_perc < 5),
+#            cloud_thres_40perc = as.integer(cloud_perc < 40)) %>%
+#     rename(Date = layer) %>% select(-value, -count)
+#   
+#   cloud_index_5perc <- which(cloud_perc_list$cloud_thres_5perc == max(cloud_perc_list$cloud_thres_5perc))
+#   cloud_index_40perc <- which(cloud_perc_list$cloud_thres_40perc == max(cloud_perc_list$cloud_thres_40perc))
+#   
+#   return(list(cloud_perc_list = cloud_perc_list, cloud_index_5perc = cloud_index_5perc, cloud_index_40perc = cloud_index_40perc))
+# }
+# 
+# process_cloud_coverage <- function(cloud_coverage, rasters_masked) {
+#   if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) > 1) {
+#     message("More than 1 raster without clouds (<5%), max mosaic made ")
+#     
+#     cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_5perc]
+#     rsrc <- sprc(cloudy_rasters_list)
+#     x <- mosaic(rsrc)
+#     names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+#     
+#   } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) == 1) {
+#     message("Only 1 raster without clouds (<5%)")
+#     
+#     x <- rast(rasters_masked[cloud_coverage$cloud_index_5perc[1]])
+#     names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+#     
+#   } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) > 1) {
+#     message("More than 1 image contains clouds, composite made of <40% cloud cover images")
+#     
+#     cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_40perc]
+#     rsrc <- sprc(cloudy_rasters_list)
+#     x <- mosaic(rsrc)
+#     names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+#     
+#   } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) == 0) {
+#     message("No cloud free images available")
+#     x <-  rast(rasters_masked[1])
+#     x[x] <- NA
+#     names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+#   }
+#   
+#   return(x)
+# }
+
+# extract_rasters_daily_func <- function(daily_vals_dir, filtered_files, fields_geojson) {
+#   extracted_files <- walk(filtered_files, extract_rasters_daily, field_geojson = fields_geojson, quadrants = TRUE, daily_vals_dir)
+#   return(extracted_files)
+# }
+
+# CI_load <- function(daily_vals_dir, grouping_variable){
+#   extracted_values <-  list.files(here(daily_vals_dir), full.names = TRUE)
+#   
+#   field_CI_values <- extracted_values %>%
+#     map_dfr(readRDS) %>%
+#     group_by(.data[[grouping_variable]]) %>%
+#     summarise(across(everything(), ~ first(na.omit(.))))
+#   return(field_CI_values)
+# }
+# 
+# CI_long <- function(field_CI_values, pivot_long_cols){
+#   field_CI_long <- field_CI_values %>%
+#     gather("Date", value, -pivot_long_cols) %>%
+#     mutate(Date = right(Date, 8),
+#            Date = ymd(Date)
+#     ) %>%
+#     drop_na(c("value","Date")) %>%
+#     mutate(value = as.numeric(value))%>%
+#     filter_all(all_vars(!is.infinite(.)))%>%
+#     rename(field = pivot_quadrant)
+#   
+#   return(field_CI_long)
+# }
+# 
+# process_year_data <- function(year, harvest_data, field_CI_long) {
+#   pivots_dates_year <- harvest_data %>% na.omit() %>% filter(year == year)
+#   pivot_select_model_year <- unique(pivots_dates_year$field)
+#   
+#   data <- map_dfr(pivot_select_model_year, ~ extract_CI_data(.x, harvest_data, field_CI_long, season = year))
+#   
+#   return(data)
+# }
+
+
+
+#functions for CI_data_prep
+# create_mask_and_crop <- function(file, pivot_sf_q) {
+#   # file <- filtered_files[5]
+#   message("starting ", file)
+#   loaded_raster <- rast(file)
+#   names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
+#   #   names(CI) <- c("green","nir")
+#   message("raster loaded")
+#   
+#   # CI <- CI[[2]] / CI[[1]] - 1
+#   CI <- loaded_raster$NIR / loaded_raster$Green - 1
+#   
+#   loaded_raster$CI <- CI
+#   #   CI <- CI$nir/CI$green-1
+#   message("CI calculated")
+#   loaded_raster <- terra::crop(loaded_raster, pivot_sf_q, mask = TRUE) #%>% CI_func()
+#   
+#   loaded_raster[loaded_raster == 0] <- NA
+#   # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI")
+#   # v_crop$CI <- v_crop$CI - 1
+#   new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
+#   writeRaster(loaded_raster, new_file, overwrite = TRUE)
+#   
+#   vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
+#   terra::vrt(new_file, vrt_file, overwrite = TRUE)
+#   
+#   # v_crop <- mask_raster(v, pivot_sf_q)
+#   return(loaded_raster)
+# }
+
+
+# extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) {
+#   # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE)
+#   # date <- date_extract(filtered_files[1])
+#   # field_geojson <- sf::st_as_sf(pivot_sf_q)
+#   field_geojson <- sf::st_as_sf(field_geojson)
+#   x <- rast(file[1])
+#   date <- date_extract(file)
+#   pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>%
+#     st_drop_geometry() %>% rename("{date}" := mean_CI)
+#   save_suffix <- if (quadrants){"quadrant"} else {"whole_field"}
+#   save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
+#   saveRDS(pivot_stats, save_path)
+# }
+
diff --git a/r_app/mosaic_creation.R b/r_app/mosaic_creation.R
new file mode 100644
index 0000000..055604b
--- /dev/null
+++ b/r_app/mosaic_creation.R
@@ -0,0 +1,216 @@
+# activeer de renv omgeving;
+
+# renv::activate('~/smartCane/r_app')
+# renv::restore()
+
+
+
+library(here)
+library(sf)
+library(terra)
+library(tidyverse)
+library(lubridate)
+# library(exactextractr)
+library(readxl)
+#funcion CI_prep package
+
+
+# Vang alle command line argumenten op
+args <- commandArgs(trailingOnly = TRUE)
+
+# Controleer of er ten minste één argument is doorgegeven
+if (length(args) == 0) {
+  stop("Geen argumenten doorgegeven aan het script")
+}
+
+# Converteer het eerste argument naar een numerieke waarde
+end_date <- as.Date(args[1])
+
+
+offset <- as.numeric(args[2])
+# Controleer of weeks_ago een geldig getal is
+if (is.na(offset)) {
+  # stop("Het argument is geen geldig getal")
+  offset <- 7
+}
+
+# Converteer het tweede argument naar een string waarde
+project_dir <- as.character(args[3])
+
+# Controleer of data_dir een geldige waarde is
+if (!is.character(project_dir)) {
+  project_dir <- "chemba"
+}
+
+
+laravel_storage_dir <- here("laravel_app/storage/app", project_dir)
+#preparing directories
+planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
+merged_final <- here(laravel_storage_dir, "merged_final_tif")
+
+new_project_question = FALSE
+planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
+merged_final <- here(laravel_storage_dir, "merged_final_tif")
+
+data_dir <- here(laravel_storage_dir, "Data")
+extracted_CI_dir <- here(data_dir, "extracted_ci")
+daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals")
+cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals")
+
+weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
+daily_vrt <- here(data_dir, "vrt")
+harvest_dir <-  here(data_dir, "HarvestData")
+
+source("parameters_project.R")
+source("utils_1.R")
+
+dir.create(here(laravel_storage_dir))
+dir.create(here(data_dir))
+dir.create(here(extracted_CI_dir))
+dir.create(here(daily_CI_vals_dir))
+dir.create(here(cumulative_CI_vals_dir))
+dir.create(here(weekly_CI_mosaic))
+dir.create(here(daily_vrt))
+dir.create(merged_final)
+dir.create(harvest_dir)
+
+# end_date <- lubridate::dmy("20-6-2024")
+week <- week(end_date)
+
+
+#weeks_ago = 0
+
+# Creating weekly mosaic
+#dates <- date_list(weeks_ago)
+dates <- date_list(end_date, offset)
+print(dates)
+#load pivot geojson
+# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
+
+raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
+
+filtered_files  <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
+  compact() %>%
+  flatten_chr()
+head(filtered_files)
+
+# filtered_files <- raster_files #for first CI extraction
+
+# create_mask_and_crop <- function(file, field_boundaries) {
+# message("starting ", file)
+#   CI <- rast(file)
+# #   names(CI) <- c("green","nir")
+# message("raster loaded")
+# 
+#   CI <- CI[[2]]/CI[[1]]-1
+# #   CI <- CI$nir/CI$green-1
+# message("CI calculated")
+#   CI <- terra::crop(CI, field_boundaries, mask = TRUE) #%>% CI_func()
+# 
+#   new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
+#   writeRaster(CI, new_file, overwrite = TRUE)
+# 
+#   vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
+#   terra::vrt(new_file, vrt_file, overwrite = TRUE)
+# 
+#   return(CI)
+# }
+
+
+vrt_list <- list()
+
+# for (file in raster_files) {
+#   v_crop <- create_mask_and_crop(file, field_boundaries)
+#   message(file, " processed")
+#   gc()
+# }
+
+for (file in filtered_files) {
+  v_crop <- create_mask_and_crop(file, field_boundaries)
+  emtpy_or_full <- global(v_crop, "notNA")
+  
+  vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
+  
+  if(emtpy_or_full[1,] >  10000){
+    vrt_list[vrt_file] <- vrt_file
+    
+  }else{
+    file.remove(vrt_file)
+    
+  }
+  
+  message(file, " processed")
+  gc()
+}
+
+# list_global <- list_global %>% flatten_chr()
+vrt_list <- vrt_list %>% flatten_chr()
+
+total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>%
+  crop(field_boundaries, mask = TRUE) %>%
+  global(., fun="notNA") #%>%
+
+layer_5_list <- purrr::map(vrt_list, function(vrt_list) {
+  rast(vrt_list[1]) %>% terra::subset(1)
+}) %>% rast()
+
+missing_pixels_count <- layer_5_list %>% global(., fun="notNA") %>%
+  mutate(
+    total_pixels = total_pix_area$notNA,
+    missing_pixels_percentage = round(100 -((notNA/total_pix_area$notNA)*100)),
+    thres_5perc = as.integer(missing_pixels_percentage < 5),
+    thres_40perc = as.integer(missing_pixels_percentage < 45)
+  )
+
+index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc)  )
+index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc))
+
+## Create mosaic
+
+if(sum(missing_pixels_count$thres_5perc)>1){
+  message("More than 1 raster without clouds (<5%), max composite made")
+  
+  cloudy_rasters_list <- vrt_list[index_5perc]
+  
+  rsrc <- sprc(cloudy_rasters_list)
+  x <- mosaic(rsrc, fun = "max")
+  
+  # names(x) <-  "CI"
+  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
+}else if(sum(missing_pixels_count$thres_5perc)==1){
+  message("Only 1 raster without clouds (<5%)")
+  
+  x <- rast(vrt_list[index_5perc[1]])
+  # names(x) <-  c("CI")
+  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
+}else if(sum(missing_pixels_count$thres_40perc)>1){
+  message("More than 1 image contains clouds, composite made of <40% cloud cover images")
+  
+  cloudy_rasters_list <- vrt_list[index_40perc]
+  
+  rsrc <- sprc(cloudy_rasters_list)
+  x <- mosaic(rsrc, fun = "max")
+  # names(x) <-  "CI"
+  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
+}else if(sum(missing_pixels_count$thres_40perc)==1){
+  message("Only 1 image available but contains clouds")
+  
+  x <- rast(vrt_list[index_40perc[1]])
+  # names(x) <-  c("CI")
+  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
+}else{
+  message("No cloud free images available, all images combined")
+  
+  rsrc <- sprc(vrt_list)
+  x <- mosaic(rsrc, fun = "max")
+  # x <-  rast(vrt_list[1]) %>% setValues(NA)
+  # names(x) <-  c("CI")
+  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
+}
+
+plot(x$CI, main = paste("CI map ", dates$week))
+plotRGB(x, main = paste("RGB map ", dates$week))
+
+file_path_tif <- here(weekly_CI_mosaic ,paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif"))
+writeRaster(x, file_path_tif, overwrite=TRUE)
+message("Raster written/made at: ", file_path_tif)
diff --git a/r_app/mosaic_creation_utils.R b/r_app/mosaic_creation_utils.R
new file mode 100644
index 0000000..b07176f
--- /dev/null
+++ b/r_app/mosaic_creation_utils.R
@@ -0,0 +1,217 @@
+# utils for mosaic creation
+
+date_list <- function(end_date, offset){
+  offset <- as.numeric(offset) - 1
+  end_date <- as.Date(end_date)
+  start_date <- end_date - lubridate::days(offset)
+  
+  week <- week(start_date)
+  year <- year(start_date)
+  days_filter <- seq(from = start_date, to = end_date, by = "day")
+  
+  return(list("week" = week, "year" = year, "days_filter" = days_filter))
+}
+
+# date_list <- function(weeks_in_the_paste){
+#   week <- week(Sys.Date()- weeks(weeks_in_the_paste) )
+#   year <- year(Sys.Date()- weeks(weeks_in_the_paste) )
+#   days_filter <- Sys.Date() - weeks(weeks_in_the_paste) - days(0:6)
+#
+#   return(c("week" = week,
+#            "year" = year,
+#            "days_filter" = list(days_filter)))
+#
+# }
+create_mask_and_crop <- function(file, field_boundaries) {
+  # file <- filtered_files[5]
+  message("starting ", file)
+  loaded_raster <- rast(file)
+  names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
+  #   names(CI) <- c("green","nir")
+  message("raster loaded")
+  
+  # CI <- CI[[2]] / CI[[1]] - 1
+  CI <- loaded_raster$NIR / loaded_raster$Green - 1
+  
+  loaded_raster$CI <- CI
+  #   CI <- CI$nir/CI$green-1
+  message("CI calculated")
+  loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE) #%>% CI_func()
+  
+  loaded_raster[loaded_raster == 0] <- NA
+  # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI")
+  # v_crop$CI <- v_crop$CI - 1
+  new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
+  writeRaster(loaded_raster, new_file, overwrite = TRUE)
+  
+  vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
+  terra::vrt(new_file, vrt_file, overwrite = TRUE)
+  
+  # v_crop <- mask_raster(v, pivot_sf_q)
+  return(loaded_raster)
+}
+
+CI_func <- function(x, drop_layers = FALSE){
+  CI <- x[[4]]/x[[2]]-1
+  add(x) <- CI
+  names(x) <- c("red", "green", "blue","nir", "cloud" ,"CI")
+  if(drop_layers == FALSE){
+    return(x)
+  }else{
+    return(x$CI)
+  }
+}
+
+mask_raster <- function(raster, fields){
+  # x <- rast(filtered_files[1])
+  x <- rast(raster) 
+  emtpy_or_full <- global(x, sum)
+  
+  if(emtpy_or_full[1,] >=  2000000){
+    names(x) <- c("red", "green", "blue","nir", "cloud")
+    cloud <- x$cloud
+    
+    cloud[cloud == 0 ] <- NA
+    x_masked <- mask(x, cloud, inverse = T) %>% crop(.,fields, mask = TRUE )
+    x_masked <-  x_masked %>% CI_func()
+    
+    message(raster, " masked")
+    return(x_masked)
+  }
+}
+
+date_extract <- function(file_path) {
+  str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}")
+}
+
+
+
+right = function(text, num_char) {
+  substr(text, nchar(text) - (num_char-1), nchar(text))
+}
+
+# load_fields <- function(geojson_path) {
+#   field_geojson <- st_read(geojson_path) %>%
+#     select(pivot, pivot_quadrant) %>%
+#     vect()
+#   return(field_geojson)
+# }
+# 
+# load_harvest_data <- function(havest_data_path){
+#   harvest_data <- readRDS(havest_data_path)
+#   return(harvest_data)
+# }
+
+# load_rasters <- function(raster_path, dates) {
+#   raster_files <- list.files(raster_path, full.names = TRUE, pattern = ".tif")
+#   filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
+#     compact() %>%
+#     flatten_chr()
+#   
+#   return(filtered_files)
+# }
+
+# mask_and_set_names <- function(filtered_files, fields) {
+#   rasters_masked <- map(filtered_files, mask_raster, fields = fields) %>% set_names(filtered_files)
+#   rasters_masked[sapply(rasters_masked, is.null)] <- NULL
+#   rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract))
+#   
+#   return(rasters_masked)
+# }
+
+# calculate_total_pix_area <- function(filtered_files, fields_geojson) {
+#   # total_pix_area <- rast(filtered_files[1]) %>%
+#   #   subset(1) %>%
+#   #   crop(fields_geojson, mask = TRUE)%>%
+#   #   global(.data, fun = "notNA")
+#   total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(fields_geojson, mask = TRUE) %>% freq(., usenames = TRUE)
+#   
+#   return(total_pix_area)
+# }
+
+# cloud_layer_extract <- function(rasters_masked){
+#   cloud_layer_rast <- map(rasters_masked, function(spatraster) {
+#     spatraster[[5]]
+#   }) %>% rast()
+#   
+#   return(cloud_layer_rast)
+# }
+
+# calculate_cloud_coverage <- function(cloud_layer_rast, total_pix_area) {
+#   cloud_perc_list <- freq(cloud_layer_rast, usenames = TRUE) %>%
+#     mutate(cloud_perc = (100 -((count/total_pix_area$count)*100)),
+#            cloud_thres_5perc = as.integer(cloud_perc < 5),
+#            cloud_thres_40perc = as.integer(cloud_perc < 40)) %>%
+#     rename(Date = layer) %>% select(-value, -count)
+#   
+#   cloud_index_5perc <- which(cloud_perc_list$cloud_thres_5perc == max(cloud_perc_list$cloud_thres_5perc))
+#   cloud_index_40perc <- which(cloud_perc_list$cloud_thres_40perc == max(cloud_perc_list$cloud_thres_40perc))
+#   
+#   return(list(cloud_perc_list = cloud_perc_list, cloud_index_5perc = cloud_index_5perc, cloud_index_40perc = cloud_index_40perc))
+# }
+
+process_cloud_coverage <- function(cloud_coverage, rasters_masked) {
+  if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) > 1) {
+    message("More than 1 raster without clouds (<5%), max mosaic made ")
+    
+    cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_5perc]
+    rsrc <- sprc(cloudy_rasters_list)
+    x <- mosaic(rsrc)
+    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+    
+  } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) == 1) {
+    message("Only 1 raster without clouds (<5%)")
+    
+    x <- rast(rasters_masked[cloud_coverage$cloud_index_5perc[1]])
+    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+    
+  } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) > 1) {
+    message("More than 1 image contains clouds, composite made of <40% cloud cover images")
+    
+    cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_40perc]
+    rsrc <- sprc(cloudy_rasters_list)
+    x <- mosaic(rsrc)
+    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+    
+  } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) == 0) {
+    message("No cloud free images available")
+    x <-  rast(rasters_masked[1])
+    x[x] <- NA
+    names(x) <-  c("red", "green", "blue", "nir", "cloud", "CI")
+  }
+  
+  return(x)
+}
+
+
+#functions for CI_data_prep
+# create_mask_and_crop <- function(file, pivot_sf_q) {
+#   # file <- filtered_files[5]
+#   message("starting ", file)
+#   loaded_raster <- rast(file)
+#   names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
+#   #   names(CI) <- c("green","nir")
+#   message("raster loaded")
+#   
+#   # CI <- CI[[2]] / CI[[1]] - 1
+#   CI <- loaded_raster$NIR / loaded_raster$Green - 1
+#   
+#   loaded_raster$CI <- CI
+#   #   CI <- CI$nir/CI$green-1
+#   message("CI calculated")
+#   loaded_raster <- terra::crop(loaded_raster, pivot_sf_q, mask = TRUE) #%>% CI_func()
+#   
+#   loaded_raster[loaded_raster == 0] <- NA
+#   # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI")
+#   # v_crop$CI <- v_crop$CI - 1
+#   new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
+#   writeRaster(loaded_raster, new_file, overwrite = TRUE)
+#   
+#   vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
+#   terra::vrt(new_file, vrt_file, overwrite = TRUE)
+#   
+#   # v_crop <- mask_raster(v, pivot_sf_q)
+#   return(loaded_raster)
+# }
+
+
diff --git a/r_app/utils_3.R b/r_app/utils_3.R
new file mode 100644
index 0000000..27a4e8e
--- /dev/null
+++ b/r_app/utils_3.R
@@ -0,0 +1,161 @@
+# utils for report 
+#functions for rmarkdown file
+
+
+create_RGB_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = TRUE, legend_is_portrait = FALSE, week, age, red = TRUE) {
+  r <- if (red) 1 else 4  # Set r based on the value of red
+  title <- if (red) paste0("RGB image of the fields") else paste0("False colour image of the fields")
+  
+  tm_shape(pivot_raster, unit = "m") + tm_rgb(r = r, g = 2, b = 3, max.value = 255) +
+    tm_layout(main.title = title, 
+              main.title.size = 1) +
+    tm_scale_bar(position = c("right", "top"), text.color = "black") +
+    tm_compass(position = c("right", "top"), text.color = "black") +
+    tm_shape(pivot_shape) + tm_borders(col = "gray") +
+    tm_text("sub_field", size = 0.6, col = "gray") +
+    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha = 0.5)
+}
+
+create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, legend_only = F){
+  tm_shape(pivot_raster, unit = "m")+
+    tm_raster(breaks =  CI_breaks, palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) +
+    tm_layout(main.title = paste0("Max CI week ", week,"\n",  age, " weeks old"),
+              main.title.size = 1, legend.show = show_legend, legend.only = legend_only) +
+    tm_shape(pivot_shape) +
+    tm_borders(lwd = 3) +  tm_text("sub_field", size = 1/2) +
+    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5) +tmap_options(check.and.fix = TRUE)
+}
+
+create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age){
+  tm_shape(pivot_raster, unit = "m")+
+    tm_raster(breaks = CI_diff_breaks, palette = "PRGn",legend.is.portrait = legend_is_portrait ,midpoint = 0, title = "CI difference") +
+    tm_layout(main.title =  paste0("CI change week ", week_1, "- week ",week_2, "\n", age," weeks old"),
+              main.title.size = 1, legend.show = show_legend) +
+    tm_shape(pivot_shape) +
+    tm_borders(lwd = 3) +  tm_text("sub_field", size = 1/2) +
+    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
+}
+
+ci_plot <- function(pivotName){
+  # pivotName = "MV2B09"
+  # pivotName = "1.1"
+  
+  pivotShape <- AllPivots_merged %>% terra::subset(field %in% pivotName)  %>% st_transform(crs(CI))
+  # age <- AllPivots %>% dplyr::filter(field %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique() %>%
+  #   mutate(Age = round(Age))
+  
+  age <- AllPivots %>%
+    group_by(field) %>%
+    filter(Season == max(Season, na.rm = TRUE), field %in% pivotName) %>%
+    dplyr::select(Age)%>% st_drop_geometry()  %>% unique() 
+  
+  AllPivots2 <- AllPivots0 %>% dplyr::filter(field %in% pivotName)
+  
+  singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape)
+  singlePivot_m1 <- CI_m1 %>% crop(., pivotShape) %>% mask(., pivotShape)
+  singlePivot_m2 <- CI_m2 %>% crop(., pivotShape) %>% mask(., pivotShape)
+  # singlePivot_m3 <- CI_m3 %>% crop(., pivotShape) %>% mask(., pivotShape)
+  
+  singlePivot_RGB <- RGB_raster %>% crop(., pivotShape) %>% mask(., pivotShape)
+  singlePivot_false <- RGB_raster_stretch %>% crop(., pivotShape) %>% mask(., pivotShape)
+  
+  abs_CI_last_week <- last_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape)
+  abs_CI_three_week <- three_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape)
+  
+  # planting_date <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% ungroup()  %>% dplyr::select(planting_date)  %>% unique()
+  
+  joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape))  %>%  dplyr::filter(field %in% pivotName) %>% st_crop(., pivotShape)
+  
+  # CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait = T, week = week_minus_2, age = age  -2)
+  Legend_map <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age  -1, legend_only = T)
+  CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age  -1)
+  CImap <- create_CI_map(singlePivot, AllPivots2,  joined_spans2, show_legend= F, legend_is_portrait = T, week = week, age = age )
+  RGBmap <- create_RGB_map(singlePivot_false, AllPivots2,  joined_spans2, show_legend= F,  week = week, age = age, red =T )
+  Falsemap <- create_RGB_map(singlePivot_false, AllPivots2,  joined_spans2, show_legend= F,  week = week, age = age, red =F )
+  
+  
+  CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week,AllPivots2,  joined_spans2,  show_legend = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_1, age = age)
+  CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2,  joined_spans2, show_legend = F, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age)
+  
+  # tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1)
+  tst <- tmap_arrange(RGBmap,Falsemap,
+                      CImap_m1, CImap, 
+                      CI_max_abs_last_week, CI_max_abs_three_week, 
+                      ncol = 2)
+  
+  cat(paste("## field", pivotName, "-", age$Age[1], "weeks after planting/harvest", "\n"))
+  # cat("\n")
+  # cat('<h2> Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest <h2>')
+  # cat(paste("# Pivot",pivots$pivot[i],"\n"))
+  print(tst)
+  
+}
+
+
+subchunkify <- function(g, fig_height=7, fig_width=5) {
+  g_deparsed <- paste0(deparse(
+    function() {g}
+  ), collapse = '')
+  
+  sub_chunk <- paste0("```{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", echo=FALSE}",
+                      "\n(",
+                      g_deparsed
+                      , ")()",
+                      "\n```
+  ")
+  
+  cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE))
+}
+
+cum_ci_plot <- function(pivotName){
+  
+  # pivotName = "2.1"
+  
+  # Check if pivotName exists in the data
+  if (!pivotName %in% unique(CI_quadrant$field)) {
+    # message("PivotName '", pivotName, "' not found. Plotting empty graph.")
+    g <- ggplot() +   labs(title = "Empty Graph - Yield dates missing")
+    return(
+      subchunkify(g, fig_height=2, fig_width = 10)
+    )
+  } else {
+    # message("PivotName '", pivotName, "'  found. Plotting normal graph.")
+    data_ci <- CI_quadrant %>% filter(field %in% pivotName)
+    
+    
+    
+    data_ci2 <- data_ci  %>% mutate(CI_rate = cumulative_CI/DOY,
+                                    week = week(Date))%>% group_by(sub_field) %>%
+      mutate(mean_rolling10 = rollapplyr(CI_rate  , width = 10, FUN = mean, partial = TRUE))
+    
+    # date_preperation_perfect_pivot <- data_ci2 %>% group_by(season)  %>% summarise(min_date = min(Date),
+    #                                                                                max_date = max(Date),
+    #                                                                                days = max_date - min_date)
+    
+    # Identify unique seasons
+    filtered_data <- data_ci2 %>%
+      group_by(season) %>%
+      mutate(rank = dense_rank(desc(season))) %>%
+      filter(rank <= 2) %>%
+      ungroup() %>%
+      dplyr::select(-rank)
+    
+    
+    # g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) +
+    g <- ggplot(data= filtered_data ) +
+      # geom_line(aes(Date, mean_rolling10, col = sub_field)) +
+      geom_line(aes(Date, CI_rate, col = sub_field)) +
+      facet_wrap(~season, scales = "free_x") +
+      # geom_line(data= perfect_pivot, aes(Date , mean_rolling10, col = "Model CI (p5.1 Data 2022, \n date x axis is fictive)"), lty="11",size=1) +
+      labs(title = paste("CI rate - field", pivotName),
+           y = "CI rate (cumulative CI / Age)")+
+      # scale_y_continuous(limits=c(0.5,3), breaks = seq(0.5, 3, 0.5))+
+      scale_x_date(date_breaks = "1 month", date_labels =  "%m-%Y") +
+      theme(axis.text.x = element_text(angle = 60, hjust = 1),
+            legend.justification=c(1,0), legend.position = c(1, 0),
+            legend.title = element_text(size = 8),
+            legend.text = element_text(size = 8)) +
+      guides(color = guide_legend(nrow = 2, byrow = TRUE))
+    subchunkify(g, fig_height=6, fig_width = 10)
+  }
+}
\ No newline at end of file
diff --git a/update_RDS.sh b/update_RDS.sh
new file mode 100755
index 0000000..6d74327
--- /dev/null
+++ b/update_RDS.sh
@@ -0,0 +1,39 @@
+#!/bin/bash
+
+end_date=$(date +'%Y-%m-%d')
+offset=7
+project_dir="chemba"
+
+# Parse command line arguments
+for arg in "$@"; do
+    case $arg in
+        --end_date=*)
+            end_date="${arg#*=}"
+            ;;
+        --offset=*)
+        offset="${arg#*=}"
+        ;;
+        --project_dir=*)
+            project_dir="${arg#*=}"
+            ;;
+        *)
+            echo "Unknown option: $arg"
+            exit 1
+            ;;
+    esac
+    shift
+done
+
+echo "end_date: $end_date"
+echo "offset: $offset"
+
+# Check if required arguments are set
+if [ -z "$end_date" ] || [ -z "$project_dir" ] || [ -z "$offset" ]; then
+    echo "Missing arguments. Use: build_RDS.sh --end_date=2024-01-01 --offset=7 --project_dir=chemba"
+    exit 1
+fi
+
+echo ci_extraction.R $end_date $offset $project_dir
+
+cd ../r_app
+Rscript ci_extraction.R $end_date $offset $project_dir
\ No newline at end of file