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

library(CIprep)


laravel_storage_dir <- here("laravel_app/storage/app")
#preparing directories
planet_tif_folder <- here(laravel_storage_dir, "chemba/merged_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(data_dir, "weekly_mosaic")
harvest_dir <-  here(data_dir, "HarvestData")

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))




# Creating weekly mosaic
dates <- date_list(0)

#load pivot geojson
pivot_sf_q <- st_read(here(geojson_dir, "pivot_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()

#load and filter raster files
raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif") #use pattern = '.tif$' or something else if you have multiple files in this folder
filtered_files  <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
  compact() %>%
  flatten_chr()

rasters_masked <- map(filtered_files, mask_raster, fields = pivot_sf_q) %>% set_names(filtered_files)

rasters_masked[sapply(rasters_masked, is.null)] <- NULL

rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract))

total_pix_area <- rast(filtered_files[1]) %>% terra::subset(1) %>%
  crop(pivot_sf_q, 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)



layer_5_list <- purrr::map(rasters_masked, function(spatraster) {
  spatraster[[5]]
}) %>% rast()

cloud_perc_list <- freq(layer_5_list, usenames = TRUE) %>%
  mutate(cloud_perc = (100 -((count/sum(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))

## Create mosaic

if(sum(cloud_perc_list$cloud_thres_5perc)>1){
  message("More than 1 raster without clouds (<5%), max composite made")

  cloudy_rasters_list <- rasters_masked[cloud_index_5perc]
  raster_list_subset <- lapply(cloudy_rasters_list, function(r) {
    subset(r, 6)})

  rsrc <- sprc(raster_list_subset)
  x <- mosaic(rsrc, fun = "max")

  # x <- rast(rasters_masked[cloud_index_5perc[1]])
  names(x) <-  "CI"
  # writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE)
  # message("raster exported")

}else if(sum(cloud_perc_list$cloud_thres_5perc)==1){
  message("Only 1 raster without clouds (<5%)")

  x <- rast(rasters_masked[cloud_index_5perc[1]])
  names(x) <-  c("red", "green", "blue","nir", "cloud" ,"CI")

  rsrc <- sprc(raster_list_subset)
  x <- mosaic(rsrc, fun = "max")

  # x <- rast(rasters_masked[cloud_index_5perc[1]])
  names(x) <-  "CI"
  # writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE)
  # message("raster exported")

}else if(sum(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_index_40perc]
  raster_list_subset <- lapply(cloudy_rasters_list, function(r) {
    subset(r, 6)
  })

  rsrc <- sprc(raster_list_subset)
  x <- mosaic(rsrc, fun = "max")
  names(x) <-  "CI"
  # writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE)
  # message("raster exported")

}else{
  message("No cloud free images available")
  x <-  rast(rasters_masked[1])
  x[x] <- NA
  names(x) <-  c("red", "green", "blue","nir", "cloud" ,"CI")
  # writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE)
  # message("raster exported")

}

plot(x$CI, main = paste("CI map", dates$week))
#plotRGB(x, main = paste("RGB image week", dates$week))

# writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE)
writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE)



# 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_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant)   %>% group_by(pivot) %>% summarise() %>% vect()

walk(filtered_files, extract_rasters_daily, field_geojson= pivot_sf_q, quadrants = TRUE, daily_CI_vals_dir)

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_dates_long <- pivots_dates0 %>%
#   select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023")) %>%
#   pivot_longer(cols = c("season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023")) %>%
#   separate(pivot_quadrant, into = c("name", "Year"), sep = "\\.")

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")) %>%
  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)
# extracted_values <-  list.files("C:\\Users\\timon\\Resilience BV\\4002 CMD App - General\\4002 CMD Team\\4002 TechnicalData\\04 WP2 technical\\DetectingSpotsR\\EcoFarm\\planet\\extracted",
#                                      pattern ="_quadrant", full.names = TRUE)
extracted_values <-  list.files(here(daily_CI_vals_dir), full.names = TRUE)

head(extracted_values)

pivot_stats <- extracted_values %>%
  map_dfr(readRDS) %>%
  group_by(pivot_quadrant) %>%
  summarise(across(everything(), ~ first(na.omit(.))))

pivots_data_present <- unique(pivots_dates0$pivot_quadrant)
quadrant_list <- pivots_data_present

# gather data into long format for easier calculation and visualisation
pivot_stats_long <- pivot_stats %>%
  tidyr::gather("Date", value, -pivot_quadrant, -pivot ) %>%
  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)

# #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)


## Extracting the correct CI values

Data_2021 <- map_dfr(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021))
Data_2022 <- map_dfr(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022))
Data_2023 <- map_dfr(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023))
#Data_2024 <- map_dfr(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024))


CI_all <- rbind(Data_2021, Data_2022, Data_2023)

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

head(CI_all)


saveRDS(CI_all, here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))