SmartCane/r_app/ci_extraction.R

# nolint start: commented_code_linter, line_length_linter,object_usage_linter.
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])
if (is.na(end_date)) {
end_date <- lubridate::dmy("28-08-2024")
}

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
}

week <- week(end_date)
# 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(here("r_app/parameters_project.R"))
# source(here("r_app/ci_extraction_utils.R"))
# source(here("r_app/mosaic_creation_utils.R"))

source("parameters_project.R")
source("ci_extraction_utils.R")
# source("mosaic_creation_utils.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)

dates <- date_list(end_date, offset)
print(dates)

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)

vrt_list <- list()

for (file in filtered_files) {
  v_crop <- create_mask_and_crop(file, field_boundaries, merged_final)
  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,] >  100){
    vrt_list[vrt_file] <- vrt_file

  }else{
    file.remove(vrt_file)
  }

  message(file, " processed")
  gc()
}

raster_files_NEW <- list.files(merged_final,full.names = T, pattern = ".tif")

# 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)) {
  # File does not exist, create it with all available data
  print("combined_CI_data.rds does not exist. Preparing combined_CI_data.rds file for all available images.")
  
  # Extract data from all raster files
  walk(raster_files_NEW, extract_rasters_daily, field_geojson = field_boundaries, quadrants = FALSE, daily_CI_vals_dir)
  
  # Combine all extracted data
  extracted_values <- list.files(here(daily_CI_vals_dir), full.names = TRUE)
  
  pivot_stats <- extracted_values %>%
    map(readRDS) %>% list_rbind() %>%
    group_by(sub_field) 
  
  # Save the combined data to the file
  saveRDS(pivot_stats, file_path)
  
  print("All CI values extracted from all historic images and saved to combined_CI_data.rds.")
  
} else {
  # File exists, add new data
  print("combined_CI_data.rds exists, adding the latest image data to the table.")
  
  # Filter and process the latest data
  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)
  
  # Extract new values
  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)
  
  # Load existing data and append new data
  combined_CI_data <- readRDS(file_path)
  pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data)
  
  # Save the updated combined data
  saveRDS(pivot_stats2, file_path)
  
  print("All CI values extracted from the latest images and added to combined_CI_data.rds.")
}