From 054cc85bdbc55d5e8fd9fb519d666e573206f98d Mon Sep 17 00:00:00 2001 From: Timon Date: Tue, 10 Feb 2026 15:38:23 +0100 Subject: [PATCH] small issues --- r_app/80_calculate_kpis.R | 1154 ++++++++++----------- r_app/80_utils_agronomic_support.R | 128 +-- r_app/80_utils_common.R | 8 +- r_app/90_CI_report_with_kpis_simple.Rmd | 2 +- r_app/DEBUG_remove_date_tiffs.R | 272 +++++ r_app/FIX_INDENTATION.R | 20 + r_app/MANUAL_PIPELINE_RUNNER.R | 3 +- r_app/extract_rds_only.R | 104 -- r_app/old_scripts/kpi_utils.R | 7 +- r_app/parameters_project_OLD.R | 1240 ----------------------- r_app/report_utils.R | 3 +- 11 files changed, 956 insertions(+), 1985 deletions(-) create mode 100644 r_app/DEBUG_remove_date_tiffs.R create mode 100644 r_app/FIX_INDENTATION.R delete mode 100644 r_app/extract_rds_only.R delete mode 100644 r_app/parameters_project_OLD.R diff --git a/r_app/80_calculate_kpis.R b/r_app/80_calculate_kpis.R index 8b76ac0..5d1c415 100644 --- a/r_app/80_calculate_kpis.R +++ b/r_app/80_calculate_kpis.R @@ -397,627 +397,649 @@ main <- function() { message("\n", strrep("-", 70)) message("PHASE 1: PER-FIELD WEEKLY ANALYSIS ") message(strrep("-", 70)) - weeks <- calculate_week_numbers(end_date) - current_week <- weeks$current_week - current_year <- weeks$current_year - previous_week <- weeks$previous_week - previous_year <- weeks$previous_year + + weeks <- calculate_week_numbers(end_date) + current_week <- weeks$current_week + current_year <- weeks$current_year + previous_week <- weeks$previous_week + previous_year <- weeks$previous_year + + message(paste("Week:", current_week, "/ Year (ISO 8601):", current_year)) + # Find per-field weekly mosaics + message("Finding per-field weekly mosaics...") - message(paste("Week:", current_week, "/ Year (ISO 8601):", current_year)) - # Find per-field weekly mosaics - message("Finding per-field weekly mosaics...") - - if (!dir.exists(weekly_mosaic)) { - stop(paste("ERROR: weekly_mosaic directory not found:", weekly_mosaic, - "\nScript 40 (mosaic creation) must be run first.")) - } - - field_dirs <- list.dirs(weekly_mosaic, full.names = FALSE, recursive = FALSE) - field_dirs <- field_dirs[field_dirs != ""] - - if (length(field_dirs) == 0) { - stop(paste("ERROR: No field subdirectories found in:", weekly_mosaic, - "\nScript 40 must create weekly_mosaic/{FIELD}/ structure.")) - } - - # Verify we have mosaics for this week - single_file_pattern <- sprintf("week_%02d_%d\\.tif", current_week, current_year) - per_field_files <- c() - for (field in field_dirs) { - field_mosaic_dir <- file.path(weekly_mosaic, field) - files <- list.files(field_mosaic_dir, pattern = single_file_pattern, full.names = TRUE) - if (length(files) > 0) { - per_field_files <- c(per_field_files, files) + if (!dir.exists(weekly_mosaic)) { + stop(paste("ERROR: weekly_mosaic directory not found:", weekly_mosaic, + "\nScript 40 (mosaic creation) must be run first.")) } - } - if (length(per_field_files) == 0) { - stop(paste("ERROR: No mosaics found for week", current_week, "year", current_year, - "\nExpected pattern:", single_file_pattern, - "\nChecked:", weekly_mosaic)) - } + field_dirs <- list.dirs(weekly_mosaic, full.names = FALSE, recursive = FALSE) + field_dirs <- field_dirs[field_dirs != ""] - message(paste(" ✓ Found", length(per_field_files), "per-field weekly mosaics")) + if (length(field_dirs) == 0) { + stop(paste("ERROR: No field subdirectories found in:", weekly_mosaic, + "\nScript 40 must create weekly_mosaic/{FIELD}/ structure.")) + } - mosaic_mode <- "per-field" - mosaic_dir <- weekly_mosaic + # Verify we have mosaics for this week + single_file_pattern <- sprintf("week_%02d_%d\\.tif", current_week, current_year) + per_field_files <- c() + for (field in field_dirs) { + field_mosaic_dir <- file.path(weekly_mosaic, field) + files <- list.files(field_mosaic_dir, pattern = single_file_pattern, full.names = TRUE) + if (length(files) > 0) { + per_field_files <- c(per_field_files, files) + } + } + + if (length(per_field_files) == 0) { + stop(paste("ERROR: No mosaics found for week", current_week, "year", current_year, + "\nExpected pattern:", single_file_pattern, + "\nChecked:", weekly_mosaic)) + } + + message(paste(" ✓ Found", length(per_field_files), "per-field weekly mosaics")) + + mosaic_mode <- "per-field" + mosaic_dir <- weekly_mosaic - # Load field boundaries - tryCatch({ - boundaries_result <- load_field_boundaries(data_dir) - - if (is.list(boundaries_result) && "field_boundaries_sf" %in% names(boundaries_result)) { - field_boundaries_sf <- boundaries_result$field_boundaries_sf - } else { - field_boundaries_sf <- boundaries_result - } - - if (nrow(field_boundaries_sf) == 0) { - stop("No fields loaded from boundaries") - } - - message(paste(" Loaded", nrow(field_boundaries_sf), "fields")) - }, error = function(e) { - stop("ERROR loading field boundaries: ", e$message) - }) - - message("Loading historical field data for trend calculations...") - # Load up to 8 weeks (max of 4-week and 8-week trend requirements) - # Function gracefully handles missing weeks and loads whatever exists - num_weeks_to_load <- max(WEEKS_FOR_FOUR_WEEK_TREND, WEEKS_FOR_CV_TREND_LONG) # Always 8 - message(paste(" Attempting to load up to", num_weeks_to_load, "weeks of historical data...")) - - # Only auto-generate on first call (not in recursive calls from within load_historical_field_data) - allow_auto_gen <- !exists("_INSIDE_AUTO_GENERATE", envir = .GlobalEnv) - - historical_data <- load_historical_field_data(project_dir, current_week, year, reports_dir, - num_weeks = num_weeks_to_load, - auto_generate = allow_auto_gen, - field_boundaries_sf = field_boundaries_sf, - daily_vals_dir = daily_vals_dir) - - # Load harvest.xlsx for planting dates (season_start) - message("\nLoading harvest data from harvest.xlsx for planting dates...") - harvest_file_path <- file.path(data_dir, "harvest.xlsx") - - harvesting_data <- tryCatch({ - if (file.exists(harvest_file_path)) { - harvest_raw <- readxl::read_excel(harvest_file_path) - harvest_raw$season_start <- as.Date(harvest_raw$season_start) - harvest_raw$season_end <- as.Date(harvest_raw$season_end) - message(paste(" ✓ Loaded harvest data:", nrow(harvest_raw), "rows")) - harvest_raw - } else { - message(paste(" WARNING: harvest.xlsx not found at", harvest_file_path)) - NULL - } - }, error = function(e) { - message(paste(" ERROR loading harvest.xlsx:", e$message)) - NULL - }) - - planting_dates <- extract_planting_dates(harvesting_data, field_boundaries_sf) - - # Validate planting_dates - if (is.null(planting_dates) || nrow(planting_dates) == 0) { - message("WARNING: No planting dates available. Using NA for all fields.") - planting_dates <- data.frame( - field_id = field_boundaries_sf$field, - date = rep(as.Date(NA), nrow(field_boundaries_sf)), - stringsAsFactors = FALSE - ) - } - - # Build per-field configuration - message("\nPreparing mosaic configuration for statistics calculation...") - message(" ✓ Using per-field mosaic architecture (1 TIFF per field)") - - # Per-field mode: each field has its own TIFF in weekly_mosaic/{FIELD}/week_*.tif - field_grid <- list( - mosaic_dir = mosaic_dir, - mode = "per-field" - ) - - message("\nUsing modular RDS-based approach for weekly statistics...") - - # Load/calculate CURRENT week stats (from tiles or RDS cache) - message("\n1. Loading/calculating CURRENT week statistics (week", current_week, ")...") - current_stats <- load_or_calculate_weekly_stats( - week_num = current_week, - year = current_iso_year, - project_dir = project_dir, - field_boundaries_sf = field_boundaries_sf, - mosaic_dir = field_grid$mosaic_dir, - reports_dir = reports_dir, - report_date = end_date - ) - - message(paste(" ✓ Loaded/calculated stats for", nrow(current_stats), "fields in current week")) - - # Load/calculate PREVIOUS week stats (from RDS cache or tiles) - message("\n2. Loading/calculating PREVIOUS week statistics (week", previous_week, ")...") - - # Calculate report date for previous week (7 days before current) - prev_report_date <- end_date - 7 - - prev_stats <- load_or_calculate_weekly_stats( - week_num = previous_week, - year = previous_year, - project_dir = project_dir, - field_boundaries_sf = field_boundaries_sf, - mosaic_dir = field_grid$mosaic_dir, - reports_dir = reports_dir, - report_date = prev_report_date - ) - - message(paste(" ✓ Loaded/calculated stats for", nrow(prev_stats), "fields in previous week")) - message(paste(" Columns in prev_stats:", paste(names(prev_stats), collapse = ", "))) - message(paste(" CV column non-NA values in prev_stats:", sum(!is.na(prev_stats$CV)))) - - # Apply trend calculations (requires both weeks) - message("\n3. Calculating trend columns...") - current_stats <- calculate_kpi_trends(current_stats, prev_stats, - project_dir = project_dir, - reports_dir = reports_dir, - current_week = current_week, - year = current_iso_year) - - message(paste(" ✓ Added Weekly_ci_change, CV_Trend_Short_Term, Four_week_trend, CV_Trend_Long_Term, nmr_of_weeks_analysed")) - - # Load weekly harvest probabilities from script 31 (if available) - # Note: Script 31 saves to reports/kpis/field_stats/ (not field_level) - message("\n4. Loading harvest probabilities from script 31...") - harvest_prob_dir <- file.path(data_dir, "..", "reports", "kpis", "field_stats") - harvest_prob_file <- file.path(harvest_prob_dir, - sprintf("%s_harvest_imminent_week_%02d_%d.csv", project_dir, current_week, year)) - message(paste(" Looking for:", harvest_prob_file)) - - imminent_prob_data <- tryCatch({ - if (file.exists(harvest_prob_file)) { - prob_df <- readr::read_csv(harvest_prob_file, show_col_types = FALSE) - message(paste(" ✓ Loaded harvest probabilities for", nrow(prob_df), "fields")) - prob_df %>% - select(field, imminent_prob, detected_prob) %>% - rename(Field_id = field, Imminent_prob_actual = imminent_prob, Detected_prob = detected_prob) - } else { - message(paste(" INFO: Harvest probabilities not available (script 31 not run)")) - NULL - } - }, error = function(e) { - message(paste(" WARNING: Could not load harvest probabilities:", e$message)) - NULL - }) - - # ============================================================================ - # CALCULATE GAP FILLING KPI (2σ method from kpi_utils.R) - # ============================================================================ - - message("\nCalculating gap filling scores (2σ method)...") - - # Try single merged mosaic first, then fall back to merging tiles - week_mosaic_file <- file.path(mosaic_dir, sprintf("week_%02d_%d.tif", current_week, current_iso_year)) - - gap_scores_df <- NULL - - if (file.exists(week_mosaic_file)) { - # Single merged mosaic exists - use it directly + # Load field boundaries tryCatch({ - current_week_raster <- terra::rast(week_mosaic_file) - current_ci_band <- current_week_raster[[5]] # CI is the 5th band - names(current_ci_band) <- "CI" - - message(paste(" Loaded single mosaic:", week_mosaic_file)) - - # Calculate gap scores for all fields - gap_result <- calculate_gap_filling_kpi(current_ci_band, field_boundaries_sf) - - # Extract field-level results (use field column directly to match current_stats Field_id) - gap_scores_df <- gap_result$field_results %>% - mutate(Field_id = field) %>% - select(Field_id, gap_score) - - message(paste(" ✓ Calculated gap scores for", nrow(gap_scores_df), "fields")) - message(paste(" Gap score range:", round(min(gap_scores_df$gap_score, na.rm=TRUE), 2), "-", round(max(gap_scores_df$gap_score, na.rm=TRUE), 2), "%")) - + boundaries_result <- load_field_boundaries(data_dir) + + if (is.list(boundaries_result) && "field_boundaries_sf" %in% names(boundaries_result)) { + field_boundaries_sf <- boundaries_result$field_boundaries_sf + } else { + field_boundaries_sf <- boundaries_result + } + + if (nrow(field_boundaries_sf) == 0) { + stop("No fields loaded from boundaries") + } + + message(paste(" Loaded", nrow(field_boundaries_sf), "fields")) }, error = function(e) { - message(paste(" WARNING: Could not calculate gap scores from single mosaic:", e$message)) - message(" Gap scores will be set to NA") - gap_scores_df <- NULL + stop("ERROR loading field boundaries: ", e$message) }) - - } else { - # Single mosaic doesn't exist - check for tiles and process per-tile - message(" Single mosaic not found. Checking for tiles...") - - # List all tiles for this week (e.g., week_04_2026_01.tif through week_04_2026_25.tif) - tile_pattern <- sprintf("week_%02d_%d_\\d{2}\\.tif$", current_week, current_iso_year) - tile_files <- list.files(mosaic_dir, pattern = tile_pattern, full.names = TRUE) - - if (length(tile_files) == 0) { - message(sprintf(" WARNING: No tiles found matching pattern: %s in %s", tile_pattern, mosaic_dir)) - message(" Gap scores will be set to NA") - - } else { + + message("Loading historical field data for trend calculations...") + # Load up to 8 weeks (max of 4-week and 8-week trend requirements) + # Function gracefully handles missing weeks and loads whatever exists + num_weeks_to_load <- max(WEEKS_FOR_FOUR_WEEK_TREND, WEEKS_FOR_CV_TREND_LONG) # Always 8 + message(paste(" Attempting to load up to", num_weeks_to_load, "weeks of historical data...")) + + # Only auto-generate on first call (not in recursive calls from within load_historical_field_data) + allow_auto_gen <- !exists("_INSIDE_AUTO_GENERATE", envir = .GlobalEnv) + + historical_data <- load_historical_field_data(project_dir, current_week, current_year, reports_dir, + num_weeks = num_weeks_to_load, + auto_generate = allow_auto_gen, + field_boundaries_sf = field_boundaries_sf, + daily_vals_dir = daily_vals_dir) + + # Load harvest.xlsx for planting dates (season_start) + message("\nLoading harvest data from harvest.xlsx for planting dates...") + harvest_file_path <- file.path(data_dir, "harvest.xlsx") + + harvesting_data <- tryCatch({ + if (file.exists(harvest_file_path)) { + harvest_raw <- readxl::read_excel(harvest_file_path) + harvest_raw$season_start <- as.Date(harvest_raw$season_start) + harvest_raw$season_end <- as.Date(harvest_raw$season_end) + message(paste(" ✓ Loaded harvest data:", nrow(harvest_raw), "rows")) + harvest_raw + } else { + message(paste(" WARNING: harvest.xlsx not found at", harvest_file_path)) + NULL + } + }, error = function(e) { + message(paste(" ERROR loading harvest.xlsx:", e$message)) + NULL + }) + + planting_dates <- extract_planting_dates(harvesting_data, field_boundaries_sf) + + # Validate planting_dates + if (is.null(planting_dates) || nrow(planting_dates) == 0) { + message("WARNING: No planting dates available. Using NA for all fields.") + planting_dates <- data.frame( + field_id = field_boundaries_sf$field, + date = rep(as.Date(NA), nrow(field_boundaries_sf)), + stringsAsFactors = FALSE + ) + } + + # Build per-field configuration + message("\nPreparing mosaic configuration for statistics calculation...") + message(" ✓ Using per-field mosaic architecture (1 TIFF per field)") + + # Per-field mode: each field has its own TIFF in weekly_mosaic/{FIELD}/week_*.tif + field_grid <- list( + mosaic_dir = mosaic_dir, + mode = "per-field" + ) + + message("\nUsing modular RDS-based approach for weekly statistics...") + + # Load/calculate CURRENT week stats (from tiles or RDS cache) + message("\n1. Loading/calculating CURRENT week statistics (week", current_week, ")...") + current_stats <- load_or_calculate_weekly_stats( + week_num = current_week, + year = current_year, + project_dir = project_dir, + field_boundaries_sf = field_boundaries_sf, + mosaic_dir = field_grid$mosaic_dir, + reports_dir = reports_dir, + report_date = end_date + ) + + message(paste(" ✓ Loaded/calculated stats for", nrow(current_stats), "fields in current week")) + + # Load/calculate PREVIOUS week stats (from RDS cache or tiles) + message("\n2. Loading/calculating PREVIOUS week statistics (week", previous_week, ")...") + + # Calculate report date for previous week (7 days before current) + prev_report_date <- end_date - 7 + + prev_stats <- load_or_calculate_weekly_stats( + week_num = previous_week, + year = previous_year, + project_dir = project_dir, + field_boundaries_sf = field_boundaries_sf, + mosaic_dir = field_grid$mosaic_dir, + reports_dir = reports_dir, + report_date = prev_report_date + ) + + message(paste(" ✓ Loaded/calculated stats for", nrow(prev_stats), "fields in previous week")) + message(paste(" Columns in prev_stats:", paste(names(prev_stats), collapse = ", "))) + message(paste(" CV column non-NA values in prev_stats:", sum(!is.na(prev_stats$CV)))) + + # Apply trend calculations (requires both weeks) + message("\n3. Calculating trend columns...") + current_stats <- calculate_kpi_trends(current_stats, prev_stats, + project_dir = project_dir, + reports_dir = reports_dir, + current_week = current_week, + year = current_year) + + message(paste(" ✓ Added Weekly_ci_change, CV_Trend_Short_Term, Four_week_trend, CV_Trend_Long_Term, nmr_of_weeks_analysed")) + + # Load weekly harvest probabilities from script 31 (if available) + # Note: Script 31 saves to reports/kpis/field_stats/ (not field_level) + message("\n4. Loading harvest probabilities from script 31...") + harvest_prob_dir <- file.path(data_dir, "..", "reports", "kpis", "field_stats") + harvest_prob_file <- file.path(harvest_prob_dir, + sprintf("%s_harvest_imminent_week_%02d_%d.csv", project_dir, current_week, current_year)) + message(paste(" Looking for:", harvest_prob_file)) + + imminent_prob_data <- tryCatch({ + if (file.exists(harvest_prob_file)) { + prob_df <- readr::read_csv(harvest_prob_file, show_col_types = FALSE, + col_types = readr::cols(.default = readr::col_character())) + message(paste(" ✓ Loaded harvest probabilities for", nrow(prob_df), "fields")) + prob_df %>% + select(field, imminent_prob, detected_prob) %>% + rename(Field_id = field, Imminent_prob_actual = imminent_prob, Detected_prob = detected_prob) + } else { + message(paste(" INFO: Harvest probabilities not available (script 31 not run)")) + NULL + } + }, error = function(e) { + message(paste(" WARNING: Could not load harvest probabilities:", e$message)) + NULL + }) + + # ============================================================================ + # CALCULATE GAP FILLING KPI (2σ method from kpi_utils.R) + # ============================================================================ + + message("\nCalculating gap filling scores (2σ method)...") + + # Try single merged mosaic first, then fall back to merging tiles + week_mosaic_file <- file.path(mosaic_dir, sprintf("week_%02d_%d.tif", current_week, current_year)) + + gap_scores_df <- NULL + + if (file.exists(week_mosaic_file)) { + # Single merged mosaic exists - use it directly tryCatch({ - message(sprintf(" Found %d tiles. Processing per-tile (memory efficient)...", length(tile_files))) - - # Process each tile separately and accumulate results - all_tile_results <- list() - - for (i in seq_along(tile_files)) { - tile_file <- tile_files[i] - - # Load tile raster - tile_raster <- terra::rast(tile_file) - tile_ci_band <- tile_raster[[5]] - names(tile_ci_band) <- "CI" - - # Calculate gap scores for fields in this tile - tile_gap_result <- calculate_gap_filling_kpi(tile_ci_band, field_boundaries_sf) - - # Store results (only keep fields with non-NA scores, use field directly to match current_stats) - if (!is.null(tile_gap_result$field_results) && nrow(tile_gap_result$field_results) > 0) { - tile_results_clean <- tile_gap_result$field_results %>% - mutate(Field_id = field) %>% - select(Field_id, gap_score) %>% - filter(!is.na(gap_score)) - - if (nrow(tile_results_clean) > 0) { - all_tile_results[[i]] <- tile_results_clean - } - } - - # Clear memory - rm(tile_raster, tile_ci_band, tile_gap_result) - gc(verbose = FALSE) + current_week_raster <- terra::rast(week_mosaic_file) + # Extract CI band by name (not assumed position) + # Extract CI band (5th band in mosaic) + ci_band_name <- "CI" + if (!(ci_band_name %in% names(current_week_raster))) { + stop(paste("ERROR: CI band not found in mosaic. Available bands:", + paste(names(current_week_raster), collapse = ", "))) } - - # Combine all tile results - if (length(all_tile_results) > 0) { - gap_scores_df <- bind_rows(all_tile_results) - - # If a field appears in multiple tiles, take the maximum gap score - gap_scores_df <- gap_scores_df %>% - group_by(Field_id) %>% - summarise(gap_score = max(gap_score, na.rm = TRUE), .groups = "drop") - - message(paste(" ✓ Calculated gap scores for", nrow(gap_scores_df), "fields across", length(all_tile_results), "tiles")) - message(paste(" Gap score range:", round(min(gap_scores_df$gap_score, na.rm=TRUE), 2), "-", round(max(gap_scores_df$gap_score, na.rm=TRUE), 2), "%")) - } else { - message(" WARNING: No gap scores calculated from any tiles") - gap_scores_df <- NULL + current_ci_band <- current_week_raster[[ci_band_name]] + names(current_ci_band) <- "CI" + if (!(ci_band_name %in% names(current_week_raster))) { + stop(paste("ERROR: CI band not found in mosaic. Available bands:", + paste(names(current_week_raster), collapse = ", "))) } - + current_ci_band <- current_week_raster[[ci_band_name]] + names(current_ci_band) <- "CI" + + message(paste(" Loaded single mosaic:", week_mosaic_file)) + + # Calculate gap scores for all fields + gap_result <- calculate_gap_filling_kpi(current_ci_band, field_boundaries_sf) + + # Extract field-level results (use field column directly to match current_stats Field_id) + gap_scores_df <- gap_result$field_results %>% + mutate(Field_id = field) %>% + select(Field_id, gap_score) + + message(paste(" ✓ Calculated gap scores for", nrow(gap_scores_df), "fields")) + message(paste(" Gap score range:", round(min(gap_scores_df$gap_score, na.rm=TRUE), 2), "-", round(max(gap_scores_df$gap_score, na.rm=TRUE), 2), "%")) + }, error = function(e) { - message(paste(" WARNING: Could not process tiles or calculate gap scores:", e$message)) + message(paste(" WARNING: Could not calculate gap scores from single mosaic:", e$message)) message(" Gap scores will be set to NA") gap_scores_df <- NULL }) - } - } - - # ============================================================================ - # Build final output dataframe with all 22 columns (including Gap_score) - # ============================================================================ - - message("\nBuilding final field analysis output...") - - # Pre-calculate acreages with geometry validation - # This avoids geometry errors during field_analysis construction - acreage_lookup <- tryCatch({ - lookup_df <- field_boundaries_sf %>% - sf::st_drop_geometry() %>% - as.data.frame() %>% - mutate( - geometry_valid = sapply(seq_len(nrow(field_boundaries_sf)), function(idx) { - tryCatch({ - sf::st_is_valid(field_boundaries_sf[idx, ]) - }, error = function(e) FALSE) - }), - area_ha = 0 - ) - # Calculate area for valid geometries - for (idx in which(lookup_df$geometry_valid)) { - tryCatch({ - area_m2 <- as.numeric(sf::st_area(field_boundaries_sf[idx, ])) - lookup_df$area_ha[idx] <- area_m2 / 10000 - }, error = function(e) { - lookup_df$area_ha[idx] <<- NA_real_ - }) - } + } else { + # Single mosaic doesn't exist - check for tiles and process per-tile + message(" Single mosaic not found. Checking for tiles...") - # Convert hectares to acres - lookup_df %>% - mutate(acreage = area_ha / 0.404686) %>% - select(field, acreage) - }, error = function(e) { - message(paste("Warning: Could not calculate acreages from geometries -", e$message)) - data.frame(field = character(0), acreage = numeric(0)) - }) - - field_analysis_df <- current_stats %>% - mutate( - # Column 2: Farm_Section (user fills) - Farm_Section = NA_character_, - # Column 3: Field_name (from GeoJSON - already have Field_id, can look up) - Field_name = Field_id, - # Column 4: Acreage (calculate from geometry) - Acreage = { - acreages_vec <- acreage_lookup$acreage[match(Field_id, acreage_lookup$field)] - if_else(is.na(acreages_vec), 0, acreages_vec) - }, - # Columns 5-6: Already in current_stats (Mean_CI, Weekly_ci_change) - # Column 7: Four_week_trend (from current_stats) - # Column 8: Last_harvest_or_planting_date (from harvest.xlsx - season_start) - Last_harvest_or_planting_date = { - planting_dates$planting_date[match(Field_id, planting_dates$field_id)] - }, - # Column 9: Age_week (calculated from report date and planting date) - Age_week = { - sapply(seq_len(nrow(current_stats)), function(idx) { - planting_dt <- Last_harvest_or_planting_date[idx] - if (is.na(planting_dt)) { - return(NA_real_) - } - round(as.numeric(difftime(end_date, planting_dt, units = "weeks")), 0) - }) - }, - # Column 10: Phase (recalculate based on updated Age_week) - Phase = { - sapply(Age_week, function(age) { - if (is.na(age)) return(NA_character_) - if (age >= 0 & age < 4) return("Germination") - if (age >= 4 & age < 17) return("Tillering") - if (age >= 17 & age < 39) return("Grand Growth") - if (age >= 39) return("Maturation") - NA_character_ - }) - }, - # Column 11: nmr_of_weeks_analysed (already in current_stats from calculate_kpi_trends) - # Column 12: Germination_progress (calculated here from CI values) - # Bin Pct_pixels_CI_gte_2 into 10% intervals: 0-10%, 10-20%, ..., 80-90%, 90-95%, 95-100% - Germination_progress = sapply(Pct_pixels_CI_gte_2, function(pct) { - if (is.na(pct)) return(NA_character_) - if (pct >= 95) return("95-100%") - else if (pct >= 90) return("90-95%") - else if (pct >= 80) return("80-90%") - else if (pct >= 70) return("70-80%") - else if (pct >= 60) return("60-70%") - else if (pct >= 50) return("50-60%") - else if (pct >= 40) return("40-50%") - else if (pct >= 30) return("30-40%") - else if (pct >= 20) return("20-30%") - else if (pct >= 10) return("10-20%") - else return("0-10%") - }), - # Column 13: Imminent_prob (from script 31 or NA if not available) - Imminent_prob = { - if (!is.null(imminent_prob_data)) { - imminent_prob_data$Imminent_prob_actual[match(Field_id, imminent_prob_data$Field_id)] - } else { - rep(NA_real_, nrow(current_stats)) - } - }, - # Column 14: Status_Alert (based on harvest probability + crop health status) - # Priority order: Ready for harvest-check → Strong decline → Harvested/bare → NA - Status_Alert = { - sapply(seq_len(nrow(current_stats)), function(idx) { - imminent_prob <- Imminent_prob[idx] - age_w <- Age_week[idx] - weekly_ci_chg <- Weekly_ci_change[idx] - mean_ci_val <- Mean_CI[idx] + # List all tiles for this week (e.g., week_04_2026_01.tif through week_04_2026_25.tif) + tile_pattern <- sprintf("week_%02d_%d_\\d{2}\\.tif$", current_week, current_year) + tile_files <- list.files(mosaic_dir, pattern = tile_pattern, full.names = TRUE) + + if (length(tile_files) == 0) { + message(sprintf(" WARNING: No tiles found matching pattern: %s in %s", tile_pattern, mosaic_dir)) + message(" Gap scores will be set to NA") + + } else { + tryCatch({ + message(sprintf(" Found %d tiles. Processing per-tile (memory efficient)...", length(tile_files))) + + # Process each tile separately and accumulate results + all_tile_results <- list() + + for (i in seq_along(tile_files)) { + tile_file <- tile_files[i] - # Priority 1: Ready for harvest-check (imminent + mature cane ≥12 months) - if (!is.na(imminent_prob) && imminent_prob > 0.5 && !is.na(age_w) && age_w >= 52) { - return("Ready for harvest-check") - } + # Load tile raster + tile_raster <- terra::rast(tile_file) + + # Extract CI band by name (not assumed position) + ci_band_name <- "CI" + if (!(ci_band_name %in% names(tile_raster))) { + stop(paste("ERROR: CI band not found in tile mosaic. Available bands:", + paste(names(tile_raster), collapse = ", "))) + } + tile_ci_band <- tile_raster[[ci_band_name]] + names(tile_ci_band) <- "CI" - # Priority 2: Strong decline in crop health (drop ≥2 points but still >1.5) - if (!is.na(weekly_ci_chg) && weekly_ci_chg <= -2.0 && !is.na(mean_ci_val) && mean_ci_val > 1.5) { - return("Strong decline in crop health") - } + # Calculate gap scores for fields in this tile + tile_gap_result <- calculate_gap_filling_kpi(tile_ci_band, field_boundaries_sf) - # Priority 3: Harvested/bare (Mean CI < 1.5) - if (!is.na(mean_ci_val) && mean_ci_val < 1.5) { - return("Harvested/bare") - } + # Store results (only keep fields with non-NA scores, use field directly to match current_stats) + if (!is.null(tile_gap_result$field_results) && nrow(tile_gap_result$field_results) > 0) { + tile_results_clean <- tile_gap_result$field_results %>% + mutate(Field_id = field) %>% + select(Field_id, gap_score) %>% + filter(!is.na(gap_score)) + + if (nrow(tile_results_clean) > 0) { + all_tile_results[[i]] <- tile_results_clean + } + } - # Fallback: no alert - NA_character_ - }) - }, - # Columns 15-16: CI-based columns already in current_stats (CI_range, CI_Percentiles) - # Column 17: Already in current_stats (CV) - # Column 18: Already in current_stats (CV_Trend_Short_Term) - # Column 19: CV_Trend_Long_Term (from current_stats - raw slope value) - # Column 19b: CV_Trend_Long_Term_Category (categorical interpretation of slope) - # 3 classes: More uniform (slope < -0.01), Stable uniformity (-0.01 to 0.01), Less uniform (slope > 0.01) - CV_Trend_Long_Term_Category = { - sapply(current_stats$CV_Trend_Long_Term, function(slope) { - if (is.na(slope)) { - return(NA_character_) - } else if (slope < -0.01) { - return("More uniform") - } else if (slope > 0.01) { - return("Less uniform") + # Clear memory + rm(tile_raster, tile_ci_band, tile_gap_result) + gc(verbose = FALSE) + } + + # Combine all tile results + if (length(all_tile_results) > 0) { + gap_scores_df <- bind_rows(all_tile_results) + + # If a field appears in multiple tiles, take the maximum gap score + gap_scores_df <- gap_scores_df %>% + group_by(Field_id) %>% + summarise(gap_score = max(gap_score, na.rm = TRUE), .groups = "drop") + + message(paste(" ✓ Calculated gap scores for", nrow(gap_scores_df), "fields across", length(all_tile_results), "tiles")) + message(paste(" Gap score range:", round(min(gap_scores_df$gap_score, na.rm=TRUE), 2), "-", round(max(gap_scores_df$gap_score, na.rm=TRUE), 2), "%")) } else { - return("Stable uniformity") + message(" WARNING: No gap scores calculated from any tiles") + gap_scores_df <- NULL } + + }, error = function(e) { + message(paste(" WARNING: Could not process tiles or calculate gap scores:", e$message)) + message(" Gap scores will be set to NA") + gap_scores_df <- NULL }) - }, - # Columns 20-21: Already in current_stats (Cloud_pct_clear, Cloud_category) - # Bin Cloud_pct_clear into 10% intervals: 0-10%, 10-20%, ..., 80-90%, 90-95%, 95-100% - Cloud_pct_clear = sapply(Cloud_pct_clear, function(pct) { - if (is.na(pct)) return(NA_character_) - if (pct >= 95) return("95-100%") - else if (pct >= 90) return("90-95%") - else if (pct >= 80) return("80-90%") - else if (pct >= 70) return("70-80%") - else if (pct >= 60) return("60-70%") - else if (pct >= 50) return("50-60%") - else if (pct >= 40) return("40-50%") - else if (pct >= 30) return("30-40%") - else if (pct >= 20) return("20-30%") - else if (pct >= 10) return("10-20%") - else return("0-10%") - }), - # Column 22: Gap_score (2σ below median - from kpi_utils.R) - Gap_score = { - if (!is.null(gap_scores_df) && nrow(gap_scores_df) > 0) { - # Debug: Print first few gap scores - message(sprintf(" Joining %d gap scores to field_analysis (first 3: %s)", - nrow(gap_scores_df), - paste(head(gap_scores_df$gap_score, 3), collapse=", "))) - message(sprintf(" First 3 Field_ids in gap_scores_df: %s", - paste(head(gap_scores_df$Field_id, 3), collapse=", "))) - message(sprintf(" First 3 Field_ids in current_stats: %s", - paste(head(current_stats$Field_id, 3), collapse=", "))) - - gap_scores_df$gap_score[match(current_stats$Field_id, gap_scores_df$Field_id)] - } else { - rep(NA_real_, nrow(current_stats)) - } } - ) %>% - select( - all_of(c("Field_id", "Farm_Section", "Field_name", "Acreage", "Status_Alert", - "Last_harvest_or_planting_date", "Age_week", "Phase", - "Germination_progress", - "Mean_CI", "Weekly_ci_change", "Four_week_trend", "CI_range", "CI_Percentiles", - "CV", "CV_Trend_Short_Term", "CV_Trend_Long_Term", "CV_Trend_Long_Term_Category", - "Imminent_prob", "Cloud_pct_clear", "Cloud_category", "Gap_score")) + } + + # ============================================================================ + # Build final output dataframe with all 22 columns (including Gap_score) + # ============================================================================ + + message("\nBuilding final field analysis output...") + + # Pre-calculate acreages with geometry validation + # This avoids geometry errors during field_analysis construction + acreage_lookup <- tryCatch({ + lookup_df <- field_boundaries_sf %>% + sf::st_drop_geometry() %>% + as.data.frame() %>% + mutate( + geometry_valid = sapply(seq_len(nrow(field_boundaries_sf)), function(idx) { + tryCatch({ + sf::st_is_valid(field_boundaries_sf[idx, ]) + }, error = function(e) FALSE) + }), + area_ha = 0 + ) + + # Calculate area for valid geometries + for (idx in which(lookup_df$geometry_valid)) { + tryCatch({ + area_m2 <- as.numeric(sf::st_area(field_boundaries_sf[idx, ])) + lookup_df$area_ha[idx] <- area_m2 / 10000 + }, error = function(e) { + lookup_df$area_ha[idx] <<- NA_real_ + }) + } + + # Convert hectares to acres + lookup_df %>% + mutate(acreage = area_ha / 0.404686) %>% + select(field, acreage) + }, error = function(e) { + message(paste("Warning: Could not calculate acreages from geometries -", e$message)) + data.frame(field = character(0), acreage = numeric(0)) + }) + + field_analysis_df <- current_stats %>% + mutate( + # Column 2: Farm_Section (user fills) + Farm_Section = NA_character_, + # Column 3: Field_name (from GeoJSON - already have Field_id, can look up) + Field_name = Field_id, + # Column 4: Acreage (calculate from geometry) + Acreage = { + acreages_vec <- acreage_lookup$acreage[match(Field_id, acreage_lookup$field)] + if_else(is.na(acreages_vec), 0, acreages_vec) + }, + # Columns 5-6: Already in current_stats (Mean_CI, Weekly_ci_change) + # Column 7: Four_week_trend (from current_stats) + # Column 8: Last_harvest_or_planting_date (from harvest.xlsx - season_start) + Last_harvest_or_planting_date = { + planting_dates$planting_date[match(Field_id, planting_dates$field_id)] + }, + # Column 9: Age_week (calculated from report date and planting date) + Age_week = { + sapply(seq_len(nrow(current_stats)), function(idx) { + planting_dt <- Last_harvest_or_planting_date[idx] + if (is.na(planting_dt)) { + return(NA_real_) + } + round(as.numeric(difftime(end_date, planting_dt, units = "weeks")), 0) + }) + }, + # Column 10: Phase (recalculate based on updated Age_week) + Phase = { + sapply(Age_week, function(age) { + if (is.na(age)) return(NA_character_) + if (age >= 0 & age < 4) return("Germination") + if (age >= 4 & age < 17) return("Tillering") + if (age >= 17 & age < 39) return("Grand Growth") + if (age >= 39) return("Maturation") + NA_character_ + }) + }, + # Column 11: nmr_of_weeks_analysed (already in current_stats from calculate_kpi_trends) + # Column 12: Germination_progress (calculated here from CI values) + # Bin Pct_pixels_CI_gte_2 into 10% intervals: 0-10%, 10-20%, ..., 80-90%, 90-95%, 95-100% + Germination_progress = sapply(Pct_pixels_CI_gte_2, function(pct) { + if (is.na(pct)) return(NA_character_) + if (pct >= 95) return("95-100%") + else if (pct >= 90) return("90-95%") + else if (pct >= 80) return("80-90%") + else if (pct >= 70) return("70-80%") + else if (pct >= 60) return("60-70%") + else if (pct >= 50) return("50-60%") + else if (pct >= 40) return("40-50%") + else if (pct >= 30) return("30-40%") + else if (pct >= 20) return("20-30%") + else if (pct >= 10) return("10-20%") + else return("0-10%") + }), + # Column 13: Imminent_prob (from script 31 or NA if not available) + Imminent_prob = { + if (!is.null(imminent_prob_data)) { + imminent_prob_data$Imminent_prob_actual[match(Field_id, imminent_prob_data$Field_id)] + } else { + rep(NA_real_, nrow(current_stats)) + } + }, + # Column 14: Status_Alert (based on harvest probability + crop health status) + # Priority order: Ready for harvest-check → Strong decline → Harvested/bare → NA + Status_Alert = { + sapply(seq_len(nrow(current_stats)), function(idx) { + imminent_prob <- Imminent_prob[idx] + age_w <- Age_week[idx] + weekly_ci_chg <- Weekly_ci_change[idx] + mean_ci_val <- Mean_CI[idx] + + # Priority 1: Ready for harvest-check (imminent + mature cane ≥12 months) + if (!is.na(imminent_prob) && imminent_prob > 0.5 && !is.na(age_w) && age_w >= 52) { + return("Ready for harvest-check") + } + + # Priority 2: Strong decline in crop health (drop ≥2 points but still >1.5) + if (!is.na(weekly_ci_chg) && weekly_ci_chg <= -2.0 && !is.na(mean_ci_val) && mean_ci_val > 1.5) { + return("Strong decline in crop health") + } + + # Priority 3: Harvested/bare (Mean CI < 1.5) + if (!is.na(mean_ci_val) && mean_ci_val < 1.5) { + return("Harvested/bare") + } + + # Fallback: no alert + NA_character_ + }) + }, + # Columns 15-16: CI-based columns already in current_stats (CI_range, CI_Percentiles) + # Column 17: Already in current_stats (CV) + # Column 18: Already in current_stats (CV_Trend_Short_Term) + # Column 19: CV_Trend_Long_Term (from current_stats - raw slope value) + # Column 19b: CV_Trend_Long_Term_Category (categorical interpretation of slope) + # 3 classes: More uniform (slope < -0.01), Stable uniformity (-0.01 to 0.01), Less uniform (slope > 0.01) + CV_Trend_Long_Term_Category = { + sapply(current_stats$CV_Trend_Long_Term, function(slope) { + if (is.na(slope)) { + return(NA_character_) + } else if (slope < -0.01) { + return("More uniform") + } else if (slope > 0.01) { + return("Less uniform") + } else { + return("Stable uniformity") + } + }) + }, + # Columns 20-21: Already in current_stats (Cloud_pct_clear, Cloud_category) + # Bin Cloud_pct_clear into 10% intervals: 0-10%, 10-20%, ..., 80-90%, 90-95%, 95-100% + Cloud_pct_clear = sapply(Cloud_pct_clear, function(pct) { + if (is.na(pct)) return(NA_character_) + if (pct >= 95) return("95-100%") + else if (pct >= 90) return("90-95%") + else if (pct >= 80) return("80-90%") + else if (pct >= 70) return("70-80%") + else if (pct >= 60) return("60-70%") + else if (pct >= 50) return("50-60%") + else if (pct >= 40) return("40-50%") + else if (pct >= 30) return("30-40%") + else if (pct >= 20) return("20-30%") + else if (pct >= 10) return("10-20%") + else return("0-10%") + }), + # Column 22: Gap_score (2σ below median - from kpi_utils.R) + Gap_score = { + if (!is.null(gap_scores_df) && nrow(gap_scores_df) > 0) { + # Debug: Print first few gap scores + message(sprintf(" Joining %d gap scores to field_analysis (first 3: %s)", + nrow(gap_scores_df), + paste(head(gap_scores_df$gap_score, 3), collapse=", "))) + message(sprintf(" First 3 Field_ids in gap_scores_df: %s", + paste(head(gap_scores_df$Field_id, 3), collapse=", "))) + message(sprintf(" First 3 Field_ids in current_stats: %s", + paste(head(current_stats$Field_id, 3), collapse=", "))) + + gap_scores_df$gap_score[match(current_stats$Field_id, gap_scores_df$Field_id)] + } else { + rep(NA_real_, nrow(current_stats)) + } + } + ) %>% + select( + all_of(c("Field_id", "Farm_Section", "Field_name", "Acreage", "Status_Alert", + "Last_harvest_or_planting_date", "Age_week", "Phase", + "Germination_progress", + "Mean_CI", "Weekly_ci_change", "Four_week_trend", "CI_range", "CI_Percentiles", + "CV", "CV_Trend_Short_Term", "CV_Trend_Long_Term", "CV_Trend_Long_Term_Category", + "Imminent_prob", "Cloud_pct_clear", "Cloud_category", "Gap_score")) + ) + + message(paste("✓ Built final output with", nrow(field_analysis_df), "fields and 22 columns (including Gap_score)")) + + export_paths <- export_field_analysis_excel( + field_analysis_df, + NULL, + project_dir, + current_week, + current_year, + reports_dir ) - message(paste("✓ Built final output with", nrow(field_analysis_df), "fields and 22 columns (including Gap_score)")) + cat("\n--- Per-field Results (first 10) ---\n") + available_cols <- c("Field_id", "Acreage", "Age_week", "Mean_CI", "Four_week_trend", "Status_Alert", "Cloud_category") + available_cols <- available_cols[available_cols %in% names(field_analysis_df)] + if (length(available_cols) > 0) { + print(head(field_analysis_df[, available_cols], 10)) + } - export_paths <- export_field_analysis_excel( - field_analysis_df, - NULL, - project_dir, - current_week, - current_iso_year, - reports_dir - ) + # ========== FARM-LEVEL KPI AGGREGATION ========== + # Aggregate the per-field analysis into farm-level summary statistics - cat("\n--- Per-field Results (first 10) ---\n") - available_cols <- c("Field_id", "Acreage", "Age_week", "Mean_CI", "Four_week_trend", "Status_Alert", "Cloud_category") - available_cols <- available_cols[available_cols %in% names(field_analysis_df)] - if (length(available_cols) > 0) { - print(head(field_analysis_df[, available_cols], 10)) - } + cat("\n=== CALCULATING FARM-LEVEL KPI SUMMARY ===\n") - # ========== FARM-LEVEL KPI AGGREGATION ========== - # Aggregate the per-field analysis into farm-level summary statistics + # Filter to only fields that have actual data (non-NA CI and valid acreage) + field_data <- field_analysis_df %>% + filter(!is.na(Mean_CI) & !is.na(Acreage)) %>% + filter(Acreage > 0) - cat("\n=== CALCULATING FARM-LEVEL KPI SUMMARY ===\n") - - # Filter to only fields that have actual data (non-NA CI and valid acreage) - field_data <- field_analysis_df %>% - filter(!is.na(Mean_CI) & !is.na(Acreage)) %>% - filter(Acreage > 0) - - if (nrow(field_data) > 0) { - if (nrow(field_data) > 0) { - # Create summary statistics - farm_summary <- list() + + if (nrow(field_data) > 0) { + # Create summary statistics + farm_summary <- list() - # 1. PHASE DISTRIBUTION - phase_dist <- field_data %>% - group_by(Phase) %>% - summarise( - num_fields = n(), - acreage = sum(Acreage, na.rm = TRUE), - .groups = 'drop' - ) %>% - rename(Category = Phase) + # 1. PHASE DISTRIBUTION + phase_dist <- field_data %>% + group_by(Phase) %>% + summarise( + num_fields = n(), + acreage = sum(Acreage, na.rm = TRUE), + .groups = 'drop' + ) %>% + rename(Category = Phase) - farm_summary$phase_distribution <- phase_dist + farm_summary$phase_distribution <- phase_dist - # 2. STATUS ALERT DISTRIBUTION - status_dist <- field_data %>% - group_by(Status_Alert) %>% - summarise( - num_fields = n(), - acreage = sum(Acreage, na.rm = TRUE), - .groups = 'drop' - ) %>% - rename(Category = Status_Alert) + # 2. STATUS ALERT DISTRIBUTION + status_dist <- field_data %>% + group_by(Status_Alert) %>% + summarise( + num_fields = n(), + acreage = sum(Acreage, na.rm = TRUE), + .groups = 'drop' + ) %>% + rename(Category = Status_Alert) - farm_summary$status_distribution <- status_dist + farm_summary$status_distribution <- status_dist - # 3. CLOUD COVERAGE DISTRIBUTION - cloud_dist <- field_data %>% - group_by(Cloud_category) %>% - summarise( - num_fields = n(), - acreage = sum(Acreage, na.rm = TRUE), - .groups = 'drop' - ) %>% - rename(Category = Cloud_category) + # 3. CLOUD COVERAGE DISTRIBUTION + cloud_dist <- field_data %>% + group_by(Cloud_category) %>% + summarise( + num_fields = n(), + acreage = sum(Acreage, na.rm = TRUE), + .groups = 'drop' + ) %>% + rename(Category = Cloud_category) - farm_summary$cloud_distribution <- cloud_dist + farm_summary$cloud_distribution <- cloud_dist - # 4. OVERALL STATISTICS - farm_summary$overall_stats <- data.frame( - total_fields = nrow(field_data), - total_acreage = sum(field_data$Acreage, na.rm = TRUE), - mean_ci = round(mean(field_data$Mean_CI, na.rm = TRUE), 2), - median_ci = round(median(field_data$Mean_CI, na.rm = TRUE), 2), - mean_cv = round(mean(field_data$CV, na.rm = TRUE), 4), - week = current_week, - year = current_iso_year, - date = as.character(end_date) - ) + # 4. OVERALL STATISTICS + farm_summary$overall_stats <- data.frame( + total_fields = nrow(field_data), + total_acreage = sum(field_data$Acreage, na.rm = TRUE), + mean_ci = round(mean(field_data$Mean_CI, na.rm = TRUE), 2), + median_ci = round(median(field_data$Mean_CI, na.rm = TRUE), 2), + mean_cv = round(mean(field_data$CV, na.rm = TRUE), 4), + week = current_week, + year = current_year, + date = as.character(end_date) + ) - # Print summaries - cat("\n--- PHASE DISTRIBUTION ---\n") - print(phase_dist) + # Print summaries + cat("\n--- PHASE DISTRIBUTION ---\n") + print(phase_dist) - cat("\n--- STATUS TRIGGER DISTRIBUTION ---\n") - print(status_dist) + cat("\n--- STATUS TRIGGER DISTRIBUTION ---\n") + print(status_dist) - cat("\n--- CLOUD COVERAGE DISTRIBUTION ---\n") - print(cloud_dist) + cat("\n--- CLOUD COVERAGE DISTRIBUTION ---\n") + print(cloud_dist) - cat("\n--- OVERALL FARM STATISTICS ---\n") - print(farm_summary$overall_stats) + cat("\n--- OVERALL FARM STATISTICS ---\n") + print(farm_summary$overall_stats) - farm_kpi_results <- farm_summary + farm_kpi_results <- farm_summary + } else { + farm_kpi_results <- NULL + } } else { farm_kpi_results <- NULL } - } else { - farm_kpi_results <- NULL - } - # ========== FINAL SUMMARY ========== + # ========== FINAL SUMMARY ========== - cat("\n", strrep("=", 70), "\n") - cat("80_CALCULATE_KPIs.R - COMPLETION SUMMARY\n") - cat(strrep("=", 70), "\n") - cat("Per-field analysis fields analyzed:", nrow(field_analysis_df), "\n") - cat("Excel export:", export_paths$excel, "\n") - cat("RDS export:", export_paths$rds, "\n") - cat("CSV export:", export_paths$csv, "\n") + cat("\n", strrep("=", 70), "\n") + cat("80_CALCULATE_KPIs.R - COMPLETION SUMMARY\n") + cat(strrep("=", 70), "\n") + cat("Per-field analysis fields analyzed:", nrow(field_analysis_df), "\n") + cat("Excel export:", export_paths$excel, "\n") + cat("RDS export:", export_paths$rds, "\n") + cat("CSV export:", export_paths$csv, "\n") - if (!is.null(farm_kpi_results)) { - cat("\nFarm-level KPIs: CALCULATED\n") - } else { - cat("\nFarm-level KPIs: SKIPPED (no valid tile data extracted)\n") - } + if (!is.null(farm_kpi_results)) { + cat("\nFarm-level KPIs: CALCULATED\n") + } else { + cat("\nFarm-level KPIs: SKIPPED (no valid tile data extracted)\n") + } - cat("\n✓ Consolidated KPI calculation complete!\n") - cat(" - Per-field data exported\n") - cat(" - Farm-level KPIs calculated\n") - cat(" - All outputs in:", reports_dir, "\n\n") + cat("\n✓ Consolidated KPI calculation complete!\n") + cat(" - Per-field data exported\n") + cat(" - Farm-level KPIs calculated\n") + cat(" - All outputs in:", reports_dir, "\n\n") - } else { - # Unknown client type - log warning and exit - warning(sprintf("Unknown client type: %s - no workflow matched", client_type)) + } else { + # Unknown client type - log warning and exit + warning(sprintf("Unknown client type: %s - no workflow matched", client_type)) cat("\n⚠️ Warning: Client type '", client_type, "' does not match any known workflow\n", sep = "") cat("Expected: 'agronomic_support' (aura) or 'cane_supply' (angata, etc.)\n") cat("Check CLIENT_TYPE_MAP in parameters_project.R\n\n") diff --git a/r_app/80_utils_agronomic_support.R b/r_app/80_utils_agronomic_support.R index c1d710a..992909c 100644 --- a/r_app/80_utils_agronomic_support.R +++ b/r_app/80_utils_agronomic_support.R @@ -353,71 +353,75 @@ calculate_weed_presence_kpi <- function(ci_pixels_by_field) { return(result) } -#' KPI 6: Calculate gap filling quality (data interpolation success) -#' -#' Measures how well cloud/missing data was interpolated during growth model -#' -#' @param ci_rds_path Path to combined CI RDS file (before/after interpolation) -#' -#' @return Data frame with gap-filling quality metrics -calculate_gap_filling_kpi <- function(ci_rds_path) { - # If ci_rds_path is NULL or not a valid path, return placeholder - if (is.null(ci_rds_path) || !is.character(ci_rds_path) || length(ci_rds_path) == 0) { - return(NULL) +#' Calculate Gap Filling Score KPI (placeholder) +#' @param ci_raster Current week CI raster +#' @param field_boundaries Field boundaries +#' @return List with summary data frame and field-level results data frame +calculate_gap_filling_kpi <- function(ci_raster, field_boundaries) { + safe_log("Calculating Gap Filling Score KPI (placeholder)") + + # Handle both sf and SpatVector inputs + if (!inherits(field_boundaries, "SpatVector")) { + field_boundaries_vect <- terra::vect(field_boundaries) + } else { + field_boundaries_vect <- field_boundaries } - - # If ci_rds_path is a directory, find the cumulative CI file - if (dir.exists(ci_rds_path)) { - ci_files <- list.files(ci_rds_path, pattern = "^All_pivots.*\\.rds$", full.names = TRUE) - if (length(ci_files) == 0) { - return(NULL) + + field_results <- data.frame() + + for (i in seq_len(nrow(field_boundaries))) { + field_name <- field_boundaries$field[i] + sub_field_name <- field_boundaries$sub_field[i] + field_vect <- field_boundaries_vect[i] + + # Extract CI values using helper function + ci_values <- extract_ci_values(ci_raster, field_vect) + valid_values <- ci_values[!is.na(ci_values) & is.finite(ci_values)] + + if (length(valid_values) > 1) { + # Gap score using 2σ below median to detect outliers + median_ci <- median(valid_values) + sd_ci <- sd(valid_values) + outlier_threshold <- median_ci - (2 * sd_ci) + low_ci_pixels <- sum(valid_values < outlier_threshold) + total_pixels <- length(valid_values) + gap_score <- (low_ci_pixels / total_pixels) * 100 + + # Classify gap severity + gap_level <- dplyr::case_when( + gap_score < 10 ~ "Minimal", + gap_score < 25 ~ "Moderate", + TRUE ~ "Significant" + ) + + field_results <- rbind(field_results, data.frame( + field = field_name, + sub_field = sub_field_name, + gap_level = gap_level, + gap_score = gap_score, + mean_ci = mean(valid_values), + outlier_threshold = outlier_threshold + )) + } else { + # Not enough valid data, fill with NA row + field_results <- rbind(field_results, data.frame( + field = field_name, + sub_field = sub_field_name, + gap_level = NA_character_, + gap_score = NA_real_, + mean_ci = NA_real_, + outlier_threshold = NA_real_ + )) } - ci_rds_path <- ci_files[1] } - - if (!file.exists(ci_rds_path)) { - return(NULL) - } - - tryCatch({ - ci_data <- readRDS(ci_rds_path) - - # ci_data should be a wide matrix: fields × weeks - # NA values = missing data before interpolation - # (Gap filling is done during growth model stage) - - result <- data.frame( - field_idx = seq_len(nrow(ci_data)), - na_percent_pre_interpolation = NA_real_, - na_percent_post_interpolation = NA_real_, - gap_filling_success = NA_character_, - stringsAsFactors = FALSE - ) - - for (field_idx in seq_len(nrow(ci_data))) { - na_count <- sum(is.na(ci_data[field_idx, ])) - na_pct <- na_count / ncol(ci_data) * 100 - - if (na_pct == 0) { - result$gap_filling_success[field_idx] <- "No gaps (100% data)" - } else if (na_pct < 10) { - result$gap_filling_success[field_idx] <- "Excellent" - } else if (na_pct < 25) { - result$gap_filling_success[field_idx] <- "Good" - } else if (na_pct < 40) { - result$gap_filling_success[field_idx] <- "Fair" - } else { - result$gap_filling_success[field_idx] <- "Poor" - } - - result$na_percent_pre_interpolation[field_idx] <- round(na_pct, 2) - } - - return(result) - }, error = function(e) { - message(paste("Error calculating gap filling KPI:", e$message)) - return(NULL) - }) + + # Summarize results + gap_summary <- field_results %>% + dplyr::group_by(gap_level) %>% + dplyr::summarise(field_count = n(), .groups = 'drop') %>% + dplyr::mutate(percent = round((field_count / sum(field_count)) * 100, 1)) + + return(list(summary = gap_summary, field_results = field_results)) } # ============================================================================ diff --git a/r_app/80_utils_common.R b/r_app/80_utils_common.R index 3c85f8a..9c9ff91 100644 --- a/r_app/80_utils_common.R +++ b/r_app/80_utils_common.R @@ -822,7 +822,8 @@ load_historical_field_data <- function(project_dir, current_week, current_year, if (file.exists(csv_path)) { tryCatch({ - data <- readr::read_csv(csv_path, show_col_types = FALSE) + data <- readr::read_csv(csv_path, show_col_types = FALSE, + col_types = readr::cols(.default = readr::col_character())) historical_data[[lookback + 1]] <- list( week = target_week, year = target_year, @@ -878,7 +879,8 @@ calculate_kpi_trends <- function(current_stats, prev_stats = NULL, analysis_files <- list.files(analysis_dir, pattern = "_field_analysis_week.*\\.csv$", full.names = TRUE) if (length(analysis_files) > 0) { recent_file <- analysis_files[which.max(file.info(analysis_files)$mtime)] - prev_field_analysis <- readr::read_csv(recent_file, show_col_types = FALSE, + prev_field_analysis <- readr::read_csv(recent_file, show_col_types = FALSE, + col_types = readr::cols(.default = readr::col_character()), col_select = c(Field_id, nmr_of_weeks_analysed, Phase)) } } @@ -1143,7 +1145,7 @@ calculate_week_numbers <- function(report_date = Sys.Date()) { return(list( current_week = current_week, previous_week = previous_week, - year = current_year, + current_year = current_year, previous_year = previous_year )) } diff --git a/r_app/90_CI_report_with_kpis_simple.Rmd b/r_app/90_CI_report_with_kpis_simple.Rmd index 0243b9f..c3a005c 100644 --- a/r_app/90_CI_report_with_kpis_simple.Rmd +++ b/r_app/90_CI_report_with_kpis_simple.Rmd @@ -704,7 +704,7 @@ The following table provides a comprehensive overview of all monitored fields wi ```{r detailed_field_table, echo=FALSE, results='asis'} # Load CI quadrant data to get field ages -CI_quadrant <- readRDS(here::here(paths$cumulative_ci_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds")) +#CI_quadrant <- readRDS(here::here(paths$cumulative_ci_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds")) # Identify the current season for each field based on report_date # The current season is the one where the report_date falls within or shortly after the season diff --git a/r_app/DEBUG_remove_date_tiffs.R b/r_app/DEBUG_remove_date_tiffs.R new file mode 100644 index 0000000..949d5a3 --- /dev/null +++ b/r_app/DEBUG_remove_date_tiffs.R @@ -0,0 +1,272 @@ +# ============================================================================== +# DEBUG_REMOVE_DATE_TIFFS.R +# ============================================================================== +# PURPOSE: +# Remove all TIFFs of a specific date from multiple storage folders. +# Useful for debugging/re-running parts of the pipeline without full re-download. +# +# USAGE: +# Rscript DEBUG_remove_date_tiffs.R [project] [date] [--dry-run] [--skip-merged] [--skip-field-tiles] [--skip-field-tiles-ci] [--skip-daily-vals] +# +# EXAMPLES: +# # Remove 2026-02-08 from all folders (WITH CONFIRMATION) +# Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 +# +# # Remove from all folders without confirmation +# Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --no-confirm +# +# # Dry run - show what WOULD be deleted without deleting +# Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --dry-run +# +# # Remove only from merged_tif and field_tiles, skip CI folders +# Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --skip-field-tiles-ci --skip-daily-vals +# +# # Remove from field_tiles_CI only +# Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --skip-merged --skip-field-tiles --skip-daily-vals +# +# ============================================================================== + +# ============================================================================== +# CONFIGURATION - TOGGLE WHICH FOLDERS TO DELETE FROM (DEFAULT: ALL) +# ============================================================================== + +# Set these to FALSE to skip deletion from that folder +DELETE_FROM_MERGED_TIF <- TRUE +DELETE_FROM_FIELD_TILES <- TRUE +DELETE_FROM_FIELD_TILES_CI <- TRUE +DELETE_FROM_DAILY_VALS <- TRUE + +# Safety settings +DRY_RUN <- FALSE # Set to TRUE to preview deletions without actually deleting +REQUIRE_CONFIRMATION <- TRUE # Set to FALSE to delete without asking + +# ============================================================================== +# MAIN FUNCTION +# ============================================================================== + +main <- function() { + # Parse command-line arguments + args <- commandArgs(trailingOnly = TRUE) + + # Validate minimum arguments + if (length(args) < 2) { + cat("\n[ERROR] Missing arguments\n") + cat("Usage: Rscript DEBUG_remove_date_tiffs.R [project] [date] [options]\n\n") + cat("Examples:\n") + cat(" Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08\n") + cat(" Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --dry-run\n") + cat(" Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --skip-field-tiles-ci\n\n") + cat("Options:\n") + cat(" --dry-run Preview deletions without actually deleting\n") + cat(" --no-confirm Delete without confirmation\n") + cat(" --skip-merged Skip merged_tif folder\n") + cat(" --skip-field-tiles Skip field_tiles folder\n") + cat(" --skip-field-tiles-ci Skip field_tiles_CI folder\n") + cat(" --skip-daily-vals Skip daily_vals folder\n\n") + quit(status = 1) + } + + # Parse positional arguments + project <- args[1] + date_str <- args[2] + + # Parse optional flags + if (length(args) >= 3) { + for (i in 3:length(args)) { + arg <- args[i] + + # Skip NA or empty arguments + if (is.na(arg) || nchar(arg) == 0) { + next + } + + if (arg == "--dry-run") { + DRY_RUN <<- TRUE + } else if (arg == "--no-confirm") { + REQUIRE_CONFIRMATION <<- FALSE + } else if (arg == "--skip-merged") { + DELETE_FROM_MERGED_TIF <<- FALSE + } else if (arg == "--skip-field-tiles") { + DELETE_FROM_FIELD_TILES <<- FALSE + } else if (arg == "--skip-field-tiles-ci") { + DELETE_FROM_FIELD_TILES_CI <<- FALSE + } else if (arg == "--skip-daily-vals") { + DELETE_FROM_DAILY_VALS <<- FALSE + } + } + } + + # Validate date format + date_obj <- tryCatch( + as.Date(date_str, format = "%Y-%m-%d"), + error = function(e) NULL + ) + + if (is.na(date_obj)) { + cat(sprintf("[ERROR] Invalid date format: %s (expected YYYY-MM-DD)\n", date_str)) + quit(status = 1) + } + + # =========================================================================== + # BUILD LIST OF FOLDERS & FILES TO DELETE + # =========================================================================== + + base_path <- file.path("laravel_app", "storage", "app", project) + + files_to_delete <- list() + + # FOLDER 1: merged_tif/{DATE}.tif + if (DELETE_FROM_MERGED_TIF) { + merged_tif_file <- file.path(base_path, "merged_tif", paste0(date_str, ".tif")) + if (file.exists(merged_tif_file)) { + files_to_delete[["merged_tif"]] <- merged_tif_file + } + } + + # FOLDER 2: field_tiles/{FIELD}/{DATE}.tif (per-field structure) + if (DELETE_FROM_FIELD_TILES) { + field_tiles_dir <- file.path(base_path, "field_tiles") + if (dir.exists(field_tiles_dir)) { + field_dirs <- list.dirs(field_tiles_dir, full.names = TRUE, recursive = FALSE) + for (field_dir in field_dirs) { + tif_file <- file.path(field_dir, paste0(date_str, ".tif")) + if (file.exists(tif_file)) { + folder_name <- basename(field_dir) + key <- paste0("field_tiles/", folder_name) + files_to_delete[[key]] <- tif_file + } + } + } + } + + # FOLDER 3: field_tiles_CI/{FIELD}/{DATE}.tif (per-field structure) + if (DELETE_FROM_FIELD_TILES_CI) { + field_tiles_ci_dir <- file.path(base_path, "field_tiles_CI") + if (dir.exists(field_tiles_ci_dir)) { + field_dirs <- list.dirs(field_tiles_ci_dir, full.names = TRUE, recursive = FALSE) + for (field_dir in field_dirs) { + tif_file <- file.path(field_dir, paste0(date_str, ".tif")) + if (file.exists(tif_file)) { + folder_name <- basename(field_dir) + key <- paste0("field_tiles_CI/", folder_name) + files_to_delete[[key]] <- tif_file + } + } + } + } + + # FOLDER 4: Data/extracted_ci/daily_vals/{SUBDIR}/{DATE}.rds (per-subdirectory structure) + if (DELETE_FROM_DAILY_VALS) { + daily_vals_dir <- file.path(base_path, "Data", "extracted_ci", "daily_vals") + if (dir.exists(daily_vals_dir)) { + subdirs <- list.dirs(daily_vals_dir, full.names = TRUE, recursive = FALSE) + for (subdir in subdirs) { + rds_file <- file.path(subdir, paste0(date_str, ".rds")) + if (file.exists(rds_file)) { + subdir_name <- basename(subdir) + key <- paste0("daily_vals/", subdir_name) + files_to_delete[[key]] <- rds_file + } + } + } + } + + # =========================================================================== + # SUMMARY & CONFIRMATION + # =========================================================================== + + cat("\n") + cat(strrep("=", 70), "\n") + cat("DELETE DATE TIFFS - SUMMARY\n") + cat(strrep("=", 70), "\n") + cat(sprintf("Project: %s\n", project)) + cat(sprintf("Date: %s\n", date_str)) + cat(sprintf("Dry run: %s\n", if (DRY_RUN) "YES" else "NO")) + cat(sprintf("Files to delete: %d\n", length(files_to_delete))) + cat("\n") + + if (length(files_to_delete) == 0) { + cat("[INFO] No files found to delete\n") + cat(strrep("=", 70), "\n\n") + quit(status = 0) + } + + # Count files by folder type for compact summary + folder_counts <- table(sapply(names(files_to_delete), function(key) strsplit(key, "/")[[1]][1])) + cat("Files to delete by folder:\n") + for (folder in names(folder_counts)) { + cat(sprintf(" %s: %d file%s\n", folder, folder_counts[folder], if (folder_counts[folder] != 1) "s" else "")) + } + cat(sprintf(" Total: %d file%s\n", length(files_to_delete), if (length(files_to_delete) != 1) "s" else "")) + cat("\n") + + # Ask for confirmation (unless --no-confirm flag was used) + if (REQUIRE_CONFIRMATION && !DRY_RUN) { + # Check if running in interactive mode + if (!interactive()) { + cat("\n[ERROR] Non-interactive mode detected (running via Rscript)\n") + cat("Cannot prompt for confirmation. Use --no-confirm flag to proceed:\n") + cat(" Rscript DEBUG_remove_date_tiffs.R angata 2026-02-08 --no-confirm\n\n") + cat(strrep("=", 70), "\n\n") + quit(status = 1) + } + + cat("⚠️ This will PERMANENTLY DELETE the above files!\n") + cat("Use --no-confirm flag to skip this prompt\n") + + # Use readline() for interactive input (only works in interactive R/RStudio) + response <- readline(prompt = "Type 'yes' to confirm, or anything else to cancel: ") + + if (tolower(response) != "yes") { + cat("[CANCELLED] No files deleted\n") + cat(strrep("=", 70), "\n\n") + quit(status = 0) + } + } + + # =========================================================================== + # DELETE OR DRY-RUN + # =========================================================================== + + deleted_count <- 0 + error_count <- 0 + + for (i in seq_along(files_to_delete)) { + folder_key <- names(files_to_delete)[i] + file_path <- files_to_delete[[i]] + + if (!DRY_RUN) { + tryCatch({ + file.remove(file_path) + deleted_count <- deleted_count + 1 + }, error = function(e) { + error_count <<- error_count + 1 + }) + } + } + + # =========================================================================== + # FINAL SUMMARY + # =========================================================================== + + cat("\n") + if (DRY_RUN) { + cat(sprintf("[DRY RUN] Would have deleted %d files\n", length(files_to_delete))) + } else { + cat(sprintf("Deleted: %d files\n", deleted_count)) + if (error_count > 0) { + cat(sprintf("Errors: %d files\n", error_count)) + } + } + cat(strrep("=", 70), "\n\n") + + quit(status = 0) +} + +# ============================================================================== +# EXECUTE +# ============================================================================== + +if (sys.nframe() == 0) { + main() +} diff --git a/r_app/FIX_INDENTATION.R b/r_app/FIX_INDENTATION.R new file mode 100644 index 0000000..6d7a47c --- /dev/null +++ b/r_app/FIX_INDENTATION.R @@ -0,0 +1,20 @@ +# Fix indentation for lines 408-1022 in 80_calculate_kpis.R +# These lines should be inside the else-if block at the CANE_SUPPLY_WORKFLOW level + +file_path <- "r_app/80_calculate_kpis.R" +lines <- readLines(file_path) + +# Lines 408-1021 (0-indexed: 407-1020) need 2 more spaces of indentation +for (i in 408:1021) { + if (i <= length(lines)) { + line <- lines[i] + # Skip empty or whitespace-only lines + if (nchar(trimws(line)) > 0) { + # Add 2 spaces + lines[i] <- paste0(" ", line) + } + } +} + +writeLines(lines, file_path) +cat("Fixed indentation for lines 408-1022\n") diff --git a/r_app/MANUAL_PIPELINE_RUNNER.R b/r_app/MANUAL_PIPELINE_RUNNER.R index 2cceb43..4850e75 100644 --- a/r_app/MANUAL_PIPELINE_RUNNER.R +++ b/r_app/MANUAL_PIPELINE_RUNNER.R @@ -73,6 +73,7 @@ # python 00_download_8band_pu_optimized.py [PROJECT] --date [DATE] --resolution 3 --cleanup # # Example: +# cd python_app # python 00_download_8band_pu_optimized.py angata --date 2026-02-04 --resolution 3 --cleanup # # COMMAND #2 - Batch Download (Multiple Dates): @@ -125,8 +126,6 @@ # # & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/10_create_per_field_tiffs.R angata # -# Example: -# & "C:\Program Files\R\R-4.4.3\bin\x64\Rscript.exe" r_app/10_create_per_field_tiffs.R angata # # COMMAND #2 - Specific Date Range: # diff --git a/r_app/extract_rds_only.R b/r_app/extract_rds_only.R deleted file mode 100644 index 61a6898..0000000 --- a/r_app/extract_rds_only.R +++ /dev/null @@ -1,104 +0,0 @@ -# EXTRACT_RDS_ONLY.R -# =================== -# Extract and combine daily CI values into combined_CI_data.rds -# Skips raster processing - assumes daily extracted files already exist -# -# Usage: Rscript r_app/extract_rds_only.R [project_dir] -# - project_dir: Project directory name (e.g., "angata", "aura", "chemba") -# -# Example: -# Rscript r_app/extract_rds_only.R angata - -suppressPackageStartupMessages({ - library(tidyverse) - library(here) -}) - -main <- function() { - # Capture command line arguments - args <- commandArgs(trailingOnly = TRUE) - - # Process project_dir argument - if (length(args) >= 1 && !is.na(args[1])) { - project_dir <- as.character(args[1]) - } else { - project_dir <- "angata" - } - - cat(sprintf("RDS Extraction: project=%s\n", project_dir)) - - # Source configuration - tryCatch({ - source("parameters_project.R") - }, error = function(e) { - warning("Default source files not found. Attempting to source from 'r_app' directory.") - tryCatch({ - source("r_app/parameters_project.R") - warning(paste("Successfully sourced files from 'r_app' directory.")) - }, error = function(e) { - stop("Failed to source parameters_project.R from both default and 'r_app' directories.") - }) - }) - - # Define paths for CI data - daily_CI_vals_dir <- file.path( - "laravel_app/storage/app", project_dir, - "Data/extracted_ci/daily_vals" - ) - - cumulative_CI_vals_dir <- file.path( - "laravel_app/storage/app", project_dir, - "Data/extracted_ci/cumulative_vals" - ) - - cat(sprintf("Daily CI values dir: %s\n", daily_CI_vals_dir)) - cat(sprintf("Cumulative CI values dir: %s\n\n", cumulative_CI_vals_dir)) - - # Check if daily CI directory exists and has files - if (!dir.exists(daily_CI_vals_dir)) { - stop(sprintf("ERROR: Daily CI directory not found: %s", daily_CI_vals_dir)) - } - - # List RDS files - files <- list.files(path = daily_CI_vals_dir, pattern = "^extracted_.*\\.rds$", full.names = TRUE) - - if (length(files) == 0) { - stop(sprintf("ERROR: No extracted CI values found in %s", daily_CI_vals_dir)) - } - - cat(sprintf("Found %d daily CI RDS files\n\n", length(files))) - - # Create cumulative directory if it doesn't exist - if (!dir.exists(cumulative_CI_vals_dir)) { - dir.create(cumulative_CI_vals_dir, recursive = TRUE) - cat(sprintf("Created directory: %s\n\n", cumulative_CI_vals_dir)) - } - - # Combine all RDS files - cat("Combining daily RDS files...\n") - combined_data <- files %>% - purrr::map(readRDS) %>% - purrr::list_rbind() %>% - dplyr::group_by(sub_field) - - # Save combined data - output_path <- file.path(cumulative_CI_vals_dir, "combined_CI_data.rds") - saveRDS(combined_data, output_path) - - cat(sprintf("✓ Combined %d daily files\n", length(files))) - cat(sprintf("✓ Total rows: %d\n", nrow(combined_data)) - cat(sprintf("✓ Saved to: %s\n\n", output_path)) - - # Summary - cat("Summary:\n") - cat(sprintf(" Fields: %d\n", n_distinct(combined_data$field, na.rm = TRUE))) - cat(sprintf(" Sub-fields: %d\n", n_distinct(combined_data$sub_field, na.rm = TRUE))) - cat(sprintf(" Total measurements: %d\n\n", nrow(combined_data))) - - cat("✓ RDS extraction complete!\n") - cat("Next: Run 02b_convert_rds_to_csv.R to convert to CSV\n") -} - -if (sys.nframe() == 0) { - main() -} diff --git a/r_app/old_scripts/kpi_utils.R b/r_app/old_scripts/kpi_utils.R index 552f6bc..be6ba36 100644 --- a/r_app/old_scripts/kpi_utils.R +++ b/r_app/old_scripts/kpi_utils.R @@ -9,11 +9,6 @@ # 5. Weed Presence Score # 6. Gap Filling Score -# Note: This file depends on functions from crop_messaging_utils.R: -# - safe_log() -# - calculate_cv() -# - calculate_spatial_autocorrelation() -# - calculate_change_percentages() # 1. Helper Functions # ----------------- @@ -676,7 +671,7 @@ calculate_growth_decline_kpi <- function(current_ci, previous_ci, field_boundari mean_change <- mean(ci_change) # Calculate spatial metrics - spatial_result <- calculate_spatial_autocorrelation(current_ci, field_vect) + spatial_result <- calculate_spatial_autocorrelation(current_field_ci, field_vect) cv_value <- calculate_cv(current_clean) # Determine risk level based on CI decline and spatial distribution diff --git a/r_app/parameters_project_OLD.R b/r_app/parameters_project_OLD.R deleted file mode 100644 index a5d9224..0000000 --- a/r_app/parameters_project_OLD.R +++ /dev/null @@ -1,1240 +0,0 @@ -# filepath: c:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane\r_app\parameters_project.R -# -# PARAMETERS_PROJECT.R -# ==================== -# This script defines project parameters, directory structures, and loads field boundaries. -# It establishes all the necessary paths and creates required directories for the SmartCane project. - -# 1. Load required libraries -# ------------------------- -suppressPackageStartupMessages({ - library(here) - library(readxl) - library(sf) - library(dplyr) - library(tidyr) - library(jsonlite) # For reading tiling_config.json -}) - -# 2. Client type mapping (for conditional script execution) -# --------------------------------------------------------- -# Maps project names to client types for pipeline control -# Client types: -# - "cane_supply": Runs Scripts 20,21,30,31,80,91 (full pipeline with Excel output) -# - "agronomic_support": Runs Scripts 80,90 only (KPI calculation + Word report) -# - "extension_service": (Future - not yet implemented) -# -# NOTE: This will eventually migrate to Laravel environment variables/database -# For now, maintain this mapping and update as projects are added -CLIENT_TYPE_MAP <- list( - "angata" = "cane_supply", - "aura" = "agronomic_support", - "chemba" = "cane_supply", - "xinavane" = "cane_supply", - "esa" = "cane_supply" -) - -get_client_type <- function(project_name) { - client_type <- CLIENT_TYPE_MAP[[project_name]] - if (is.null(client_type)) { - warning(sprintf("Project '%s' not in CLIENT_TYPE_MAP - defaulting to 'cane_supply'", project_name)) - return("cane_supply") - } - return(client_type) -} - -# 2b. Client-specific KPI configurations -# ---------------------------------------- -# Defines which KPIs and outputs are required for each client type -# This enables Script 80 to conditionally calculate only relevant metrics -# -# Structure: -# - kpi_calculations: Vector of KPI types to calculate for this client -# - outputs: Vector of output formats to generate (determines RDS/Excel naming) -# - requires_harvest_data: Boolean - whether Script 31 harvest predictions are needed -# - script_90_compatible: Boolean - whether output should match Script 90 expectations -# - script_91_compatible: Boolean - whether output should match Script 91 expectations -# -CLIENT_TYPE_CONFIGS <- list( - - # Aura (agronomic_support): Farm-level KPI summaries for weekly reports to agronomists - "agronomic_support" = list( - client_type = "agronomic_support", - description = "Farm-level KPI summaries for agronomic decision support", - kpi_calculations = c( - "field_uniformity", - "area_change", - "tch_forecasted", - "growth_decline", - "weed_presence", - "gap_filling" - ), - outputs = c( - "kpi_summary_tables", # Summary statistics for Script 90 report front page - "field_details" # Detailed field table for Script 90 report end section - ), - requires_harvest_data = FALSE, # Script 31 predictions not used - script_90_compatible = TRUE, # Output format matches Script 90 expectations - script_91_compatible = FALSE - ), - - # Cane Supply (cane_supply): Per-field analysis with harvest timing prediction - "cane_supply" = list( - client_type = "cane_supply", - description = "Per-field analysis with harvest prediction and phase assignment", - kpi_calculations = c( - "per_field_analysis", # Use 80_weekly_stats_utils.R for field-level statistics - "phase_assignment", # Assign growth phases (Germination, Tillering, Grand Growth, Maturation) - "harvest_prediction", # Include Script 31 harvest age predictions if available - "status_triggers" # Calculate field status (Normal, Monitor, Alert, Urgent) - ), - outputs = c( - "field_analysis_excel", # Excel file with per-field metrics - "field_analysis_summary" # Summary RDS for Script 91 report - ), - requires_harvest_data = TRUE, # harvest.xlsx is required for phase assignment - script_90_compatible = FALSE, - script_91_compatible = TRUE - ) -) - -#' Get KPI configuration for a specific client type -#' @param client_type Character string of client type (e.g., "agronomic_support", "cane_supply") -#' @return List containing configuration for that client type -get_client_kpi_config <- function(client_type) { - config <- CLIENT_TYPE_CONFIGS[[client_type]] - - if (is.null(config)) { - warning(sprintf("Client type '%s' not in CLIENT_TYPE_CONFIGS - defaulting to 'cane_supply'", client_type)) - return(CLIENT_TYPE_CONFIGS[["cane_supply"]]) - } - - return(config) -} - -# 3. Smart detection for tile-based vs single-file mosaic approach -# ---------------------------------------------------------------- -detect_tile_structure_from_merged_final <- function(merged_final_tif_dir, daily_tiles_split_dir = NULL) { - # PRIORITY 1: Check for tiling_config.json metadata file from script 10 - # This is the most reliable source since script 10 explicitly records its decision - - if (!is.null(daily_tiles_split_dir) && dir.exists(daily_tiles_split_dir)) { - # Try to find tiling_config.json in any grid-size subfolder - config_files <- list.files(daily_tiles_split_dir, - pattern = "tiling_config\\.json$", - recursive = TRUE, - full.names = TRUE) - - if (length(config_files) > 0) { - # Found a config file - use the most recent one - config_file <- config_files[which.max(file.info(config_files)$mtime)] - - tryCatch({ - config_json <- jsonlite::read_json(config_file) - return(list( - has_tiles = config_json$has_tiles %||% TRUE, - detected_tiles = character(), - total_files = 0, - source = "tiling_config.json", - grid_size = config_json$grid_size %||% "unknown" - )) - }, error = function(e) { - warning("Error reading tiling_config.json: ", e$message) - # Fall through to file-based detection - }) - } - } - - # PRIORITY 2: File-based detection (fallback if metadata not found) - # Check if merged_final_tif/ contains tile-named files OR grid-size subdirectories - - if (!dir.exists(merged_final_tif_dir)) { - return(list( - has_tiles = FALSE, - detected_tiles = character(), - total_files = 0, - source = "directory_not_found" - )) - } - - # First check if there are grid-size subdirectories (5x5, 10x10, etc.) - # This indicates the tiles are organized: merged_final_tif/{grid_size}/{DATE}/{DATE}_XX.tif - grid_subfolders <- list.dirs(merged_final_tif_dir, full.names = FALSE, recursive = FALSE) - grid_patterns <- grep("^\\d+x\\d+$", grid_subfolders, value = TRUE) - - if (length(grid_patterns) > 0) { - # Found grid-size subdirectories - tiles exist! - grid_size <- grid_patterns[1] - grid_dir <- file.path(merged_final_tif_dir, grid_size) - - # List sample tile files from the grid directory - sample_tiles <- list.files(grid_dir, pattern = "\\.tif$", recursive = TRUE)[1:3] - - return(list( - has_tiles = TRUE, - detected_tiles = sample_tiles, - total_files = length(sample_tiles), - source = "grid_subdirectory_detection", - grid_size = grid_size, - grid_path = grid_dir - )) - } - - # Fall back to checking for tile-named files directly in merged_final_tif - # List all .tif files in merged_final_tif - tif_files <- list.files(merged_final_tif_dir, pattern = "\\.tif$", full.names = FALSE) - - if (length(tif_files) == 0) { - return(list( - has_tiles = FALSE, - detected_tiles = character(), - total_files = 0, - source = "no_files_found" - )) - } - - # Check if ANY file matches tile naming pattern: *_XX.tif (where XX is 2 digits) - # Tile pattern examples: 2025-11-27_00.tif, 2025-11-27_01.tif, week_50_2024_00.tif - tile_pattern <- "_(\\d{2})\\.tif$" - tile_files <- tif_files[grepl(tile_pattern, tif_files)] - - has_tiles <- length(tile_files) > 0 - - return(list( - has_tiles = has_tiles, - detected_tiles = tile_files, - total_files = length(tif_files), - source = "file_pattern_detection" - )) -} - -# 4. Define project directory structure -# ----------------------------------- -# ============================================================================== -# CENTRALIZED PATH MANAGEMENT - setup_project_directories() -# ============================================================================== -# This function is the single source of truth for ALL file paths used across the pipeline. -# All scripts should call this function once at startup and use returned paths. -# This eliminates ~88 hardcoded file.path() calls scattered across 8 scripts. -# -# USAGE: -# paths <- setup_project_directories(project_dir) -# merged_tif_dir <- paths$merged_tif_folder -# daily_ci_dir <- paths$daily_ci_vals_dir -# kpi_output_dir <- paths$kpi_reports_dir -# -# TIERS (8-layer directory structure): -# Tier 1: Raw data (merged_tif) -# Tier 2: Per-field TIFFs (field_tiles, field_tiles_CI) -# Tier 3: CI Extraction (daily_ci_vals, cumulative_ci_vals) -# Tier 4: Growth Model (growth_model_interpolated) -# Tier 5: Mosaics (weekly_mosaic, weekly_tile_max) -# Tier 6: KPI & Reporting (kpi_reports_dir, kpi_field_stats_dir) -# Tier 7: Support (data, vrt, harvest, logs) -# Tier 8: Config & Metadata (field_boundaries_path, tiling_config_path) -# -# BENEFITS: -# ✓ Single source of truth (eliminates ~88 hardcoded file.path() calls) -# ✓ Auto-creates all directories (no scattered dir.create() calls) -# ✓ Easy to update storage structure globally -# ✓ Consistent naming across all 8 scripts -# ============================================================================== -setup_project_directories <- function(project_dir, data_source = "merged_tif") { - # =========================================================================== - # BASE DIRECTORIES (Foundation for all paths) - # =========================================================================== - laravel_storage_dir <- here("laravel_app", "storage", "app", project_dir) - - # =========================================================================== - # TIER 1: RAW DATA & INPUT PATHS (Script 00 - Python download output) - # =========================================================================== - merged_tif_folder <- here(laravel_storage_dir, "merged_tif") # 4-band raw GeoTIFFs from Planet - - # =========================================================================== - # TIER 2: TILING PATHS (Script 10 - Per-field tiff creation) - # =========================================================================== - # Per-field TIFF structure: field_tiles/{FIELD_NAME}/{YYYY-MM-DD}.tif - field_tiles_dir <- here(laravel_storage_dir, "field_tiles") - - # Per-field CI TIFFs (pre-computed, used by Script 40): field_tiles_CI/{FIELD_NAME}/{YYYY-MM-DD}.tif - field_tiles_ci_dir <- here(laravel_storage_dir, "field_tiles_CI") - - # Legacy tiling (for backward compatibility): daily_tiles_split/{grid_size}/{YYYY-MM-DD}/{YYYY-MM-DD}_XX.tif - daily_tiles_split_dir <- here(laravel_storage_dir, "daily_tiles_split") - - # =========================================================================== - # TIER 3: CI EXTRACTION PATHS (Script 20 - Canopy Index calculation) - # =========================================================================== - extracted_ci_base_dir <- here(laravel_storage_dir, "Data", "extracted_ci") - - # Daily CI values (cumulative RDS): combined_CI_data.rds - daily_ci_vals_dir <- here(extracted_ci_base_dir, "daily_vals") - - # Cumulative CI across time: All_pivots_Cumulative_CI_quadrant_year_v2.rds - cumulative_ci_vals_dir <- here(extracted_ci_base_dir, "cumulative_vals") - - # Per-field CI data for Python harvest prediction (Script 21): ci_data_for_python.csv - ci_for_python_dir <- here(extracted_ci_base_dir, "ci_data_for_python") - - # =========================================================================== - # TIER 4: GROWTH MODEL PATHS (Script 30 - Interpolation & smoothing) - # =========================================================================== - growth_model_interpolated_dir <- here(laravel_storage_dir, "growth_model_interpolated") - - # =========================================================================== - # TIER 5: MOSAIC PATHS (Script 40 - Weekly mosaics) - # =========================================================================== - # Per-field weekly mosaics (per-field architecture): weekly_mosaic/{FIELD}/{week_XX_YYYY}.tif - weekly_mosaic_dir <- here(laravel_storage_dir, "weekly_mosaic") - - # Tile-based weekly max (legacy): weekly_tile_max/{grid_size}/week_XX_YYYY.tif - weekly_tile_max_dir <- here(laravel_storage_dir, "weekly_tile_max") - - # =========================================================================== - # TIER 6: KPI & REPORTING PATHS (Scripts 80, 90, 91) - # =========================================================================== - reports_dir <- here(laravel_storage_dir, "reports") - kpi_reports_dir <- here(reports_dir, "kpis") # Where Script 80 outputs KPI CSV/RDS files - kpi_field_stats_dir <- here(kpi_reports_dir, "field_stats") # Per-field KPI details - kpi_field_analysis_dir <- here(kpi_reports_dir, "field_analysis") # Field-level analysis for Script 91 - - # =========================================================================== - # TIER 7: SUPPORT PATHS (Data, VRT, Harvest) - # =========================================================================== - data_dir <- here(laravel_storage_dir, "Data") - vrt_dir <- here(data_dir, "vrt") # Virtual Raster files created during CI extraction - harvest_dir <- here(data_dir, "HarvestData") # Harvest schedule data - log_dir <- here(laravel_storage_dir, "logs") # Log files - - # =========================================================================== - # TIER 8: CONFIG & METADATA PATHS - # =========================================================================== - # Field boundaries GeoJSON (same across all scripts) - field_boundaries_path <- here(data_dir, "pivot.geojson") - - # Tiling configuration metadata from Script 10 - tiling_config_path <- here(daily_tiles_split_dir, "tiling_config.json") - - # =========================================================================== - # CREATE ALL DIRECTORIES (once per pipeline run) - # =========================================================================== - all_dirs <- c( - # Tier 1 - merged_tif_folder, - # Tier 2 - field_tiles_dir, field_tiles_ci_dir, daily_tiles_split_dir, - # Tier 3 - extracted_ci_base_dir, daily_ci_vals_dir, cumulative_ci_vals_dir, ci_for_python_dir, - # Tier 4 - growth_model_interpolated_dir, - # Tier 5 - weekly_mosaic_dir, weekly_tile_max_dir, - # Tier 6 - reports_dir, kpi_reports_dir, kpi_field_stats_dir, kpi_field_analysis_dir, - # Tier 7 - data_dir, vrt_dir, harvest_dir, log_dir - ) - - for (dir_path in all_dirs) { - dir.create(dir_path, showWarnings = FALSE, recursive = TRUE) - } - - # =========================================================================== - # RETURN COMPREHENSIVE PATH LIST - # Scripts should source parameters_project.R and receive paths object like: - # paths <- setup_project_directories(project_dir) - # Then use: paths$merged_tif_folder, paths$daily_ci_vals_dir, etc. - # =========================================================================== - return(list( - # PROJECT ROOT - laravel_storage_dir = laravel_storage_dir, - - # TIER 1: Raw data - merged_tif_folder = merged_tif_folder, - - # TIER 2: Per-field TIFFs - field_tiles_dir = field_tiles_dir, - field_tiles_ci_dir = field_tiles_ci_dir, - daily_tiles_split_dir = daily_tiles_split_dir, - - # TIER 3: CI Extraction - extracted_ci_base_dir = extracted_ci_base_dir, - daily_ci_vals_dir = daily_ci_vals_dir, - cumulative_ci_vals_dir = cumulative_ci_vals_dir, - ci_for_python_dir = ci_for_python_dir, - - # TIER 4: Growth Model - growth_model_interpolated_dir = growth_model_interpolated_dir, - - # TIER 5: Mosaics - weekly_mosaic_dir = weekly_mosaic_dir, - weekly_tile_max_dir = weekly_tile_max_dir, - - # TIER 6: KPI & Reporting - reports_dir = reports_dir, - kpi_reports_dir = kpi_reports_dir, - kpi_field_stats_dir = kpi_field_stats_dir, - kpi_field_analysis_dir = kpi_field_analysis_dir, - - # TIER 7: Support - data_dir = data_dir, - vrt_dir = vrt_dir, - harvest_dir = harvest_dir, - log_dir = log_dir, - - # TIER 8: Config & Metadata - field_boundaries_path = field_boundaries_path, - tiling_config_path = tiling_config_path - )) -} - -# ============================================================================== -# TIER-BY-TIER PATH REFERENCE (for setup_project_directories output) -# ============================================================================== -# -# TIER 1: RAW DATA (Script 00 - Python download) -# paths$merged_tif_folder -# └─ {YYYY-MM-DD}.tif (4-band uint16 GeoTIFFs from Planet API) -# -# TIER 2: PER-FIELD TIFFS (Script 10) -# paths$field_tiles_dir/{FIELD_NAME}/{YYYY-MM-DD}.tif -# paths$field_tiles_ci_dir/{FIELD_NAME}/{YYYY-MM-DD}.tif -# paths$daily_tiles_split_dir/{grid_size}/{YYYY-MM-DD}/{YYYY-MM-DD}_XX.tif (legacy) -# -# TIER 3: CI EXTRACTION (Script 20) -# paths$daily_ci_vals_dir/combined_CI_data.rds -# paths$cumulative_ci_vals_dir/All_pivots_Cumulative_CI_quadrant_year_v2.rds -# paths$ci_for_python_dir/ci_data_for_python.csv (Script 21 output) -# -# TIER 4: GROWTH MODEL (Script 30) -# paths$growth_model_interpolated_dir/ (RDS files with interpolated CI) -# -# TIER 5: MOSAICS (Script 40) -# paths$weekly_mosaic_dir/{FIELD_NAME}/week_XX_YYYY.tif -# paths$weekly_tile_max_dir/{grid_size}/week_XX_YYYY_00.tif (legacy) -# -# TIER 6: KPI & REPORTING (Scripts 80, 90, 91) -# paths$kpi_reports_dir/ (KPI outputs from Script 80) -# paths$kpi_field_stats_dir/ (Per-field KPI RDS) -# paths$kpi_field_analysis_dir/ (Analysis RDS for Script 91) -# paths$reports_dir/ (Word/HTML reports) -# -# TIER 7: SUPPORT (Various scripts) -# paths$data_dir/pivot.geojson (Field boundaries) -# paths$data_dir/harvest.xlsx (Harvest schedule) -# paths$vrt_dir/ (Virtual raster files) -# paths$harvest_dir/ (Harvest predictions from Python) -# paths$log_dir/ (Pipeline logs) -# -# TIER 8: CONFIG & METADATA -# paths$field_boundaries_path (Full path to pivot.geojson) -# paths$tiling_config_path (Metadata from Script 10) -# -# ============================================================================== - -# 5. Utility Functions -# ---------------------- -# NOTE: load_field_boundaries() and load_harvesting_data() are defined in 00_common_utils.R -# to avoid duplication. They are sourced after parameters_project.R and used by all scripts. - -# 6. Load harvesting data -# --------------------- -load_harvesting_data <- function(data_dir) { - harvest_file <- here(data_dir, "harvest.xlsx") - - if (!file.exists(harvest_file)) { - warning(paste("Harvest data file not found at path:", harvest_file)) - return(NULL) - } - - # Helper function to parse dates with multiple format detection - parse_flexible_date <- function(x) { - if (is.na(x) || is.null(x)) return(NA_real_) - if (inherits(x, "Date")) return(x) - if (inherits(x, "POSIXct")) return(as.Date(x)) - - # If it's numeric (Excel date serial), convert directly - if (is.numeric(x)) { - return(as.Date(x, origin = "1899-12-30")) - } - - # Try character conversion with multiple formats - x_char <- as.character(x) - - # Try common formats: YYYY-MM-DD, DD/MM/YYYY, MM/DD/YYYY, YYYY-MM-DD HH:MM:SS - formats <- c("%Y-%m-%d", "%d/%m/%Y", "%m/%d/%Y", "%Y-%m-%d %H:%M:%S") - - for (fmt in formats) { - result <- suppressWarnings(as.Date(x_char, format = fmt)) - if (!is.na(result)) return(result) - } - - # If all else fails, return NA - return(NA) - } - - tryCatch({ - harvesting_data <- read_excel(harvest_file) %>% - dplyr::select( - c( - "field", - "sub_field", - "year", - "season_start", - "season_end", - "age", - "sub_area", - "tonnage_ha" - ) - ) %>% - mutate( - field = as.character(field), - sub_field = as.character(sub_field), - year = as.numeric(year), - season_start = sapply(season_start, parse_flexible_date), - season_end = sapply(season_end, parse_flexible_date), - season_start = as.Date(season_start, origin = "1970-01-01"), - season_end = as.Date(season_end, origin = "1970-01-01"), - age = as.numeric(age), - sub_area = as.character(sub_area), - tonnage_ha = as.numeric(tonnage_ha) - ) %>% - mutate( - season_end = case_when( - season_end > Sys.Date() ~ Sys.Date(), - is.na(season_end) ~ Sys.Date(), - TRUE ~ season_end - ), - age = round(as.numeric(season_end - season_start) / 7, 0) - ) - - return(harvesting_data) - }, error = function(e) { - warning(paste("Error loading harvesting data:", e$message)) - return(NULL) - }) -} - -# 5. Define logging functions globally first -# --------------------------------------- -# Create a simple default log function in case setup_logging hasn't been called yet -log_message <- function(message, level = "INFO") { - timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") - formatted_message <- paste0("[", level, "] ", timestamp, " - ", message) - cat(formatted_message, "\n") -} - -log_head <- function(list, level = "INFO") { - log_message(paste(capture.output(str(head(list))), collapse = "\n"), level) -} - -# 8. Set up full logging system with file output -# ------------------------------------------- -setup_logging <- function(log_dir) { - log_file <- here(log_dir, paste0(format(Sys.Date(), "%Y%m%d"), ".log")) - - # Create enhanced log functions - log_message <- function(message, level = "INFO") { - timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S") - formatted_message <- paste0("[", level, "] ", timestamp, " - ", message) - cat(formatted_message, "\n", file = log_file, append = TRUE) - - # Also print to console for debugging - if (level %in% c("ERROR", "WARNING")) { - cat(formatted_message, "\n") - } - } - - log_head <- function(list, level = "INFO") { - log_message(paste(capture.output(str(head(list))), collapse = "\n"), level) - } - - # Update the global functions with the enhanced versions - assign("log_message", log_message, envir = .GlobalEnv) - assign("log_head", log_head, envir = .GlobalEnv) - - return(list( - log_file = log_file, - log_message = log_message, - log_head = log_head - )) -} - -# 8. HELPER FUNCTIONS FOR COMMON CALCULATIONS -# ----------------------------------------------- -# Centralized functions to reduce duplication across scripts - -# Get ISO week and year from a date -get_iso_week <- function(date) { - as.numeric(format(date, "%V")) -} - -get_iso_year <- function(date) { - as.numeric(format(date, "%G")) -} - -# Get both ISO week and year as a list -get_iso_week_year <- function(date) { - list( - week = as.numeric(format(date, "%V")), - year = as.numeric(format(date, "%G")) - ) -} - -# Format week/year into a readable label -format_week_label <- function(date, separator = "_") { - wwy <- get_iso_week_year(date) - sprintf("week%02d%s%d", wwy$week, separator, wwy$year) -} - -# Auto-detect mosaic mode -# For per-field architecture, always returns "single-file" (weekly_mosaic/{FIELD}/week_*.tif) -detect_mosaic_mode <- function(project_dir) { - # Per-field architecture uses single-file mosaics organized per-field - weekly_mosaic <- file.path("laravel_app", "storage", "app", project_dir, "weekly_mosaic") - if (dir.exists(weekly_mosaic)) { - return("single-file") # Per-field structure - } - return("unknown") -} - -# Auto-detect grid size from tile directory structure -# For per-field architecture, returns "unknown" since grid-based organization is legacy -detect_grid_size <- function(project_dir) { - # Per-field architecture doesn't use grid-based organization anymore - return("unknown") -} - -# Build storage paths consistently across all scripts -get_project_storage_path <- function(project_dir, subdir = NULL) { - base <- file.path("laravel_app", "storage", "app", project_dir) - if (!is.null(subdir)) file.path(base, subdir) else base -} - -get_mosaic_dir <- function(project_dir, mosaic_mode = "auto") { - # Per-field architecture always uses weekly_mosaic (single-file, per-field organization) - get_project_storage_path(project_dir, "weekly_mosaic") -} - -get_kpi_dir <- function(project_dir, client_type) { - subdir <- if (client_type == "agronomic_support") "field_level" else "field_analysis" - get_project_storage_path(project_dir, file.path("reports", "kpis", subdir)) -} - -# Logging functions moved to 00_common_utils.R -# - smartcane_log() — Main logging function with level prefix -# - smartcane_debug() — Conditional debug logging -# - smartcane_warn() — Warning wrapper -# Import with: source("r_app/00_common_utils.R") - -# ============================================================================ -# PHASE 3 & 4: OPTIMIZATION & DOCUMENTATION -# ============================================================================ - -# System Constants -# ---------------- -# Define once, use everywhere - -RSCRIPT_PATH <- "C:\\Program Files\\R\\R-4.4.3\\bin\\x64\\Rscript.exe" -# Used in run_full_pipeline.R for calling R scripts via system() - -# Data Source Documentation -# --------------------------- -# Explains the two satellite data formats and when to use each -# -# SmartCane uses PlanetScope imagery from Planet Labs API in two formats: -# -# 1. merged_tif (4-band): -# - Standard format: Red, Green, Blue, Near-Infrared -# - Size: ~150-200 MB per date -# - Use case: Agronomic support, general crop health monitoring -# - Projects: aura, xinavane -# - Cloud handling: Basic cloud masking from Planet metadata -# -# 2. merged_tif_8b (8-band with cloud confidence): -# - Enhanced format: 4-band imagery + 4-band UDM2 cloud mask -# - UDM2 bands: Clear, Snow, Shadow, Light Haze -# - Size: ~250-350 MB per date -# - Use case: Harvest prediction, supply chain optimization -# - Projects: angata, chemba, esa (cane_supply clients) -# - Cloud handling: Per-pixel cloud confidence from Planet UDM2 -# - Why: Cane supply chains need precise confidence to predict harvest dates -# (don't want to predict based on cloudy data) -# -# The system auto-detects which is available via detect_data_source() - -# Mosaic Mode Documentation -# -------------------------- -# SmartCane supports two ways to store and process weekly mosaics: -# -# 1. Single-file mosaic ("single-file"): -# - One GeoTIFF per week: weekly_mosaic/week_02_2026.tif -# - 5 bands per file: R, G, B, NIR, CI (Canopy Index) -# - Size: ~300-500 MB per week -# - Pros: Simpler file management, easier full-field visualization -# - Cons: Slower for field-specific queries, requires loading full raster -# - Best for: Agronomic support (aura) with <100 fields -# - Script 04 output: 5-band single-file mosaic -# -# 2. Tiled mosaic ("tiled"): -# - Grid of tiles per week: weekly_tile_max/5x5/week_02_2026_{TT}.tif -# - Example: 25 files (5×5 grid) × 5 bands = 125 individual tiffs -# - Size: ~15-20 MB per tile, organized in folders -# - Pros: Parallel processing, fast field lookups, scales to 1000+ fields -# - Cons: More file I/O, requires tile-to-field mapping metadata -# - Best for: Cane supply (angata, chemba) with 500+ fields -# - Script 04 output: Per-tile tiff files in weekly_tile_max/{grid}/ -# - Tile assignment: Field boundaries mapped to grid coordinates -# -# The system auto-detects which is available via detect_mosaic_mode() - -# Client Type Documentation -# -------------------------- -# SmartCane runs different analysis pipelines based on client_type: -# -# CLIENT_TYPE: cane_supply -# Purpose: Optimize sugar mill supply chain (harvest scheduling) -# Scripts run: 20 (CI), 21 (RDS to CSV), 30 (Growth), 31 (Harvest pred), 40 (Mosaic), 80 (KPI), 91 (Excel) -# Outputs: -# - Per-field analysis: field status, growth phase, harvest readiness -# - Excel reports (Script 91): Detailed metrics for logistics planning -# - KPI directory: reports/kpis/field_analysis/ (one RDS per week) -# Harvest data: Required (harvest.xlsx - planting dates for phase assignment) -# Data source: merged_tif_8b (uses cloud confidence for confidence) -# Mosaic mode: tiled (scales to 500+ fields) -# Projects: angata, chemba, xinavane, esa -# -# CLIENT_TYPE: agronomic_support -# Purpose: Provide weekly crop health insights to agronomists -# Scripts run: 80 (KPI), 90 (Word report) -# Outputs: -# - Farm-level KPI summaries (no per-field breakdown) -# - Word reports (Script 90): Charts and trends for agronomist decision support -# - KPI directory: reports/kpis/field_level/ (one RDS per week) -# Harvest data: Not used -# Data source: merged_tif (simpler, smaller) -# Mosaic mode: single-file (100-200 fields) -# Projects: aura -# - -# Detect data source (merged_tif vs merged_tif_8b) based on availability -# Returns the first available source; defaults to merged_tif_8b if neither exists -detect_data_source <- function(project_dir) { - # Data source is always merged_tif for consistency - return("merged_tif") -} - -# Check KPI completeness for a reporting period -# Returns: List with kpis_df (data.frame), missing_count, and all_complete (boolean) -# This replaces duplicate KPI checking logic in run_full_pipeline.R (lines ~228-270, ~786-810) -check_kpi_completeness <- function(project_dir, client_type, end_date, reporting_weeks_needed) { - kpi_dir <- get_kpi_dir(project_dir, client_type) - - kpis_needed <- data.frame() - - for (weeks_back in 0:(reporting_weeks_needed - 1)) { - check_date <- end_date - (weeks_back * 7) - wwy <- get_iso_week_year(check_date) - - # Build week pattern and check if it exists - week_pattern <- sprintf("week%02d_%d", wwy$week, wwy$year) - files_this_week <- list.files(kpi_dir, pattern = week_pattern) - has_kpis <- length(files_this_week) > 0 - - # Track missing weeks - kpis_needed <- rbind(kpis_needed, data.frame( - week = wwy$week, - year = wwy$year, - date = check_date, - has_kpis = has_kpis, - pattern = week_pattern, - file_count = length(files_this_week) - )) - - # Debug logging - smartcane_debug(sprintf( - "Week %02d/%d (%s): %s (%d files)", - wwy$week, wwy$year, format(check_date, "%Y-%m-%d"), - if (has_kpis) "✓ FOUND" else "✗ MISSING", - length(files_this_week) - )) - } - - # Summary statistics - missing_count <- sum(!kpis_needed$has_kpis) - all_complete <- missing_count == 0 - - return(list( - kpis_df = kpis_needed, - kpi_dir = kpi_dir, - missing_count = missing_count, - missing_weeks = kpis_needed[!kpis_needed$has_kpis, ], - all_complete = all_complete - )) -} - -# ============================================================================== -# HELPER FUNCTIONS FOR run_full_pipeline.R PATH VERIFICATION (SC-116) -# ============================================================================== -# These functions replace hardcoded file.path() calls in run_full_pipeline.R -# with centralized, testable helper functions. Each function verifies a specific -# output directory for a pipeline stage. - -#' Get verification path for Script 31 harvest output -#' -#' @param project_dir Character. Project name (e.g., "angata", "aura") -#' @param week_num Integer. ISO week number (01-53) -#' @param year_num Integer. Year (e.g., 2026) -#' @return Character. Full path to expected harvest imminent CSV file -#' @details -#' Script 31 generates: {project}_{project}_harvest_imminent_week_{WW}_{YYYY}.csv -#' Location: laravel_app/storage/app/{project}/reports/kpis/field_stats/ -#' -get_harvest_output_path <- function(project_dir, week_num, year_num) { - file.path( - "laravel_app", "storage", "app", project_dir, "reports", "kpis", "field_stats", - sprintf("%s_harvest_imminent_week_%02d_%d.csv", project_dir, week_num, year_num) - ) -} - -#' Check if harvest output file exists for a given week -#' -#' @param project_dir Character. Project name -#' @param week_num Integer. ISO week number -#' @param year_num Integer. Year -#' @return Logical. TRUE if file exists, FALSE otherwise -#' -check_harvest_output_exists <- function(project_dir, week_num, year_num) { - path <- get_harvest_output_path(project_dir, week_num, year_num) - file.exists(path) -} - -#' Get expected output directory for a mosaic verification based on mode -#' -#' @param project_dir Character. Project name -#' @param mosaic_mode Character. Either "tiled" or "single-file" -#' @return Character. Full path to mosaic directory -#' -#' @details -#' Tiled: laravel_app/storage/app/{project}/weekly_tile_max/ -#' Single-file: laravel_app/storage/app/{project}/weekly_mosaic/ -#' -get_mosaic_verification_dir <- function(project_dir, mosaic_mode) { - base <- file.path("laravel_app", "storage", "app", project_dir) - - if (mosaic_mode == "tiled") { - file.path(base, "weekly_tile_max") - } else { - # Default to single-file (single-file is standard for per-field architecture) - file.path(base, "weekly_mosaic") - } -} - -#' Check if mosaic files exist for a specific week -#' -#' @param project_dir Character. Project name -#' @param week_num Integer. ISO week number -#' @param year_num Integer. Year -#' @param mosaic_mode Character. "tiled" or "single-file" -#' @return List with created (logical), file_count (int), and sample_files (char vector) -#' -check_mosaic_exists <- function(project_dir, week_num, year_num, mosaic_mode) { - mosaic_dir <- get_mosaic_verification_dir(project_dir, mosaic_mode) - - if (!dir.exists(mosaic_dir)) { - return(list(created = FALSE, file_count = 0, sample_files = character())) - } - - week_pattern <- sprintf("week_%02d_%d\\.tif", week_num, year_num) - # Search recursively for per-field architecture support - mosaic_files <- list.files(mosaic_dir, pattern = week_pattern, recursive = TRUE, full.names = FALSE) - - list( - created = length(mosaic_files) > 0, - file_count = length(mosaic_files), - sample_files = head(mosaic_files, 3) # First 3 files as sample - ) -} - -# 9. Initialize the project -# ---------------------- -# Export project directories and settings -initialize_project <- function(project_dir, data_source = "merged_tif") { - # Set up directory structure, passing data_source to select TIF folder - dirs <- setup_project_directories(project_dir, data_source = data_source) - - # Set up logging - logging <- setup_logging(dirs$log_dir) - - # Load field boundaries - boundaries <- load_field_boundaries(dirs$data_dir) - - # Load harvesting data - harvesting_data <- load_harvesting_data(dirs$data_dir) - - # Return all initialized components - return(c( - dirs, - list( - logging = logging, - field_boundaries = boundaries$field_boundaries, - field_boundaries_sf = boundaries$field_boundaries_sf, - harvesting_data = harvesting_data - ) - )) -} - -# When script is sourced, initialize with the global project_dir variable if it exists -if (exists("project_dir")) { - # Now we can safely log before initialization - log_message(paste("Initializing project with directory:", project_dir)) - - # Use data_source if it exists (passed from 02_ci_extraction.R), otherwise use default - data_src <- if (exists("data_source")) data_source else "merged_tif" - log_message(paste("Using data source directory:", data_src)) - - project_config <- initialize_project(project_dir, data_source = data_src) - - # Expose all variables to the global environment - list2env(project_config, envir = .GlobalEnv) - - # Log project initialization completion with tile mode info - log_message(paste("Project initialized with directory:", project_dir)) - if (exists("use_tile_mosaic")) { - mosaic_mode <- if (use_tile_mosaic) "TILE-BASED" else "SINGLE-FILE" - log_message(paste("Mosaic mode detected:", mosaic_mode)) - if (exists("tile_detection_info") && !is.null(tile_detection_info)) { - log_message(paste(" - Detection source:", tile_detection_info$detected_source)) - log_message(paste(" - Grid size:", tile_detection_info$grid_size)) - log_message(paste(" - Detected files in storage:", tile_detection_info$detected_count)) - if (length(tile_detection_info$sample_tiles) > 0) { - log_message(paste(" - Sample tile files:", paste(tile_detection_info$sample_tiles, collapse = ", "))) - } - } - } -} else { - warning("project_dir variable not found. Please set project_dir before sourcing parameters_project.R") -} - - - -#' Safe Logging Function -#' -#' Generic logging with [LEVEL] prefix. Works standalone without any framework. -#' Consistent with SmartCane logging standard. -#' -#' @param message The message to log -#' @param level The log level (default: "INFO"). Options: "INFO", "WARNING", "ERROR", "DEBUG" -#' @return NULL (invisible, used for side effects) -#' -#' @examples -#' safe_log("Processing started", "INFO") -#' safe_log("Check input file", "WARNING") -#' safe_log("Failed to load data", "ERROR") -#' -safe_log <- function(message, level = "INFO") { - prefix <- sprintf("[%s]", level) - cat(sprintf("%s %s\n", prefix, message)) -} - -#' SmartCane Debug Logging (Conditional) -#' -#' Logs DEBUG-level messages only if verbose=TRUE or SMARTCANE_DEBUG env var is set. -#' Useful for development/troubleshooting without cluttering normal output. -#' -#' @param message The message to log -#' @param verbose Whether to output regardless of SMARTCANE_DEBUG (default: FALSE) -#' @return NULL (invisible, used for side effects) -#' -#' @examples -#' smartcane_debug("Processing field 1", verbose = FALSE) # Only if SMARTCANE_DEBUG=TRUE -#' smartcane_debug("Detailed state info", verbose = TRUE) # Always outputs -#' -smartcane_debug <- function(message, verbose = FALSE) { - if (!verbose && Sys.getenv("SMARTCANE_DEBUG") != "TRUE") { - return(invisible(NULL)) - } - safe_log(message, level = "DEBUG") -} - -#' SmartCane Warning Logging -#' -#' Logs WARN-level messages. Convenience wrapper around safe_log(). -#' -#' @param message The message to log -#' @return NULL (invisible, used for side effects) -#' -#' @examples -#' smartcane_warn("Check data format before proceeding") -#' -smartcane_warn <- function(message) { - safe_log(message, level = "WARN") -} - -#' Load Field Boundaries from GeoJSON -#' -#' Loads field polygon geometries from GeoJSON file (pivot.geojson or pivot_2.geojson). -#' Handles CRS validation and column standardization. -#' -#' @param data_dir Directory containing GeoJSON file -#' @return List with elements: -#' - field_boundaries_sf: sf (Simple Features) object -#' - field_boundaries: terra SpatVect object (if conversion successful, else sf fallback) -#' -#' @details -#' Automatically selects pivot_2.geojson for ESA project during CI extraction, -#' otherwise uses pivot.geojson. Handles both multi-polygon and simple polygon geometries. -#' -#' @examples -#' boundaries <- load_field_boundaries("laravel_app/storage/app/angata") -#' head(boundaries$field_boundaries_sf) -#' -load_field_boundaries <- function(data_dir) { - # Choose field boundaries file based on project and script type - # ESA project uses pivot_2.geojson ONLY for scripts 02-03 (CI extraction & growth model) - # All other scripts (including 04-mosaic, 09-KPIs, 10-reports) use pivot.geojson - use_pivot_2 <- exists("project_dir") && project_dir == "esa" && - exists("ci_extraction_script") # ci_extraction_script flag set by scripts 02-03 - - if (use_pivot_2) { - field_boundaries_path <- file.path(data_dir, "pivot_2.geojson") - } else { - field_boundaries_path <- file.path(data_dir, "pivot.geojson") - } - - if (!file.exists(field_boundaries_path)) { - stop(paste("Field boundaries file not found at path:", field_boundaries_path)) - } - - tryCatch({ - # Read GeoJSON with explicit CRS handling - field_boundaries_sf <- st_read(field_boundaries_path, quiet = TRUE) - - # Remove OBJECTID column immediately if it exists - if ("OBJECTID" %in% names(field_boundaries_sf)) { - field_boundaries_sf <- field_boundaries_sf %>% select(-OBJECTID) - } - - # **CRITICAL**: Repair invalid geometries (degenerate vertices, self-intersections, etc.) - # This must happen BEFORE any spatial operations (CRS transform, intersect, crop, etc.) - # to prevent S2 geometry validation errors during downstream processing - field_boundaries_sf <- repair_geojson_geometries(field_boundaries_sf) - - # Validate and fix CRS if needed - tryCatch({ - # Simply assign WGS84 if not already set (safe approach) - if (is.na(sf::st_crs(field_boundaries_sf)$epsg)) { - st_crs(field_boundaries_sf) <- 4326 - warning("CRS was missing, assigned WGS84 (EPSG:4326)") - } - }, error = function(e) { - tryCatch({ - st_crs(field_boundaries_sf) <<- 4326 - }, error = function(e2) { - warning(paste("Could not set CRS:", e2$message)) - }) - }) - - # Handle column names - accommodate optional sub_area column - if ("sub_area" %in% names(field_boundaries_sf)) { - field_boundaries_sf <- field_boundaries_sf %>% - dplyr::select(field, sub_field, sub_area) %>% - sf::st_set_geometry("geometry") - } else { - field_boundaries_sf <- field_boundaries_sf %>% - dplyr::select(field, sub_field) %>% - sf::st_set_geometry("geometry") - } - - # Convert to terra vector if possible, otherwise use sf - field_boundaries <- tryCatch({ - field_boundaries_terra <- terra::vect(field_boundaries_sf) - crs_value <- tryCatch(terra::crs(field_boundaries_terra), error = function(e) NULL) - crs_str <- if (!is.null(crs_value)) as.character(crs_value) else "" - - if (is.null(crs_value) || length(crs_value) == 0 || nchar(crs_str) == 0) { - terra::crs(field_boundaries_terra) <- "EPSG:4326" - warning("Terra object CRS was empty, assigned WGS84 (EPSG:4326)") - } - field_boundaries_terra - - }, error = function(e) { - warning(paste("Terra conversion failed, using sf object instead:", e$message)) - field_boundaries_sf - }) - - return(list( - field_boundaries_sf = field_boundaries_sf, - field_boundaries = field_boundaries - )) - }, error = function(e) { - cat("[DEBUG] Error in load_field_boundaries:\n") - cat(" Message:", e$message, "\n") - cat(" Call:", deparse(e$call), "\n") - stop(paste("Error loading field boundaries:", e$message)) - }) -} - - - -#' Generate a Sequence of Dates for Processing -#' -#' Creates a date range from start_date to end_date and extracts week/year info. -#' Used by Scripts 20, 30, 40 to determine data processing windows. -#' -#' @param end_date The end date for the sequence (Date object or "YYYY-MM-DD" string) -#' @param offset Number of days to look back from end_date (e.g., 7 for one week) -#' @return A list containing: -#' - week: ISO week number of start_date -#' - year: ISO year of start_date -#' - days_filter: Vector of dates in "YYYY-MM-DD" format -#' - start_date: Start date as Date object -#' - end_date: End date as Date object -#' -#' @details -#' IMPORTANT: Uses `lubridate::week()` and `lubridate::year()` which return -#' ISO week numbers (week 1 starts on Monday). For ISO week-based calculations, -#' use `lubridate::isoweek()` and `lubridate::isoyear()` instead. -#' -#' @examples -#' dates <- date_list(as.Date("2025-01-15"), offset = 7) -#' # Returns: week=2, year=2025, days_filter = c("2025-01-09", ..., "2025-01-15") -#' -#' dates <- date_list("2025-12-31", offset = 14) -#' # Handles string input and returns 14 days of data -#' -date_list <- function(end_date, offset) { - # Input validation - if (!lubridate::is.Date(end_date)) { - end_date <- as.Date(end_date) - if (is.na(end_date)) { - stop("Invalid end_date provided. Expected a Date object or a string convertible to Date.") - } - } - - offset <- as.numeric(offset) - if (is.na(offset) || offset < 1) { - stop("Invalid offset provided. Expected a positive number.") - } - - # Calculate date range - offset <- offset - 1 # Adjust offset to include end_date - start_date <- end_date - lubridate::days(offset) - - # Extract ISO week and year information (from END date for reporting period) - week <- lubridate::isoweek(end_date) - year <- lubridate::isoyear(end_date) - - # Generate sequence of dates - days_filter <- seq(from = start_date, to = end_date, by = "day") - days_filter <- format(days_filter, "%Y-%m-%d") # Format for consistent filtering - - # Log the date range - safe_log(paste("Date range generated from", start_date, "to", end_date)) - - return(list( - "week" = week, - "year" = year, - "days_filter" = days_filter, - "start_date" = start_date, - "end_date" = end_date - )) -} - -# ============================================================================== -#' Repair Invalid GeoJSON Geometries -#' -#' Fixes common geometry issues in GeoJSON/sf objects: -#' - Degenerate vertices (duplicate points) -#' - Self-intersecting polygons -#' - Invalid ring orientation -#' - Empty or NULL geometries -#' -#' Uses sf::st_make_valid() with buffer trick as fallback. -#' -#' @param sf_object sf object (GeoDataFrame) with potentially invalid geometries -#' @return sf object with repaired geometries -#' -#' @details -#' **Why this matters:** -#' Pivot GeoJSON files sometimes contain degenerate vertices or self-intersecting -#' rings from manual editing or GIS data sources. These cause errors when using -#' S2 geometry (strict validation) during cropping operations. -#' -#' **Repair strategy (priority order):** -#' 1. Try st_make_valid() - GEOS-based repair (most reliable) -#' 2. Fallback: st_union() + buffer(0) - Forces polygon validity -#' 3. Last resort: Silently keep original if repair fails -#' -#' @examples -#' \dontrun{ -#' fields <- st_read("pivot.geojson") -#' fields_fixed <- repair_geojson_geometries(fields) -#' cat(paste("Fixed geometries: before=", -#' nrow(fields[!st_is_valid(fields), ]), -#' ", after=", -#' nrow(fields_fixed[!st_is_valid(fields_fixed), ]))) -#' } -#' -repair_geojson_geometries <- function(sf_object) { - if (!inherits(sf_object, "sf")) { - stop("Input must be an sf (Simple Features) object") - } - - # Count invalid geometries BEFORE repair - invalid_before <- sum(!sf::st_is_valid(sf_object), na.rm = TRUE) - - if (invalid_before == 0) { - safe_log("All geometries already valid - no repair needed", "INFO") - return(sf_object) - } - - safe_log(paste("Found", invalid_before, "invalid geometries - attempting repair"), "WARNING") - - # STRATEGY: Apply st_make_valid() to entire sf object (most reliable for GEOS) - # This handles degenerate vertices, self-intersections, invalid rings while preserving all features - repaired <- tryCatch({ - # st_make_valid() on entire sf object preserves all features and attributes - repaired_geom <- sf::st_make_valid(sf_object) - - # Verify we still have the same number of rows - if (nrow(repaired_geom) != nrow(sf_object)) { - warning("st_make_valid() changed number of features - attempting row-wise repair") - - # Fallback: Repair row-by-row to maintain original structure - repaired_geom <- sf_object - for (i in seq_len(nrow(sf_object))) { - tryCatch({ - if (!sf::st_is_valid(sf_object[i, ])) { - repaired_geom[i, ] <- sf::st_make_valid(sf_object[i, ]) - } - }, error = function(e) { - safe_log(paste("Could not repair row", i, "-", e$message), "WARNING") - }) - } - } - - safe_log("✓ st_make_valid() successfully repaired geometries", "INFO") - repaired_geom - }, error = function(e) { - safe_log(paste("st_make_valid() failed:", e$message), "WARNING") - NULL - }) - - # If repair failed, keep original - if (is.null(repaired)) { - safe_log(paste("Could not repair", invalid_before, "invalid geometries - keeping original"), - "WARNING") - return(sf_object) - } - - # Count invalid geometries AFTER repair - invalid_after <- sum(!sf::st_is_valid(repaired), na.rm = TRUE) - safe_log(paste("Repair complete: before =", invalid_before, ", after =", invalid_after), "INFO") - - return(repaired) -} - -# ============================================================================== -# END 00_COMMON_UTILS.R -# ============================================================================== diff --git a/r_app/report_utils.R b/r_app/report_utils.R index cc1a9aa..0d6ffe3 100644 --- a/r_app/report_utils.R +++ b/r_app/report_utils.R @@ -244,7 +244,8 @@ ci_plot <- function(pivotName, # Create spans for borders joined_spans2 <- field_boundaries %>% - sf::st_transform(sf::st_crs(pivotShape)) %>% dplyr::filter(field %in% pivotName) + sf::st_transform(sf::st_crs(pivotShape)) %>% + dplyr::filter(field %in% pivotName) # Create the maps for different timepoints CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2,