1. removed gap score 2.fixed morans I to calculate per-field

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
+  # Initialize result data frame
-  result$interpretation <- NA_character_
+  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))) {
+  if (is.null(previous_stats) || nrow(previous_stats) == 0) {
-    change <- result$mean_ci_pct_change[i]
+    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
+  # DETECT STRUCTURE FIRST - before any use of is_per_field
-  message("Loading current week mosaic...")
+  week_file <- sprintf("week_%02d_%d.tif", current_week, current_year)
-  current_mosaic <- load_weekly_ci_mosaic(current_week, current_year, current_mosaic_dir)
+  field_dirs <- list.dirs(current_mosaic_dir, full.names = FALSE, recursive = FALSE)
+  field_dirs <- field_dirs[field_dirs != ""]
   
-  if (is.null(current_mosaic)) {
+  is_per_field <- length(field_dirs) > 0 && 
-    stop("Could not load current week mosaic")
+                  file.exists(file.path(current_mosaic_dir, field_dirs[1], week_file))
-  }
   
-  # Extract field statistics
+  if (is_per_field) {
-  message("Extracting field statistics from current mosaic...")
+    message("Detected per-field mosaic structure...")
-  current_stats <- extract_field_statistics_from_ci(current_mosaic, field_boundaries_sf)
+    message("Using field-by-field extraction (similar to cane supply workflow)...")
-  #Extract CI pixels for each field individually
+    
-  ci_pixels_by_field <- list()
+    # Use the same extraction method as cane supply
-  for (i in seq_len(nrow(field_boundaries_sf))) {
+    current_stats <- calculate_field_statistics(
-    field_vect <- terra::vect(field_boundaries_sf[i, ])
+      field_boundaries_sf, 
-    ci_pixels_by_field[[i]] <- extract_ci_values(current_mosaic, field_vect)
+      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)) {
+    if (is_per_field) {
-      previous_stats <- extract_field_statistics_from_ci(previous_mosaic, field_boundaries_sf)
+      # Try loading previous week from the same directory structure
-    } else {
+      prev_week_file <- sprintf("week_%02d_%d.tif", target_prev$week, target_prev$year)
-      message("Previous week mosaic not available - skipping area change KPI")
+      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)
+  # Build list of per-field files for this week
-
+  per_field_files <- c()
-  # Add field_idx to gap filling results
+  for (field_name in field_dirs) {
-  gap_filling_kpi <- gap_filling_result$field_results %>%
+    field_mosaic_path <- file.path(current_mosaic_dir, field_name, week_file)
-  mutate(field_idx = row_number()) %>%
+    if (file.exists(field_mosaic_path)) {
-  select(field_idx, gap_score, gap_level, mean_ci, outlier_threshold)
+      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(