parameterised harvestdata etc

r_app/2_CI_data_prep.R

@@ -11,6 +11,7 @@ library(lubridate)
 library(exactextractr)
 
 library(CIprep)
+source(here("r_app", "parameters_project.R"))
 
 # Vang alle command line argumenten op
 args <- commandArgs(trailingOnly = TRUE)
@@ -55,6 +56,8 @@ 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")
+
+dir.create(here(laravel_storage_dir))
 dir.create(here(data_dir))
 dir.create(here(extracted_CI_dir))
 dir.create(here(daily_CI_vals_dir))
@@ -70,7 +73,7 @@ dir.create(harvest_dir)
 dates <- date_list(weeks_ago)
 
 #load pivot geojson
-pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
+# 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")
 head(raster_files)
@@ -79,65 +82,72 @@ filtered_files  <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x =
   compact() %>%
   flatten_chr()
 
+# filtered_files <- raster_files #for first CI extraction
 
-#rasters_masked <- map(filtered_files, mask_raster, fields = pivot_sf_q) %>% set_names(filtered_files)
-
-
-# rasters_masked <- list()  # Creëer een lege lijst om de resultaten op te slaan
-
-# for (i in seq_along(filtered_files[1])) {
-#     file_name <- filtered_files[i]
-#     result <- mask_raster(file_name, fields = pivot_sf_q)
-#     rasters_masked[[file_name]] <- result
+# 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, pivot_sf_q) {
-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, pivot_sf_q, mask = TRUE) #%>% CI_func()
-
-  # v_crop[v_crop == 0] <- NA
+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(CI, new_file, overwrite = TRUE)
-
+  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(CI)
+  return(loaded_raster)
 }
 
-# rasters_masked <- map(filtered_files, create_mask_and_crop, pivot_sf_q)
-# list_global <- list()
 vrt_list <- list()
 
 for (file in filtered_files) {
-  v_crop <- create_mask_and_crop(file, pivot_sf_q)
+  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){
-    # list_global[file] <- file
     vrt_list[vrt_file] <- vrt_file
 
   }else{
     file.remove(vrt_file)
-    # file.remove(file)
-    # message(file, " removed")
+
   }
 
-  # Save the processed raster to a new file
-  # output_file <- here(data_dir, "vrt", paste0(basename(file), "_processed"))
-  # terra::writeRaster(v_crop, output_file, overwrite = TRUE)
   message(file, " processed")
   gc()
 }
@@ -145,38 +155,9 @@ for (file in filtered_files) {
 # list_global <- list_global %>% flatten_chr()
 vrt_list <- vrt_list %>% flatten_chr()
 
-#testing writing raster
-# for(i in seq_along(filtered_files)){
-#   message("starting ", i)
-#   x <- mask_raster(filtered_files[i], fields = pivot_sf_q)
-#   writeRaster(x, filtered_files[i], overwrite=TRUE)
-#   message("writing ", i)
-# }
-
-
-# rasters_masked[sapply(rasters_masked, is.null)] <- NULL
-
-# rasters_masked <- setNames(list_global, map_chr(names(list_global), date_extract))
-
-# })
-
-
 total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>%
-  crop(pivot_sf_q, mask = TRUE) %>%
+  crop(field_boundaries, mask = TRUE) %>%
   global(., fun="notNA") #%>%
-# as.matrix() %>%
-# `[`(1, 1)
-
-# total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(pivot_sf_q, mask = TRUE) %>% freq(., usenames = TRUE)
-
-# list_global
-# rast(list_global[1])[[5]] %>% plot()
-
-# vrt_files <- list.files(here(data_dir, "vrt"),full.names = T)
-# vrt_days_filter <- tools::file_path_sans_ext(basename(list_global))
-# vrt_list <- map(vrt_days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>%
-#   compact() %>%
-#   flatten_chr()
 
 layer_5_list <- purrr::map(vrt_list, function(vrt_list) {
   rast(vrt_list[1]) %>% terra::subset(1)
@@ -190,11 +171,6 @@ missing_pixels_count <- layer_5_list %>% global(., fun="notNA") %>%
     thres_40perc = as.integer(missing_pixels_percentage < 45)
     )
 
-# cloud_perc_list <- freq(layer_5_list, usenames = TRUE) %>%
-#   mutate(cloud_perc = (100 -((count/sum(total_pix_area$notNA))*100)),
-#          cloud_thres_5perc = as.integer(cloud_perc < 5),
-#          cloud_thres_40perc = as.integer(cloud_perc < 40)) %>%
-#   rename(Date = layer) %>% select(-value, -count)
 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))
 
@@ -208,14 +184,14 @@ if(sum(missing_pixels_count$thres_5perc)>1){
   rsrc <- sprc(cloudy_rasters_list)
   x <- mosaic(rsrc, fun = "max")
   
-  names(x) <-  "CI"
-  
+  # 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("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")
   
@@ -223,26 +199,26 @@ if(sum(missing_pixels_count$thres_5perc)>1){
   
   rsrc <- sprc(cloudy_rasters_list)
   x <- mosaic(rsrc, fun = "max")
-  names(x) <-  "CI"
-  
+  # 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("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("CI")
+  names(x) <-  c("Red", "Green", "Blue", "NIR", "CI")
 }
 
-plot(x$CI, main = paste("CI map", dates$week))
-#plotRGB(x, main = paste("RGB image week", dates$week))
+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)
@@ -257,7 +233,7 @@ extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_di
   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, field_geojson, fun = "mean"), 2)) %>%
+  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"))
@@ -265,33 +241,31 @@ extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_di
 }
 
 # 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")
+# 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)
+# 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")
+# )
 
-pivots_data_present <- unique(pivots_dates0$pivot_quadrant)
-quadrant_list <- pivots_data_present
+# 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.
 if (new_project_question == TRUE) {
 
   print("combined_CI_data.rds does not exist. Preparing combined_CI_data.rds file for all available images.")
 
-  walk(raster_files_NEW, extract_rasters_daily, field_geojson= pivot_sf_q, quadrants = TRUE, daily_CI_vals_dir)
+  walk(raster_files_NEW, extract_rasters_daily, field_geojson= field_boundaries, quadrants = F, daily_CI_vals_dir)
 
   extracted_values <-  list.files(here(daily_CI_vals_dir), full.names = TRUE)
 
@@ -313,7 +287,7 @@ if (new_project_question == TRUE) {
     compact() %>%
     flatten_chr()
 
-  walk(filtered_files, extract_rasters_daily, field_geojson= pivot_sf_q, quadrants = TRUE, daily_CI_vals_dir)
+  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)
 
@@ -326,34 +300,13 @@ if (new_project_question == TRUE) {
     group_by(pivot_quadrant) %>%
     summarise(across(everything(), ~ first(na.omit(.))))
 
-  combined_CI_data <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) %>% drop_na(pivot_quadrant)
+  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.")
 
 }
 
-
-
-
-#extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>%
-#  compact() %>%
-#  flatten_chr()
-
-#combine them into one df
-
-
-
- #saveRDS(pivot_stats, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file
-
-
-#load historic CI data and update it with last week of CI data
-# 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 <- purrr::map(list.files(here(daily_CI_vals_dir), full.names = TRUE, pattern = "quadrant"), readRDS) %>% list_rbind() %>% group_by(pivot_quadrant) %>%
-#   summarise(across(everything(), ~ first(na.omit(.))))
-
-
-
 # gather data into long format for easier calculation and visualisation
 pivot_stats_long <- pivot_stats2 %>%
   tidyr::gather("Date", value, -pivot_quadrant, -pivot ) %>%
@@ -362,36 +315,73 @@ pivot_stats_long <- pivot_stats2 %>%
   ) %>%
   drop_na(c("value","Date")) %>%
   mutate(value = as.numeric(value))%>%
-  filter_all(all_vars(!is.infinite(.)))%>%
-  rename(Field = pivot_quadrant) %>%
+  filter_all(all_vars(!is.infinite(.))) %>%
+  rename(Field = pivot_quadrant,
+         subField = Field) %>%
   unique()
 
-# #2021
-pivots_dates_Data_2021 <- pivots_dates0 %>% filter(!is.na(season_start_2021))
-pivot_select_model_Data_2021 <- unique(pivots_dates_Data_2021$pivot_quadrant)
 
-# #2022
-pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022))
-pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant )
-# #2023
-pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023))
-pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant)
-# #2024
-pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024))
-pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant)
+# #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) %>% pull(Field)
+
+pivot_select_model_Data_2024 <- harvesting_data %>% filter(Year == 2024) %>% pull(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.1A"
+  #  harvesting_data = harvesting_data
+  # field_CI_data = pivot_stats_long
+  #  season= 2023
+  
+  filtered_harvesting_data <- harvesting_data %>%
+    filter(Year == season, subField %in% field_names)
+  
+  filtered_field_CI_data <- field_CI_data %>%
+    filter(subField %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,
+           subField = 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_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_2022, Data_2023, Data_2024)
+CI_all <- rbind(Data_2023, Data_2024)
 message('CI_all created')
 #CI_all <-  Data_2023
 
@@ -403,4 +393,4 @@ 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')
+message('rds saved')

r_app/CI_report_dashboard_planet.Rmd

@@ -43,13 +43,17 @@ library(rsample)
 library(caret)
 library(randomForest)
 library(CAST)
+
+
+
 ```
 
 ```{r directories, message=FALSE, warning=FALSE, include=FALSE}
 laravel_storage_dir <- here("laravel_app/storage/app/",params$data_dir)
+laravel_storage_dir <- here("laravel_app/storage/app",project_dir)
 
 data_dir_project <- here(laravel_storage_dir, "Data")
-message('DATA_DIR',data_dir_project)
+# message('DATA_DIR',data_dir_project)
 extracted_CI_dir <- here(data_dir_project, "extracted_ci")
 daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals")
 cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals")
@@ -57,8 +61,8 @@ harvest_dir <-  here(data_dir_project, "HarvestData")
 
 weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
 
-
-s2_dir <- "C:/Users/timon/Resilience BV/4002 CMD App - General/4002 CMD Team/4002 TechnicalData/04 WP2 technical/python/chemba_S2/"
+source(here("r_app", "parameters_project.R"))
+# s2_dir <- "C:/Users/timon/Resilience BV/4002 CMD App - General/4002 CMD Team/4002 TechnicalData/04 WP2 technical/python/chemba_S2/"
 ```
 
 
@@ -119,21 +123,22 @@ CI_m3 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_3, "_", year_3,
 last_week_dif_raster_abs <- (CI - CI_m1)
 three_week_dif_raster_abs <- (CI - CI_m3)
 
-AllPivots0 <-st_read(here(data_dir_project, "pivot.geojson"))
+# AllPivots0 <-st_read(here(data_dir_project, "pivot.geojson"))
+
 # AllPivots0$pivot  <- factor(AllPivots0$pivot, levels = 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"))
-AllPivots0
+AllPivots0 <- field_boundaries_sf
 
-joined_spans <-st_read(here(data_dir_project, "span.geojson")) %>% st_transform(crs(AllPivots0))
+# joined_spans <-st_read(here(data_dir_project, "span.geojson")) %>% st_transform(crs(AllPivots0))
 
-pivots_dates <- readRDS(here(harvest_dir, "harvest_data_new"))
+# pivots_dates <- readRDS(here(harvest_dir, "harvest_data_new"))
 # pivots_dates$pivot  <- factor(pivots_dates$pivot, levels = 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"))
 
-AllPivots <- merge(AllPivots0, pivots_dates, by = "pivot_quadrant") %>%
-  rename(pivot = pivot.x) #%>% select(-pivot.y)
+AllPivots <- merge(AllPivots0, harvesting_data, by = c("Field", "subField")) #%>%
+  # rename(pivot = pivot.x) #%>% select(-pivot.y)
   head(AllPivots)
 
-AllPivots_merged <- AllPivots %>%
-   group_by(pivot) %>% summarise()
+AllPivots_merged <- AllPivots %>% 
+   group_by(Field) %>% summarise()
 
 AllPivots_merged <- st_transform(AllPivots_merged, crs = proj4string(CI))
 
@@ -208,7 +213,7 @@ create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend =
   tm_layout(main.title = paste0("\nMax CI week ", week,"\n",  age, " weeks old"),
             main.title.size = 0.7, legend.show = show_legend) +
   tm_shape(pivot_shape) +
-  tm_borders(lwd = 3) +  tm_text("pivot_quadrant", size = 1/2) +
+  tm_borders(lwd = 3) +  tm_text("subField", size = 1/2) +
   tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
 }
 
@@ -218,16 +223,16 @@ create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_lege
   tm_layout(main.title =  paste0("CI change week ", week_1, "- week ",week_2, "\n", age," weeks old"),
             main.title.size = 0.7, legend.show = show_legend) +
   tm_shape(pivot_shape) +
-  tm_borders(lwd = 3) +  tm_text("pivot_quadrant", size = 1/2) +
+  tm_borders(lwd = 3) +  tm_text("subField", size = 1/2) +
   tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
 }
 
 ci_plot <- function(pivotName){
    # pivotName = "1.1"
-   pivotShape <- AllPivots_merged %>% terra::subset(pivot %in% pivotName)  %>% st_transform(crs(CI))
-   age <- AllPivots %>% dplyr::filter(pivot %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique()
+   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()
 
-   AllPivots2 <- AllPivots %>% dplyr::filter(pivot %in% pivotName)
+   AllPivots2 <- AllPivots %>% dplyr::filter(Field %in% pivotName)
 
    singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape)
 
@@ -238,9 +243,9 @@ ci_plot <- function(pivotName){
    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 <- pivots_dates %>% dplyr::filter(pivot %in% pivotName) %>% ungroup()  %>% dplyr::select(planting_date)  %>% unique()
+   planting_date <- pivots_dates %>% dplyr::filter(Field %in% pivotName) %>% ungroup()  %>% dplyr::select(planting_date)  %>% unique()
 
-   joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape))  %>%  dplyr::filter(pivot %in% pivotName) %>% st_crop(., pivotShape)
+   joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape))  %>%  dplyr::filter(Field %in% pivotName) %>% unique() %>% 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)
    CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= F, legend_is_portrait = F, week = week_minus_1, age = age  -1)
@@ -266,7 +271,7 @@ ci_plot <- function(pivotName){
 cum_ci_plot <- function(pivotName){
 
   # pivotName = "1.17"
-  data_ci <- CI_quadrant %>% filter(pivot == pivotName)
+  data_ci <- CI_quadrant %>% filter(Field == pivotName)
   data_ci2 <- data_ci  %>% mutate(CI_rate = cumulative_CI/DOY,
                                   week = week(Date))%>% group_by(Field) %>%
     mutate(mean_rolling10 = rollapplyr(CI_rate  , width = 10, FUN = mean, partial = TRUE)) #%>%
@@ -275,38 +280,38 @@ cum_ci_plot <- function(pivotName){
                                                                                 max_date = max(Date),
                                                                                 days = max_date - min_date)
 
-  perfect_pivot_raw <- CI_quadrant %>% group_by(pivot) %>% filter(pivot == "5.1" & season == "2022") %>%
-    group_by(DOY) %>% summarise(cumulative_CI = mean(cumulative_CI)) %>% mutate(CI_rate = cumulative_CI/DOY) %>%
-    mutate(mean_rolling10 = rollapplyr(CI_rate  , width = 10, FUN = mean, partial = TRUE))
+  # perfect_pivot_raw <- CI_quadrant %>% group_by(Field) %>% filter(Field == "5.1" & season == "2022") %>%
+  #   group_by(DOY) %>% summarise(cumulative_CI = mean(cumulative_CI)) %>% mutate(CI_rate = cumulative_CI/DOY) %>%
+  #   mutate(mean_rolling10 = rollapplyr(CI_rate  , width = 10, FUN = mean, partial = TRUE))
 
   unique_seasons <- unique(date_preperation_perfect_pivot$season)
 
-   if(length(unique_seasons) == 3) {
-     unique_seasons <- unique_seasons[c(2,3)]
-   } else if(length(unique_seasons) == 4) {
-     unique_seasons <- unique_seasons[c(3,4)]
-   } else {
-     unique_seasons <- unique_seasons
-   }
+   # if(length(unique_seasons) == 3) {
+   #   unique_seasons <- unique_seasons[c(2,3)]
+   # } else if(length(unique_seasons) == 4) {
+   #   unique_seasons <- unique_seasons[c(3,4)]
+   # } else {
+   #   unique_seasons <- unique_seasons
+   # }
 
-  perfect_pivot <- perfect_pivot_raw
+  # perfect_pivot <- perfect_pivot_raw
+  # 
+  # for (s in unique_seasons) {
+  #   season_dates <- seq(from = date_preperation_perfect_pivot$min_date[date_preperation_perfect_pivot$season == s],
+  #                       to = date_preperation_perfect_pivot$min_date[date_preperation_perfect_pivot$season == s] + nrow(perfect_pivot_raw) - 1,
+  #                       by = "1 day")
+  #   col_name <- as.character(s)
+  #   perfect_pivot <- dplyr::bind_cols(perfect_pivot, tibble(!!col_name := season_dates))
+  # }
+  # 
+  # perfect_pivot <- perfect_pivot %>%
+  #   pivot_longer(cols = -c("DOY", "cumulative_CI", "CI_rate", "mean_rolling10"),
+  #                names_to = "season", values_to = "Date")
 
-  for (s in unique_seasons) {
-    season_dates <- seq(from = date_preperation_perfect_pivot$min_date[date_preperation_perfect_pivot$season == s],
-                        to = date_preperation_perfect_pivot$min_date[date_preperation_perfect_pivot$season == s] + nrow(perfect_pivot_raw) - 1,
-                        by = "1 day")
-    col_name <- as.character(s)
-    perfect_pivot <- dplyr::bind_cols(perfect_pivot, tibble(!!col_name := season_dates))
-  }
-
-  perfect_pivot <- perfect_pivot %>%
-    pivot_longer(cols = -c("DOY", "cumulative_CI", "CI_rate", "mean_rolling10"),
-                 names_to = "season", values_to = "Date")
-
-  g <- ggplot() +
+  g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) +
           facet_wrap(~season, scales = "free_x") +
-          geom_line(data= data_ci2 %>% filter(season %in% unique_seasons), aes(Date, mean_rolling10, col = Field)) +
-          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) +
+          geom_line( aes(Date, mean_rolling10, col = Field)) +
+          # 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("14 day rolling MEAN CI rate - Pivot ", pivotName))+
           # 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") +
@@ -348,7 +353,7 @@ tm_shape(CI, unit = "m")+
 
    tm_compass(position = c("right", "top"), text.color = "black") +
   tm_shape(AllPivots)+ tm_borders( col = "black") +
-    tm_text("pivot_quadrant", size = .6, col = "black")
+    tm_text("subField", size = .6, col = "black")
 
 ```
 \newpage
@@ -362,7 +367,7 @@ tm_shape(CI, unit = "m")+
 
     tm_compass(position = c("right", "top"), text.color = "black") +
    tm_shape(AllPivots)+ tm_borders( col = "black") +
-     tm_text("pivot_quadrant", size = .6, col = "black")
+     tm_text("subField", size = .6, col = "black")
 
 ```
 
@@ -371,7 +376,9 @@ tm_shape(CI, unit = "m")+
 
 ```{r plots_ci_estate,  echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis'}
 # # pivots <- AllPivots_merged %>% filter(pivot != c("1.1", "1.17"))
-pivots_estate <- AllPivots_merged # %>% 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" , "6.1", "6.2", "DL1.1", "DL1.3")) %>% filter(pivot != "1.17")
+pivots_estate <- AllPivots_merged  %>% 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" , "6.1", "6.2", "DL1.1", "DL1.3")) %>% filter(pivot != "1.17")
+
+
 # pivots <- AllPivots_merged %>% filter(pivot != c("1.1", "1.17"))
 # pivots_estate <- AllPivots_merged %>% filter(pivot %in% c("1.1", "1.2", "1.7")) %>% filter(pivot != "1.17")
 
@@ -383,6 +390,22 @@ walk(pivots_estate$pivot, ~ {
 ```
 
 
+```{r looping over sub_area,  echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis'}
+pivots_grouped <- AllPivots_merged %>% 
+  group_by(sub_area) %>%
+  arrange(sub_area)  # Optional: arrange the groups alphabetically by sub_area
+
+# Iterate over each subgroup
+for (subgroup in unique(pivots_grouped$sub_area)) {
+  cat("\n## Subgroup:", subgroup, "\n")  # Add a title for the subgroup
+  subset_data <- filter(pivots_grouped, sub_area == subgroup)
+  walk(subset_data$Field, ~ {
+    cat("\n")  # Add an empty line for better spacing
+    ci_plot(.x)
+    cum_ci_plot(.x)
+  })
+}
+```
 
 ```{r eval=FALSE, fig.height=10, fig.width=14, include=FALSE}
 CI_all2 <- readRDS(here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_whole_pivot_year.rds")) %>%

r_app/parameters_project.R

@@ -0,0 +1,93 @@
+
+
+#chemba
+if(project_dir == "chemba"){
+  message("Yield data for Chemba")
+
+  
+  field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant)%>% 
+    mutate(sub_area = case_when(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" , "6.1", "6.2", "DL1.1", "DL1.3") ~ "estate",
+                                TRUE ~ "Cooperative")) 
+  names(field_boundaries_sf) <- c("Field", "subField", "geometry", "sub_area")
+  
+  field_boundaries <- field_boundaries_sf %>% vect()
+  names(field_boundaries) <- c("Field", "subField", "sub_area")
+  
+  # field_boundaries_merged <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant)   %>% group_by(pivot) %>% summarise() %>% vect()
+  
+  joined_spans <-st_read(here(data_dir_project, "span.geojson")) %>% st_transform(crs(field_boundaries_sf))
+  names(joined_spans) <- c("Field", "area", "radius", "spans", "span", "geometry")
+  
+  
+  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 %>%
+    dplyr::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", "Age")) %>%
+    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(subField = pivot_quadrant) %>%
+    mutate(Field = substr(subField, 1, 3))
+  
+  
+  
+  
+  
+  
+  
+  
+} else if (project_dir == "xinavane"){
+  library(readxl)
+  message("Yield data for Xinavane")
+  
+  field_boundaries <- st_read(here(data_dir, "Xinavane_demo.geojson")) %>% vect() 
+  field_boundaries_sf <- st_read(here(data_dir, "Xinavane_demo.geojson")) 
+  
+  joined_spans <- field_boundaries
+  
+  pivots_dates0 <- read_excel(here(harvest_dir, "Yield data.xlsx"), 
+                              skip = 3,
+                              col_types = c("numeric", "text", "skip", "text", "numeric", "numeric", "numeric", "numeric", "date", 
+                              "numeric", "skip", "numeric")) %>%
+    rename(
+      Year = 1,
+      Field = 2,
+      sub_area = 3,
+      hectares = 4,
+      tons = 5,
+      tcha = 6,
+      tchy = 7, 
+      Season_end = 8,
+      age = 9,
+      ratoon = 10
+    ) %>%
+    mutate(Season_end = ymd(Season_end),
+           Season_start = as.Date(Season_end - (duration(months = age))),
+           subField = Field) #don't forget to add 2022 as a year for the 'curent' season
+  
+  pivots_dates0 <- pivots_dates0 %>%
+    mutate(Year = Year + 6)
+  
+  # Add 6 years to Season_end column
+  pivots_dates0 <- pivots_dates0 %>%
+    mutate(Season_end = Season_end + years(6))
+  
+  # Add 6 years to Season_start column
+  pivots_dates0 <- pivots_dates0 %>%
+    mutate(Season_start = Season_start + years(6))
+  
+  
+  harvesting_data <- pivots_dates0 %>% select(c("Field","subField", "Year", "Season_start","Season_end", "Age" ))
+
+  
+} else {
+  message("No yield data supplied")
+}
+