# CI_EXTRACTION.R # ============== # This script processes satellite imagery to extract Canopy Index (CI) values for agricultural fields. # It handles image processing, masking, and extraction of statistics by field/sub-field. # # Usage: Rscript ci_extraction.R [end_date] [offset] [project_dir] # - end_date: End date for processing (YYYY-MM-DD format) # - offset: Number of days to look back from end_date # - project_dir: Project directory name (e.g., "chemba") # # 1. Load required packages # ----------------------- suppressPackageStartupMessages({ library(sf) library(terra) library(tidyverse) library(lubridate) library(exactextractr) library(readxl) library(here) }) # 2. Process command line arguments # ------------------------------ main <- function() { # Capture command line arguments args <- commandArgs(trailingOnly = TRUE) # Process end_date argument if (length(args) >= 1 && !is.na(args[1])) { end_date <- as.Date(args[1]) if (is.na(end_date)) { warning("Invalid end_date provided. Using default (current date).") end_date <- Sys.Date() #end_date <- "2023-10-01" } } else { end_date <- Sys.Date() #end_date <- "2023-10-01" } # Process offset argument if (length(args) >= 2 && !is.na(args[2])) { offset <- as.numeric(args[2]) if (is.na(offset) || offset <= 0) { warning("Invalid offset provided. Using default (7 days).") offset <- 7 } } else { offset <- 7 } # Process project_dir argument if (length(args) >= 3 && !is.na(args[3])) { project_dir <- as.character(args[3]) } else { project_dir <- "chemba" } # Make project_dir available globally so parameters_project.R can use it assign("project_dir", project_dir, envir = .GlobalEnv) # 3. Initialize project configuration # -------------------------------- new_project_question <- FALSE tryCatch({ source("parameters_project.R") source("ci_extraction_utils.R") }, error = function(e) { warning("Default source files not found. Attempting to source from 'r_app' directory.") tryCatch({ source(here::here("r_app", "parameters_project.R")) source(here::here("r_app", "ci_extraction_utils.R")) warning(paste("Successfully sourced files from 'r_app' directory.")) }, error = function(e) { stop("Failed to source required files from both default and 'r_app' directories.") }) }) # 4. Generate date list for processing # --------------------------------- dates <- date_list(end_date, offset) log_message(paste("Processing data for week", dates$week, "of", dates$year)) # 5. Find and filter raster files by date # ----------------------------------- log_message("Searching for raster files") tryCatch({ # Use the new utility function to find satellite images existing_files <- find_satellite_images(planet_tif_folder, dates$days_filter) log_message(paste("Found", length(existing_files), "raster files for processing")) # 6. Process raster files and create VRT # ----------------------------------- # Use the new utility function for batch processing vrt_list <- process_satellite_images(existing_files, field_boundaries, merged_final, daily_vrt) # 7. Process and combine CI values # ------------------------------ # Call the process_ci_values function from utils with all required parameters process_ci_values(dates, field_boundaries, merged_final, field_boundaries_sf, daily_CI_vals_dir, cumulative_CI_vals_dir) }, error = function(e) { log_message(paste("Error in main processing:", e$message), level = "ERROR") stop(e$message) }) } if (sys.nframe() == 0) { main() }