#
# YOUNG_FIELD_ANALYSIS.R
# ======================
# Specialized analysis for young sugarcane fields (0-12 months)
# Focuses on:
# 1. Germination issues (gap detection, uneven emergence)
# 2. Weed pressure detection (both patchy and uniform weeds)
# 3. Age-specific intervention recommendations
#
# Usage: Rscript young_field_analysis.R [current_week] [previous_week] [project_dir]
#

# 1. Load required packages
# -----------------------
suppressPackageStartupMessages({
  library(sf)
  library(terra)
  library(tidyverse)
  library(lubridate)
  library(here)
  library(readxl)  # For reading harvest_data.xlsx
  library(spdep)   # For spatial statistics
})

# 2. Young field analysis configuration
# -----------------------------------
# Age thresholds (months)
GERMINATION_PHASE_MAX <- 4    # 0-4 months: critical germination period
WEED_CRITICAL_PHASE_MAX <- 8  # 4-8 months: weeds most competitive
YOUNG_FIELD_MAX <- 12         # 0-12 months: all young field analysis

# Detection thresholds
WEED_CI_CHANGE_THRESHOLD <- 1.5      # Weekly CI increase indicating possible weeds
SEVERE_WEED_CI_CHANGE <- 2.0         # Severe weed pressure
GERMINATION_GAP_CV_THRESHOLD <- 0.30  # High CV in young fields = poor germination
LOW_CI_GERMINATION_THRESHOLD <- 0.5   # Very low CI = poor emergence

# 3. Enhanced spectral indices for young crop detection
# ---------------------------------------------------

#' Calculate enhanced vegetation indices from RGBNIR bands
#' @param red Red band raster
#' @param green Green band raster  
#' @param blue Blue band raster
#' @param nir Near-infrared band raster
#' @return List of vegetation index rasters
calculate_enhanced_indices <- function(red, green, blue, nir) {
  
  cat("Calculating enhanced vegetation indices...\n")
  
  # 1. Standard NDVI (baseline)
  ndvi <- (nir - red) / (nir + red)
  names(ndvi) <- "NDVI"
  
  # 2. Green NDVI (more sensitive to early vegetation)
  gndvi <- (nir - green) / (nir + green)
  names(gndvi) <- "GNDVI"
  
  # 3. Excess Green Index (early vegetation detection)
  exg <- 2*green - red - blue
  names(exg) <- "ExG"
  
  # 4. Visible Atmospherically Resistant Index
  vari <- (green - red) / (green + red + blue)
  names(vari) <- "VARI"
  
  # 5. Green Ratio Vegetation Index
  grvi <- green / red
  names(grvi) <- "GRVI"
  
  # 6. Chlorophyll Index (CI = NIR / Green - 1, NOT NIR/Red)
  # *** CRITICAL: Correct formula uses GREEN band, not RED ***
  ci <- nir / green - 1
  names(ci) <- "CI"
  
  return(list(
    NDVI = ndvi,
    GNDVI = gndvi,
    ExG = exg,
    VARI = vari,
    GRVI = grvi,
    CI = ci
  ))
}

#' Create visual maps of different indices for comparison
#' @param indices_list List of index rasters
#' @param field_boundaries SF object with field boundaries
#' @param output_dir Directory to save maps
#' @param week_num Week number for filename
create_index_comparison_maps <- function(indices_list, field_boundaries, output_dir, week_num) {
  
  if (!dir.exists(output_dir)) {
    dir.create(output_dir, recursive = TRUE)
  }
  
  cat("Creating visual comparison maps and TIF exports...\n")
  
  # Create individual maps for each index
  for (index_name in names(indices_list)) {
    
    tryCatch({
      index_raster <- indices_list[[index_name]]
      
      # Create filenames
      map_filename <- file.path(output_dir, paste0("week_", sprintf("%02d", week_num), "_", index_name, "_map.png"))
      tif_filename <- file.path(output_dir, paste0("week_", sprintf("%02d", week_num), "_", index_name, ".tif"))
      
      # Export TIF file for QGIS
      terra::writeRaster(index_raster, tif_filename, overwrite = TRUE)
      cat("- Exported TIF:", tif_filename, "\n")
      
      # Create PNG map without field boundaries
      png(map_filename, width = 1200, height = 800, res = 150)
      
      # Plot raster only (no field boundaries)
      plot(index_raster, main = paste("Week", week_num, "-", index_name))
      
      dev.off()
      
      cat("- Saved map:", map_filename, "\n")
      
    }, error = function(e) {
      warning(paste("Error creating outputs for", index_name, ":", e$message))
    })
  }
  
  cat("✓ Index comparison maps and TIF files created\n")
}

#' Load field age data from existing harvest_data object (loaded in parameters_project.R)
#' @param harvest_data Data frame with harvest/field information
#' @param field_boundaries SF object with field boundaries
#' @return Field boundaries with age information added
#' # harvest_data = harvesting_data
#' # field_boundaries = field_boundaries_sf
load_field_ages_from_existing <- function(harvest_data, field_boundaries) {
  
  cat("Processing field age data from existing harvest_data object...\n")
  
  tryCatch({
    # Display structure of harvest data
    cat("- Harvest data columns:", paste(colnames(harvest_data), collapse = ", "), "\n")
    cat("- Harvest data rows:", nrow(harvest_data), "\n")
    
    # Display field boundaries structure
    cat("- Field boundaries columns:", paste(colnames(field_boundaries), collapse = ", "), "\n")
    cat("- Field boundaries rows:", nrow(field_boundaries), "\n")
    
    # Join harvest data with field boundaries on field and sub_field
    # Convert field and sub_field to character to ensure consistent matching
    # Filter for current season only (where season_end is today's date)
    current_date <- Sys.Date()
    
    harvest_data_clean <- harvest_data %>%
      mutate(
        field = as.character(field),
        sub_field = as.character(sub_field)
      ) %>%
      # Filter for current season only
      filter(season_end == current_date | (season_end >= current_date - 7 & season_end <= current_date)) %>%
      select(field, sub_field, age, season_start, season_end, tonnage_ha)
    
    cat("- Harvest data after filtering for current season:", nrow(harvest_data_clean), "rows\n")
    
    if (nrow(harvest_data_clean) == 0) {
      cat("⚠️ No current season data found in harvest data\n")
      cat("- Looking for season_end near:", current_date, "\n")
      
      # Show available season_end dates to help debug
      available_dates <- harvest_data %>%
        select(season_end) %>%
        distinct() %>%
        arrange(season_end)
      
      cat("- Available season_end dates in harvest data:\n")
      print(available_dates)
      
      return(field_boundaries)
    }
    
    field_boundaries_clean <- field_boundaries %>%
      mutate(
        field = as.character(field),
        sub_field = as.character(sub_field)
      )
    
    # Perform the join
    field_boundaries_with_age <- field_boundaries_clean %>%
      left_join(harvest_data_clean, by = c("field", "sub_field"))
    
    # Check join success
    matched_fields <- sum(!is.na(field_boundaries_with_age$age))
    total_fields <- nrow(field_boundaries_with_age)
    
    cat("✓ Successfully joined harvest data\n")
    cat("- Fields with age data:", matched_fields, "out of", total_fields, "\n")
    
    if (matched_fields > 0) {
      age_summary <- field_boundaries_with_age %>%
        filter(!is.na(age)) %>%
        pull(age) %>%
        summary()
      
      cat("- Age range (weeks):", paste(names(age_summary), "=", round(age_summary, 1), collapse = ", "), "\n")
      
      # Convert age from weeks to months for analysis
      field_boundaries_with_age <- field_boundaries_with_age %>%
        mutate(age_months = round(age / 4.33, 1))  # 4.33 weeks per month average
      
      cat("- Age range (months):", paste(round(range(field_boundaries_with_age$age_months, na.rm = TRUE), 1), collapse = " to "), "\n")
    }
    
    return(field_boundaries_with_age)
    
  }, error = function(e) {
    warning(paste("Error processing harvest data:", e$message))
    cat("Returning original field boundaries without age information\n")
    return(field_boundaries)
  })
}

#' Detect germination issues in young fields
#' @param indices_list List of vegetation index rasters
#' @param field_boundaries SF object with field boundaries (with age info)
#' @param young_fields_only Logical: analyze only young fields?
#' @return List with germination analysis results
detect_germination_issues <- function(indices_list, field_boundaries, young_fields_only = TRUE) {
  
  cat("=== GERMINATION ANALYSIS ===\n")
  
  germination_results <- list()
  field_boundaries_vect <- terra::vect(field_boundaries)
  
  for (i in seq_len(nrow(field_boundaries))) {
    field_name <- field_boundaries$field[i]
    sub_field_name <- field_boundaries$sub_field[i]
    field_id <- paste0(field_name, "_", sub_field_name)
    
    # Get field age if available
    field_age_months <- if ("age_months" %in% colnames(field_boundaries) && !is.na(field_boundaries$age_months[i])) {
      field_boundaries$age_months[i]
    } else {
      NA
    }
    
    cat("Analyzing germination for field:", field_id)
    if (!is.na(field_age_months)) {
      cat(" (Age:", field_age_months, "months)")
    }
    cat("\n")
    
    # Extract field boundary
    field_vect <- field_boundaries_vect[i]
    
    # Analyze each index for germination patterns
    field_analysis <- list()
    
    for (index_name in names(indices_list)) {
      index_raster <- indices_list[[index_name]]
      
      # Extract values for this field
      field_values <- terra::extract(index_raster, field_vect, fun = NULL)
      valid_values <- unlist(field_values)
      valid_values <- valid_values[!is.na(valid_values) & is.finite(valid_values)]
      
      if (length(valid_values) > 10) {
        # Calculate germination-relevant metrics
        mean_val <- mean(valid_values)
        cv_val <- sd(valid_values) / mean_val
        min_val <- min(valid_values)
        
        # Only check for germination issues in fields ≤ 4 months old
        check_germination <- is.na(field_age_months) || field_age_months <= GERMINATION_PHASE_MAX
        
        if (check_germination) {
          # Detect potential germination issues
          high_variation <- cv_val > GERMINATION_GAP_CV_THRESHOLD
          low_values <- mean_val < LOW_CI_GERMINATION_THRESHOLD
          very_low_areas <- sum(valid_values < (mean_val - 2*sd(valid_values))) / length(valid_values) * 100
        } else {
          # Field too old for germination analysis
          high_variation <- FALSE
          low_values <- FALSE
          very_low_areas <- 0
        }
        
        field_analysis[[index_name]] <- list(
          mean = mean_val,
          cv = cv_val,
          min = min_val,
          high_variation = high_variation,
          low_values = low_values,
          very_low_areas_pct = very_low_areas,
          n_pixels = length(valid_values),
          age_months = field_age_months,
          germination_analysis_applied = check_germination
        )
      }
    }
    
    germination_results[[field_id]] <- field_analysis
  }
  
  return(germination_results)
}

#' Detect weed pressure using change analysis
#' @param current_indices List of current week indices
#' @param previous_indices List of previous week indices  
#' @param field_boundaries SF object with field boundaries
#' @return List with weed detection results
detect_weed_pressure <- function(current_indices, previous_indices, field_boundaries) {
  
  cat("=== WEED PRESSURE DETECTION ===\n")
  
  weed_results <- list()
  field_boundaries_vect <- terra::vect(field_boundaries)
  
  for (i in seq_len(nrow(field_boundaries))) {
    field_name <- field_boundaries$field[i]
    sub_field_name <- field_boundaries$sub_field[i]
    field_id <- paste0(field_name, "_", sub_field_name)
    
    cat("Analyzing weed pressure for field:", field_id, "\n")
    
    field_vect <- field_boundaries_vect[i]
    field_weed_analysis <- list()
    
    # Analyze change in each index
    for (index_name in names(current_indices)) {
      
      if (index_name %in% names(previous_indices)) {
        
        current_raster <- current_indices[[index_name]]
        previous_raster <- previous_indices[[index_name]]
        
        # Extract values for both weeks
        current_values <- unlist(terra::extract(current_raster, field_vect, fun = NULL))
        previous_values <- unlist(terra::extract(previous_raster, field_vect, fun = NULL))
        
        # Clean values
        valid_idx <- !is.na(current_values) & !is.na(previous_values) & 
                     is.finite(current_values) & is.finite(previous_values)
        current_clean <- current_values[valid_idx]
        previous_clean <- previous_values[valid_idx]
        
        if (length(current_clean) > 10) {
          
          # Calculate change metrics
          change_values <- current_clean - previous_clean
          mean_change <- mean(change_values)
          change_cv <- sd(change_values) / abs(mean(current_clean))
          
          # Weed detection criteria
          # 1. Significant positive change (weeds growing faster than cane)
          rapid_increase <- mean_change >= WEED_CI_CHANGE_THRESHOLD
          severe_increase <- mean_change >= SEVERE_WEED_CI_CHANGE
          
          # 2. Percentage of field with rapid increase
          rapid_increase_pct <- sum(change_values >= WEED_CI_CHANGE_THRESHOLD) / length(change_values) * 100
          
          # 3. Patchy vs uniform weed patterns
          patchy_weeds <- change_cv > 0.5 && rapid_increase_pct > 10 && rapid_increase_pct < 80
          uniform_weeds <- change_cv < 0.3 && rapid_increase_pct > 60  # Whole field weeds
          
          # 4. Current vegetation level (high CI + rapid change = likely weeds)
          current_mean <- mean(current_clean)
          high_current_ci <- current_mean > 2.0  # Adjust threshold as needed
          
          field_weed_analysis[[index_name]] <- list(
            mean_change = mean_change,
            change_cv = change_cv,
            current_mean = current_mean,
            rapid_increase = rapid_increase,
            severe_increase = severe_increase,
            rapid_increase_pct = rapid_increase_pct,
            patchy_weeds = patchy_weeds,
            uniform_weeds = uniform_weeds,
            high_current_ci = high_current_ci,
            n_pixels = length(current_clean)
          )
        }
      }
    }
    
    weed_results[[field_id]] <- field_weed_analysis
  }
  
  return(weed_results)
}

#' Generate intervention recommendations based on detected issues
#' @param field_id Field identifier
#' @param germination_analysis Germination analysis results
#' @param weed_analysis Weed detection results
#' @param field_age_months Field age in months (if available)
#' @return List with recommendations
generate_young_field_recommendations <- function(field_id, germination_analysis, weed_analysis, field_age_months = NA) {
  
  recommendations <- list()
  
  # Priority level: 1=urgent, 2=high, 3=moderate, 4=monitoring
  priority <- 4
  messages <- c()
  interventions <- c()
  
  # Check for germination issues (if CI analysis available and field is young enough)
  if ("CI" %in% names(germination_analysis)) {
    ci_analysis <- germination_analysis[["CI"]]
    
    # Only apply germination analysis for fields ≤ 4 months old
    if (!is.null(ci_analysis$germination_analysis_applied) && ci_analysis$germination_analysis_applied) {
      
      if (ci_analysis$high_variation && ci_analysis$low_values) {
        priority <- min(priority, 1)
        messages <- c(messages, "🚨 URGENT: Poor germination detected - high variation with low CI values")
        interventions <- c(interventions, "Immediate replanting in gap areas", "Check seed quality and planting conditions")
      } else if (ci_analysis$high_variation) {
        priority <- min(priority, 2)
        messages <- c(messages, "⚠️ Alert: Uneven emergence detected - moderate intervention needed")
        interventions <- c(interventions, "Monitor gap areas for potential replanting")
      }
      
      if (ci_analysis$very_low_areas_pct > 15) {
        priority <- min(priority, 2)
        messages <- c(messages, paste0("⚠️ Alert: ", round(ci_analysis$very_low_areas_pct, 1), "% of field has very low vegetation"))
        interventions <- c(interventions, "Investigate cause of low-performing areas")
      }
      
    } else {
      # Field too old for germination analysis - check for other issues
      if (ci_analysis$cv > 0.5) {
        priority <- min(priority, 3)
        messages <- c(messages, "📊 Info: Uneven crop growth detected (field too old for replanting)")
        interventions <- c(interventions, "Monitor for potential management issues")
      }
    }
  }
  
  # Check for weed pressure (if CI analysis available)
  if ("CI" %in% names(weed_analysis)) {
    ci_weed <- weed_analysis[["CI"]]
    
    if (ci_weed$severe_increase) {
      priority <- min(priority, 1)
      messages <- c(messages, paste0("🚨 CRITICAL: Severe weed pressure detected - CI increased by ", round(ci_weed$mean_change, 2)))
      interventions <- c(interventions, "Immediate weed control required", "Consider herbicide application")
    } else if (ci_weed$rapid_increase) {
      priority <- min(priority, 2)
      messages <- c(messages, paste0("⚠️ Alert: Weed pressure detected - CI increased by ", round(ci_weed$mean_change, 2)))
      interventions <- c(interventions, "Schedule weed control within 1-2 weeks")
    }
    
    if (ci_weed$uniform_weeds) {
      priority <- min(priority, 2)
      messages <- c(messages, paste0("🔶 Pattern: Uniform weed emergence across ", round(ci_weed$rapid_increase_pct, 1), "% of field"))
      interventions <- c(interventions, "Broad-scale weed management needed")
    } else if (ci_weed$patchy_weeds) {
      priority <- min(priority, 3)
      messages <- c(messages, paste0("🔶 Pattern: Patchy weed emergence in ", round(ci_weed$rapid_increase_pct, 1), "% of field"))
      interventions <- c(interventions, "Targeted spot treatment recommended")
    }
  }
  
  # Age-specific recommendations
  if (!is.na(field_age_months)) {
    if (field_age_months <= GERMINATION_PHASE_MAX) {
      interventions <- c(interventions, "Critical germination phase - maintain optimal moisture", "Monitor for pest damage")
    } else if (field_age_months <= WEED_CRITICAL_PHASE_MAX) {
      interventions <- c(interventions, "Peak weed competition period - prioritize weed control")
    }
  }
  
  # Default monitoring if no issues
  if (length(messages) == 0) {
    messages <- c("✅ No significant issues detected")
    interventions <- c("Continue routine monitoring")
  }
  
  return(list(
    field_id = field_id,
    priority = priority,
    messages = messages,
    interventions = interventions,
    requires_action = priority <= 3
  ))
}

# 4. Main analysis function
# -----------------------
main_young_field_analysis <- function() {
  
  # Get command line arguments
  args <- commandArgs(trailingOnly = TRUE)
  
  current_week <- if (length(args) >= 1) as.numeric(args[1]) else 30
  previous_week <- if (length(args) >= 2) as.numeric(args[2]) else 29
  project_dir <- if (length(args) >= 3) args[3] else "simba"
  
  cat("=== YOUNG FIELD ANALYSIS SYSTEM ===\n")
  cat("Project:", project_dir, "\n")
  cat("Analyzing weeks:", previous_week, "→", current_week, "\n\n")
  
  # Load project configuration
  assign("project_dir", project_dir, envir = .GlobalEnv)
  
  tryCatch({
    source("parameters_project.R")
    cat("✓ Project configuration loaded\n")
  }, error = function(e) {
    tryCatch({
      source(here::here("r_app", "parameters_project.R"))
      cat("✓ Project configuration loaded from r_app directory\n")
    }, error = function(e) {
      stop("Failed to load project configuration")
    })
  })
  
  # Verify required variables
  if (!exists("weekly_CI_mosaic") || !exists("field_boundaries_sf")) {
    stop("Required project variables not found")
  }
  
  # Check if harvest data is available from parameters_project.R
  if (exists("harvesting_data") && !is.null(harvesting_data)) {
    cat("✓ Harvest data loaded from parameters_project.R\n")
    field_boundaries_with_age <- load_field_ages_from_existing(harvesting_data, field_boundaries_sf)
  } else {
    cat("⚠️ No harvest data found in parameters_project.R\n")
    field_boundaries_with_age <- field_boundaries_sf
  }
  
  # Construct mosaic file paths
  year <- 2025
  current_week_file <- sprintf("week_%02d_%d.tif", current_week, year)
  previous_week_file <- sprintf("week_%02d_%d.tif", previous_week, year)
  
  current_mosaic_path <- file.path(weekly_CI_mosaic, current_week_file)
  previous_mosaic_path <- file.path(weekly_CI_mosaic, previous_week_file)
  
  cat("Current week mosaic:", current_mosaic_path, "\n")
  cat("Previous week mosaic:", previous_mosaic_path, "\n")
  
  # Check if files exist
  if (!file.exists(current_mosaic_path)) {
    stop("Current week mosaic not found: ", current_mosaic_path)
  }
  if (!file.exists(previous_mosaic_path)) {
    stop("Previous week mosaic not found: ", previous_mosaic_path)
  }
  
  # Load mosaics and extract bands
  cat("\nLoading mosaics and extracting RGBNIR bands...\n")
  
  current_mosaic <- terra::rast(current_mosaic_path)
  previous_mosaic <- terra::rast(previous_mosaic_path)
  
  cat("Current mosaic bands:", nlyr(current_mosaic), "\n")
  cat("Previous mosaic bands:", nlyr(previous_mosaic), "\n")
  
  # Extract RGBNIR bands (assuming standard order: R=1, G=2, B=3, NIR=4)
  current_red <- current_mosaic[[1]]
  current_green <- current_mosaic[[2]]
  current_blue <- current_mosaic[[3]]
  current_nir <- current_mosaic[[4]]
  
  previous_red <- previous_mosaic[[1]]
  previous_green <- previous_mosaic[[2]]
  previous_blue <- previous_mosaic[[3]]
  previous_nir <- previous_mosaic[[4]]
  
  # Calculate enhanced indices for both weeks
  cat("\nCalculating vegetation indices...\n")
  current_indices <- calculate_enhanced_indices(current_red, current_green, current_blue, current_nir)
  previous_indices <- calculate_enhanced_indices(previous_red, previous_green, previous_blue, previous_nir)
  
  # Create output directory for maps
  output_dir <- file.path(dirname(weekly_CI_mosaic), "young_field_analysis")
  
  # Create visual comparison maps
  cat("\nCreating visual maps...\n")
# Check if maps already exist before creating them
current_week_maps_exist <- all(sapply(names(current_indices), function(index_name) {
    map_file <- file.path(output_dir, paste0("week_", sprintf("%02d", current_week), "_", index_name, "_map.png"))
    tif_file <- file.path(output_dir, paste0("week_", sprintf("%02d", current_week), "_", index_name, ".tif"))
    file.exists(map_file) && file.exists(tif_file)
}))

previous_week_maps_exist <- all(sapply(names(previous_indices), function(index_name) {
    map_file <- file.path(output_dir, paste0("week_", sprintf("%02d", previous_week), "_", index_name, "_map.png"))
    tif_file <- file.path(output_dir, paste0("week_", sprintf("%02d", previous_week), "_", index_name, ".tif"))
    file.exists(map_file) && file.exists(tif_file)
}))

if (!current_week_maps_exist) {
    cat("Creating current week maps (week", current_week, ")...\n")
    create_index_comparison_maps(current_indices, field_boundaries_sf, output_dir, current_week)
} else {
    cat("Current week maps already exist, skipping creation...\n")
}

if (!previous_week_maps_exist) {
    cat("Creating previous week maps (week", previous_week, ")...\n")
    create_index_comparison_maps(previous_indices, field_boundaries_sf, output_dir, previous_week)
} else {
    cat("Previous week maps already exist, skipping creation...\n")
}
  
  # Perform germination analysis
  cat("\nPerforming germination analysis...\n")
  germination_results <- detect_germination_issues(current_indices, field_boundaries_with_age)
  
  # Perform weed detection analysis
  cat("\nPerforming weed pressure analysis...\n")
  weed_results <- detect_weed_pressure(current_indices, previous_indices, field_boundaries_with_age)
  
  # Generate recommendations for each field
  cat("\n=== FIELD RECOMMENDATIONS ===\n")
  
  field_recommendations <- list()
  action_needed_fields <- 0
  
  for (field_id in names(germination_results)) {
    
    germination_analysis <- germination_results[[field_id]]
    weed_analysis <- if (field_id %in% names(weed_results)) weed_results[[field_id]] else list()
    
    # Get field age if available
    field_row <- field_boundaries_with_age[field_boundaries_with_age$field == strsplit(field_id, "_")[[1]][1] & 
                                          field_boundaries_with_age$sub_field == strsplit(field_id, "_")[[1]][2], ]
    field_age_months <- if (nrow(field_row) > 0 && !is.na(field_row$age_months[1])) field_row$age_months[1] else NA
    
    recommendations <- generate_young_field_recommendations(
      field_id, 
      germination_analysis, 
      weed_analysis,
      field_age_months
    )
    
    field_recommendations[[field_id]] <- recommendations
    
    # Print recommendations
    cat("\nFIELD:", field_id, "\n")
    if (!is.na(field_age_months)) {
      cat("Age:", field_age_months, "months\n")
    }
    cat("Priority Level:", recommendations$priority, "\n")
    
    for (message in recommendations$messages) {
      cat("-", message, "\n")
    }
    
    if (length(recommendations$interventions) > 0) {
      cat("Recommended Actions:\n")
      for (intervention in recommendations$interventions) {
        cat("  •", intervention, "\n")
      }
    }
    
    if (recommendations$requires_action) {
      action_needed_fields <- action_needed_fields + 1
    }
    
    # Show index comparison for this field
    if ("CI" %in% names(germination_analysis)) {
      ci_stats <- germination_analysis[["CI"]]
      cat("CI Stats: Mean =", round(ci_stats$mean, 3), 
          ", CV =", round(ci_stats$cv, 3),
          ", Low areas =", round(ci_stats$very_low_areas_pct, 1), "%\n")
    }
    
    if ("CI" %in% names(weed_analysis)) {
      weed_stats <- weed_analysis[["CI"]]
      cat("Change Stats: CI Δ =", round(weed_stats$mean_change, 3),
          ", Rapid increase =", round(weed_stats$rapid_increase_pct, 1), "%\n")
    }
  }
  
  # Summary
  cat("\n=== ANALYSIS SUMMARY ===\n")
  cat("Total fields analyzed:", length(field_recommendations), "\n")
  cat("Fields requiring action:", action_needed_fields, "\n")
  cat("Maps saved to:", output_dir, "\n")
  
  cat("\nIndex comparison maps created for visual inspection:\n")
  cat("- NDVI: Standard vegetation index\n")
  cat("- GNDVI: Green NDVI (sensitive to early vegetation)\n")
  cat("- ExG: Excess Green (early detection)\n")
  cat("- VARI: Atmospherically resistant\n")
  cat("- GRVI: Green ratio\n")
  cat("- CI: Chlorophyll Index (current standard)\n")
  
  cat("\n✓ Young field analysis complete\n")
  
  return(list(
    germination_results = germination_results,
    weed_results = weed_results,
    recommendations = field_recommendations,
    indices_calculated = names(current_indices),
    maps_created = file.path(output_dir)
  ))
}

# Run analysis if script called directly
if (sys.nframe() == 0) {
  main_young_field_analysis()
}