diff --git a/r_app/80_utils_agronomic_support.R b/r_app/80_utils_agronomic_support.R index aecc837..b874ee9 100644 --- a/r_app/80_utils_agronomic_support.R +++ b/r_app/80_utils_agronomic_support.R @@ -80,7 +80,8 @@ prepare_predictions <- function(harvest_model, field_data, scenario = "optimisti #' @param ci_band Raster band with CI values #' #' @return Data frame with field_idx, cv_value, morans_i, uniformity_score, interpretation -calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_sf, ci_band = NULL) { +calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_sf, ci_band = NULL, + mosaic_dir = NULL, week_file = NULL) { result <- data.frame( field_idx = integer(), cv_value = numeric(), @@ -90,6 +91,9 @@ calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_ stringsAsFactors = FALSE ) + # Determine if we're using per-field structure + is_per_field <- !is.null(mosaic_dir) && !is.null(week_file) + for (field_idx in seq_len(nrow(field_boundaries_sf))) { ci_pixels <- ci_pixels_by_field[[field_idx]] @@ -107,10 +111,31 @@ calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_ cv_val <- calculate_cv(ci_pixels) + # Calculate Moran's I morans_i <- NA_real_ - if (!is.null(ci_band) && inherits(ci_band, "SpatRaster")) { + if (is_per_field) { + # Load individual field raster for per-field structure + field_name <- field_boundaries_sf$field[field_idx] + field_mosaic_path <- file.path(mosaic_dir, field_name, week_file) + + if (file.exists(field_mosaic_path)) { + tryCatch({ + field_raster <- terra::rast(field_mosaic_path)[["CI"]] + single_field <- field_boundaries_sf[field_idx, ] + morans_result <- calculate_spatial_autocorrelation(field_raster, single_field) + + if (is.list(morans_result)) { + morans_i <- morans_result$morans_i + } else { + morans_i <- morans_result + } + }, error = function(e) { + message(paste(" Warning: Spatial autocorrelation failed for field", field_name, ":", e$message)) + }) + } + } else if (!is.null(ci_band) && inherits(ci_band, "SpatRaster")) { + # Use single raster for single-file structure tryCatch({ - # Get single field geometry single_field <- field_boundaries_sf[field_idx, ] morans_result <- calculate_spatial_autocorrelation(ci_band, single_field) @@ -121,12 +146,11 @@ calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_ } }, error = function(e) { message(paste(" Warning: Spatial autocorrelation failed for field", field_idx, ":", e$message)) - morans_i <<- NA_real_ }) } # Normalize CV (0-1 scale, invert so lower CV = higher score) - cv_normalized <- min(cv_val / 0.3, 1) # 0.3 = threshold for CV + cv_normalized <- min(cv_val / 0.3, 1) cv_score <- 1 - cv_normalized # Normalize Moran's I (-1 to 1 scale, shift to 0-1) @@ -174,23 +198,54 @@ calculate_field_uniformity_kpi <- function(ci_pixels_by_field, field_boundaries_ #' #' @return Data frame with field-level CI changes calculate_area_change_kpi <- function(current_stats, previous_stats) { - result <- calculate_change_percentages(current_stats, previous_stats) - # Add interpretation - result$interpretation <- NA_character_ + # Initialize result data frame + result <- data.frame( + field_idx = seq_len(nrow(current_stats)), + mean_ci_pct_change = NA_real_, + interpretation = NA_character_, + stringsAsFactors = FALSE + ) - for (i in seq_len(nrow(result))) { - change <- result$mean_ci_pct_change[i] + if (is.null(previous_stats) || nrow(previous_stats) == 0) { + result$interpretation <- "No previous data" + return(result) + } + + # Match fields between current and previous stats + for (i in seq_len(nrow(current_stats))) { + field_id <- current_stats$Field_id[i] - if (is.na(change)) { + # Find matching field in previous stats + prev_idx <- which(previous_stats$Field_id == field_id) + + if (length(prev_idx) == 0) { result$interpretation[i] <- "No previous data" - } else if (change > 15) { + next + } + + prev_idx <- prev_idx[1] # Take first match + + current_ci <- current_stats$Mean_CI[i] + previous_ci <- previous_stats$Mean_CI[prev_idx] + + if (is.na(current_ci) || is.na(previous_ci) || previous_ci == 0) { + result$interpretation[i] <- "No previous data" + next + } + + # Calculate percentage change + pct_change <- ((current_ci - previous_ci) / previous_ci) * 100 + result$mean_ci_pct_change[i] <- round(pct_change, 2) + + # Add interpretation + if (pct_change > 15) { result$interpretation[i] <- "Rapid growth" - } else if (change > 5) { + } else if (pct_change > 5) { result$interpretation[i] <- "Positive growth" - } else if (change > -5) { + } else if (pct_change > -5) { result$interpretation[i] <- "Stable" - } else if (change > -15) { + } else if (pct_change > -15) { result$interpretation[i] <- "Declining" } else { result$interpretation[i] <- "Rapid decline" @@ -210,9 +265,37 @@ calculate_area_change_kpi <- function(current_stats, previous_stats) { #' #' @return Data frame with field-level TCH forecasts calculate_tch_forecasted_kpi <- function(field_statistics, harvesting_data = NULL, field_boundaries_sf = NULL) { + + # Handle both naming conventions (Field_id/Mean_CI vs field_idx/mean_ci) + if ("Field_id" %in% names(field_statistics)) { + # Add field_idx to match field_boundaries row numbers + field_statistics <- field_statistics %>% + mutate(field_idx = match(Field_id, field_boundaries_sf$field)) + mean_ci_col <- "Mean_CI" + } else { + mean_ci_col <- "mean_ci" + } + + # Filter out any fields without a match + field_statistics <- field_statistics %>% + filter(!is.na(field_idx)) + + if (nrow(field_statistics) == 0) { + warning("No fields matched between statistics and boundaries") + return(data.frame( + field_idx = integer(), + mean_ci = numeric(), + tch_forecasted = numeric(), + tch_lower_bound = numeric(), + tch_upper_bound = numeric(), + confidence = character(), + stringsAsFactors = FALSE + )) + } + result <- data.frame( field_idx = field_statistics$field_idx, - mean_ci = field_statistics$mean_ci, + mean_ci = field_statistics[[mean_ci_col]], tch_forecasted = NA_real_, tch_lower_bound = NA_real_, tch_upper_bound = NA_real_, @@ -360,73 +443,6 @@ calculate_weed_presence_kpi <- function(ci_pixels_by_field) { return(result) } -#' 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) { - # Handle both sf and SpatVector inputs - if (!inherits(field_boundaries, "SpatVector")) { - field_boundaries_vect <- terra::vect(field_boundaries) - } else { - field_boundaries_vect <- field_boundaries - } - - field_results <- data.frame() - - for (i in seq_len(nrow(field_boundaries))) { - field_name <- if ("field" %in% names(field_boundaries)) field_boundaries$field[i] else NA_character_ - sub_field_name <- if ("sub_field" %in% names(field_boundaries)) field_boundaries$sub_field[i] else NA_character_ - 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 <- round((low_ci_pixels / total_pixels) * 100, 2) - - # 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_ - )) - } - } - # 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)) -} # ============================================================================ # KPI ORCHESTRATOR AND REPORTING @@ -627,41 +643,135 @@ calculate_all_field_analysis_agronomic_support <- function( message("\n============ AURA KPI CALCULATION (6 KPIs) ============") - # Load current week mosaic - message("Loading current week mosaic...") - current_mosaic <- load_weekly_ci_mosaic(current_week, current_year, current_mosaic_dir) + # DETECT STRUCTURE FIRST - before any use of is_per_field + week_file <- sprintf("week_%02d_%d.tif", current_week, current_year) + field_dirs <- list.dirs(current_mosaic_dir, full.names = FALSE, recursive = FALSE) + field_dirs <- field_dirs[field_dirs != ""] - if (is.null(current_mosaic)) { - stop("Could not load current week mosaic") - } + is_per_field <- length(field_dirs) > 0 && + file.exists(file.path(current_mosaic_dir, field_dirs[1], week_file)) - # Extract field statistics - message("Extracting field statistics from current mosaic...") - current_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf) - #Extract CI pixels for each field individually - ci_pixels_by_field <- list() - for (i in seq_len(nrow(field_boundaries_sf))) { - field_vect <- terra::vect(field_boundaries_sf[i, ]) - ci_pixels_by_field[[i]] <- extract_ci_values(current_mosaic, field_vect) + if (is_per_field) { + message("Detected per-field mosaic structure...") + message("Using field-by-field extraction (similar to cane supply workflow)...") + + # Use the same extraction method as cane supply + current_stats <- calculate_field_statistics( + field_boundaries_sf, + current_week, + current_year, + current_mosaic_dir, + report_date = Sys.Date() + ) + + # Extract CI pixels for each field from their individual mosaics + ci_pixels_by_field <- list() + for (i in seq_len(nrow(field_boundaries_sf))) { + field_name <- field_boundaries_sf$field[i] + field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file) + + if (file.exists(field_mosaic_path)) { + tryCatch({ + field_raster <- terra::rast(field_mosaic_path) + ci_band <- field_raster[["CI"]] + field_vect <- terra::vect(field_boundaries_sf[i, ]) + ci_pixels_by_field[[i]] <- extract_ci_values(ci_band, field_vect) + }, error = function(e) { + message(paste(" Warning: Could not extract CI for field", field_name, ":", e$message)) + ci_pixels_by_field[[i]] <- NULL + }) + } else { + ci_pixels_by_field[[i]] <- NULL + } + } + + # For uniformity calculations that need a reference raster, load first available + current_mosaic <- NULL + for (field_name in field_dirs) { + field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file) + if (file.exists(field_mosaic_path)) { + tryCatch({ + current_mosaic <- terra::rast(field_mosaic_path)[["CI"]] + break + }, error = function(e) { + next + }) + } + } + + } else { + # Single-file mosaic (original behavior) + message("Loading current week mosaic...") + current_mosaic <- load_weekly_ci_mosaic(current_week, current_year, current_mosaic_dir) + + if (is.null(current_mosaic)) { + stop("Could not load current week mosaic") + } + + message("Extracting field statistics from current mosaic...") + current_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf) + + # Extract CI pixels for each field individually + ci_pixels_by_field <- list() + for (i in seq_len(nrow(field_boundaries_sf))) { + field_vect <- terra::vect(field_boundaries_sf[i, ]) + ci_pixels_by_field[[i]] <- extract_ci_values(current_mosaic, field_vect) + } } # Load previous week mosaic (if available) previous_stats <- NULL - if (!is.null(previous_mosaic_dir)) { + if (!is.null(previous_mosaic_dir) || is_per_field) { target_prev <- calculate_target_week_and_year(current_week, current_year, offset_weeks = 1) message(paste("Loading previous week mosaic (week", target_prev$week, target_prev$year, ")...")) - previous_mosaic <- load_weekly_ci_mosaic(target_prev$week, target_prev$year, previous_mosaic_dir) - if (!is.null(previous_mosaic)) { - previous_stats <- extract_field_statistics_from_ci(previous_mosaic, field_boundaries_sf) - } else { - message("Previous week mosaic not available - skipping area change KPI") + if (is_per_field) { + # Try loading previous week from the same directory structure + prev_week_file <- sprintf("week_%02d_%d.tif", target_prev$week, target_prev$year) + prev_field_exists <- any(sapply(field_dirs, function(field) { + file.exists(file.path(current_mosaic_dir, field, prev_week_file)) + })) + + if (prev_field_exists) { + message(" Found previous week per-field mosaics, calculating statistics...") + previous_stats <- calculate_field_statistics( + field_boundaries_sf, + target_prev$week, + target_prev$year, + current_mosaic_dir, + report_date = Sys.Date() - 7 + ) + } else { + message(" Previous week mosaic not available - skipping area change KPI") + } + } else if (!is.null(previous_mosaic_dir)) { + previous_mosaic <- load_weekly_ci_mosaic(target_prev$week, target_prev$year, previous_mosaic_dir) + + if (!is.null(previous_mosaic)) { + previous_stats <- extract_field_statistics_from_ci(previous_mosaic, field_boundaries_sf) + } else { + message(" Previous week mosaic not available - skipping area change KPI") + } } } # Calculate 6 KPIs message("\nCalculating KPI 1: Field Uniformity...") - uniformity_kpi <- calculate_field_uniformity_kpi(ci_pixels_by_field, field_boundaries_sf, current_mosaic) + if (is_per_field) { + uniformity_kpi <- calculate_field_uniformity_kpi( + ci_pixels_by_field, + field_boundaries_sf, + ci_band = NULL, + mosaic_dir = current_mosaic_dir, + week_file = week_file + ) + } else { + uniformity_kpi <- calculate_field_uniformity_kpi( + ci_pixels_by_field, + field_boundaries_sf, + current_mosaic + ) + } message("Calculating KPI 2: Area Change...") if (!is.null(previous_stats)) { @@ -679,19 +789,49 @@ calculate_all_field_analysis_agronomic_support <- function( message("Calculating KPI 4: Growth Decline...") growth_decline_kpi <- calculate_growth_decline_kpi( - list(ci_pixels_by_field) # Would need historical data for real trend + ci_pixels_by_field ) message("Calculating KPI 5: Weed Presence...") weed_kpi <- calculate_weed_presence_kpi(ci_pixels_by_field) message("Calculating KPI 6: Gap Filling...") - gap_filling_result <- calculate_gap_filling_kpi(current_mosaic, field_boundaries_sf) - - # Add field_idx to gap filling results - gap_filling_kpi <- gap_filling_result$field_results %>% - mutate(field_idx = row_number()) %>% - select(field_idx, gap_score, gap_level, mean_ci, outlier_threshold) + # Build list of per-field files for this week + per_field_files <- c() + for (field_name in field_dirs) { + field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file) + if (file.exists(field_mosaic_path)) { + per_field_files <- c(per_field_files, field_mosaic_path) + } + } + + if (length(per_field_files) > 0) { + # Use the common wrapper function (same as cane supply) + gap_scores_result <- calculate_gap_scores(per_field_files, field_boundaries_sf) + + # Convert to the format expected by orchestrator + gap_filling_kpi <- gap_scores_result %>% + mutate(field_idx = match(Field_id, field_boundaries_sf$field)) %>% + select(field_idx, gap_score) %>% + mutate( + gap_level = dplyr::case_when( + gap_score < 10 ~ "Minimal", + gap_score < 25 ~ "Moderate", + TRUE ~ "Significant" + ), + mean_ci = NA_real_, + outlier_threshold = NA_real_ + ) + } else { + # Fallback: no per-field files + gap_filling_kpi <- data.frame( + field_idx = seq_len(nrow(field_boundaries_sf)), + gap_score = NA_real_, + gap_level = NA_character_, + mean_ci = NA_real_, + outlier_threshold = NA_real_ + ) + } # Compile results all_kpis <- list(