# Enhanced NDVI with quality filtering
ndvi_enhanced <- calculate_ndvi_enhanced(
  red_data = red_band,
  nir_data = nir_band,
  quality_filter = TRUE,
  temporal_smoothing = FALSE,  # Single date
  verbose = TRUE
)

# Visualize enhanced NDVI
create_ndvi_map(
  ndvi_data = ndvi_enhanced,
  title = "Enhanced NDVI with Quality Filtering",
  ndvi_classes = "none"
)

print("Enhanced NDVI Analysis:")
ndvi_values <- terra::values(ndvi_enhanced, mat = FALSE)
valid_ndvi <- ndvi_values[!is.na(ndvi_values)]
print(paste("Mean NDVI:", round(mean(valid_ndvi), 3)))
print(paste("NDVI range:", paste(round(range(valid_ndvi), 3), collapse = " - ")))

# Calculate multiple indices relevant to agriculture
agricultural_indices <- calculate_multiple_indices(
  red = red_band,
  nir = nir_band,
  indices = c("NDVI", "SAVI", "MSAVI", "DVI", "RVI"),
  output_stack = TRUE,
  verbose = TRUE
)

print("Agricultural Indices Calculated:")
print(names(agricultural_indices))

# Visualize key indices
if (terra::nlyr(agricultural_indices) >= 2) {
  # NDVI visualization
  create_spatial_map(
    spatial_data = agricultural_indices[[1]],
    title = "NDVI for Agricultural Assessment",
    color_scheme = "ndvi"
  )
  
  # SAVI visualization (soil-adjusted)
  if ("SAVI" %in% names(agricultural_indices)) {
    plot_raster_fast(
      agricultural_indices[["SAVI"]],
      title = "SAVI (Soil Adjusted Vegetation Index)",
      color_scheme = "viridis"
    )
  }
}

# Extract yield analysis if available
if (!is.null(crop_analysis$analysis_results$yield_analysis)) {
  yield_results <- crop_analysis$analysis_results$yield_analysis
  
  print("Yield Potential Analysis:")
  if (!"error" %in% names(yield_results)) {
    print(paste("Composite yield index:", round(yield_results$composite_yield_index, 3)))
    print(paste("Yield potential class:", yield_results$yield_potential_class))
    print(paste("Classification confidence:", round(yield_results$classification_confidence, 3)))
    
    print("Index contributions to yield assessment:")
    for (idx in names(yield_results$index_contributions)) {
      contrib <- yield_results$index_contributions[[idx]]
      print(paste("  ", idx, "- Normalized:", round(contrib$mean_normalized, 3),
                  "Raw mean:", round(contrib$raw_mean, 3)))
    }
  }
}

# Create yield potential map visualization
yield_potential_raster <- agricultural_indices[[1]]  # Use NDVI as proxy
terra::values(yield_potential_raster) <- ifelse(terra::values(yield_potential_raster) > 0.6, 3,
                                               ifelse(terra::values(yield_potential_raster) > 0.4, 2, 1))
names(yield_potential_raster) <- "Yield_Potential"

plot_raster_fast(
  yield_potential_raster,
  title = "Yield Potential Classification",
  color_scheme = "terrain"
)

# Calculate comprehensive vegetation statistics
if (inherits(agricultural_indices, "SpatRaster")) {
  veg_stats <- list()
  
  for (i in 1:terra::nlyr(agricultural_indices)) {
    index_name <- names(agricultural_indices)[i]
    values <- terra::values(agricultural_indices[[i]], mat = FALSE)
    valid_values <- values[!is.na(values)]
    
    if (length(valid_values) > 0) {
      veg_stats[[index_name]] <- list(
        mean = mean(valid_values),
        median = median(valid_values),
        sd = sd(valid_values),
        cv = sd(valid_values) / mean(valid_values),
        min = min(valid_values),
        max = max(valid_values),
        q25 = quantile(valid_values, 0.25),
        q75 = quantile(valid_values, 0.75)
      )
    }
  }
  
  print("Crop Performance Metrics:")
  for (index in names(veg_stats)) {
    stats <- veg_stats[[index]]
    print(paste(index, "- Mean:", round(stats$mean, 3), 
                "CV:", round(stats$cv, 3),
                "Range:", paste(round(c(stats$min, stats$max), 3), collapse = "-")))
  }
}

# Create sample field boundaries
field_1 <- sf::st_polygon(list(matrix(c(
  -83.5, 40.0, -83.3, 40.0, -83.3, 40.2, -83.5, 40.2, -83.5, 40.0
), ncol = 2, byrow = TRUE)))

field_2 <- sf::st_polygon(list(matrix(c(
  -83.3, 40.0, -83.1, 40.0, -83.1, 40.2, -83.3, 40.2, -83.3, 40.0
), ncol = 2, byrow = TRUE)))

fields_sf <- sf::st_sf(
  field_id = c("Field_1", "Field_2"),
  crop_type = c("Corn", "Soybeans"),
  geometry = sf::st_sfc(field_1, field_2, crs = 4326)
)

# Extract vegetation statistics by field using spatial join
field_analysis <- universal_spatial_join(
  source_data = agricultural_indices,
  target_data = fields_sf,
  method = "zonal",
  summary_function = "mean",
  verbose = TRUE
)

print("Field-Level Analysis Results:")
if (inherits(field_analysis, "sf")) {
  print(sf::st_drop_geometry(field_analysis))
}

# Calculate coefficient of variation for management zones
if (inherits(agricultural_indices, "SpatRaster") && terra::nlyr(agricultural_indices) > 0) {
  # Use NDVI for variability analysis
  ndvi_layer <- agricultural_indices[[1]]
  
  # Calculate local variability using focal statistics
  local_cv <- terra::focal(ndvi_layer, w = matrix(1, 3, 3), 
                           fun = function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE))
  names(local_cv) <- "Local_CV"
  
  # Classify variability zones
  variability_zones <- terra::classify(local_cv, 
                                       matrix(c(0, 0.1, 1,     # Low variability
                                               0.1, 0.2, 2,    # Medium variability  
                                               0.2, 1, 3),     # High variability
                                              ncol = 3, byrow = TRUE))
  names(variability_zones) <- "Variability_Zones"
  
  # Visualize variability
  plot_raster_fast(
    local_cv,
    title = "Spatial Variability (Coefficient of Variation)",
    color_scheme = "plasma"
  )
  
  plot_raster_fast(
    variability_zones,
    title = "Management Zones (1=Low, 2=Medium, 3=High Variability)",
    color_scheme = "categorical"
  )
}

# Create stress detection workflow
stress_indices <- c("NDVI", "SAVI", "DVI")

# Calculate stress thresholds based on crop type
corn_thresholds <- list(
  healthy = list(NDVI = c(0.6, 1.0), SAVI = c(0.4, 0.8)),
  stressed = list(NDVI = c(0.3, 0.6), SAVI = c(0.2, 0.4)),
  severely_stressed = list(NDVI = c(0.0, 0.3), SAVI = c(0.0, 0.2))
)

# Apply stress classification
if (inherits(agricultural_indices, "SpatRaster") && "NDVI" %in% names(agricultural_indices)) {
  ndvi_values <- terra::values(agricultural_indices[["NDVI"]], mat = FALSE)
  
  # Classify stress levels
  stress_classification <- ifelse(ndvi_values >= 0.6, "Healthy",
                                 ifelse(ndvi_values >= 0.3, "Moderate Stress", "Severe Stress"))
  
  # Create stress map
  stress_raster <- agricultural_indices[["NDVI"]]
  stress_numeric <- ifelse(ndvi_values >= 0.6, 3,
                          ifelse(ndvi_values >= 0.3, 2, 1))
  terra::values(stress_raster) <- stress_numeric
  names(stress_raster) <- "Stress_Level"
  
  plot_raster_fast(
    stress_raster,
    title = "Crop Stress Classification (1=Severe, 2=Moderate, 3=Healthy)",
    color_scheme = "terrain"
  )
  
  # Calculate stress statistics
  stress_table <- table(stress_classification)
  stress_percent <- round(prop.table(stress_table) * 100, 1)
  
  print("Crop Stress Distribution:")
  for (level in names(stress_table)) {
    print(paste(level, ":", stress_table[[level]], "pixels (", stress_percent[[level]], "%)"))
  }
}

# Define warning thresholds
warning_thresholds <- list(
  drought_stress = 0.4,    # NDVI below this indicates drought stress
  disease_risk = 0.3,      # NDVI below this indicates disease risk
  optimal_growth = 0.7     # NDVI above this indicates optimal conditions
)

# Generate alerts
if (exists("ndvi_values")) {
  alerts <- list()
  
  # Calculate percentages in each category
  drought_pixels <- sum(ndvi_values < warning_thresholds$drought_stress, na.rm = TRUE)
  disease_pixels <- sum(ndvi_values < warning_thresholds$disease_risk, na.rm = TRUE)
  optimal_pixels <- sum(ndvi_values > warning_thresholds$optimal_growth, na.rm = TRUE)
  total_pixels <- sum(!is.na(ndvi_values))
  
  alerts$drought_risk <- round((drought_pixels / total_pixels) * 100, 1)
  alerts$disease_risk <- round((disease_pixels / total_pixels) * 100, 1)
  alerts$optimal_conditions <- round((optimal_pixels / total_pixels) * 100, 1)
  
  print("Agricultural Alert System:")
  print(paste("Drought stress risk:", alerts$drought_risk, "% of field"))
  print(paste("Disease risk:", alerts$disease_risk, "% of field"))
  print(paste("Optimal conditions:", alerts$optimal_conditions, "% of field"))
  
  # Generate recommendations
  if (alerts$drought_risk > 20) {
    print("RECOMMENDATION: Consider irrigation scheduling")
  }
  if (alerts$disease_risk > 10) {
    print("RECOMMENDATION: Increase disease monitoring")
  }
  if (alerts$optimal_conditions > 70) {
    print("STATUS: Crop conditions are generally favorable")
  }
}

# Simulate time series data for growing season
growth_stages <- c("Planting", "V6", "V12", "R1", "R3", "R6")
ndvi_progression <- c(0.2, 0.4, 0.7, 0.8, 0.75, 0.6)

# Create time series visualization concept
seasonal_data <- data.frame(
  Stage = growth_stages,
  NDVI = ndvi_progression,
  DOY = c(120, 150, 180, 200, 220, 260)  # Day of year
)

print("Seasonal NDVI Progression:")
print(seasonal_data)

# Create conceptual growth curve
if (requireNamespace("ggplot2", quietly = TRUE)) {
  library(ggplot2)
  
  growth_plot <- ggplot(seasonal_data, aes(x = DOY, y = NDVI)) +
    geom_line(color = "darkgreen", size = 1.2) +
    geom_point(color = "red", size = 3) +
    geom_text(aes(label = Stage), vjust = -0.5) +
    labs(title = "Typical Corn NDVI Progression",
         x = "Day of Year",
         y = "NDVI Value") +
    theme_minimal() +
    ylim(0, 1)
  
  print(growth_plot)
}

# Define maturity indicators based on NDVI decline
maturity_thresholds <- list(
  corn = list(
    early_maturity = 0.7,     # NDVI starts declining
    harvest_ready = 0.5,      # Optimal harvest window
    post_harvest = 0.3        # Past optimal timing
  ),
  soybeans = list(
    early_maturity = 0.6,
    harvest_ready = 0.4,
    post_harvest = 0.25
  )
)

# Calculate harvest readiness
if (exists("ndvi_values")) {
  # Use corn thresholds for demonstration
  thresholds <- maturity_thresholds$corn
  
  harvest_assessment <- list(
    early_maturity = sum(ndvi_values <= thresholds$early_maturity & 
                        ndvi_values > thresholds$harvest_ready, na.rm = TRUE),
    harvest_ready = sum(ndvi_values <= thresholds$harvest_ready & 
                       ndvi_values > thresholds$post_harvest, na.rm = TRUE),
    post_optimal = sum(ndvi_values <= thresholds$post_harvest, na.rm = TRUE),
    still_growing = sum(ndvi_values > thresholds$early_maturity, na.rm = TRUE)
  )
  
  total_pixels <- sum(!is.na(ndvi_values))
  
  print("Harvest Timing Assessment:")
  for (stage in names(harvest_assessment)) {
    percentage <- round((harvest_assessment[[stage]] / total_pixels) * 100, 1)
    print(paste(stage, ":", harvest_assessment[[stage]], "pixels (", percentage, "%)"))
  }
}

# Calculate water-related indices
green_band <- red_band
terra::values(green_band) <- runif(2500, 0.1, 0.2)
names(green_band) <- "Green"

# Calculate NDWI for water stress detection
ndwi <- calculate_water_index(
  green = green_band,
  nir = nir_band,
  index_type = "NDWI",
  verbose = TRUE
)

# Irrigation needs based on NDWI
ndwi_values <- terra::values(ndwi, mat = FALSE)
valid_ndwi <- ndwi_values[!is.na(ndwi_values)]

if (length(valid_ndwi) > 0) {
  irrigation_zones <- ifelse(valid_ndwi < -0.3, "High Need",
                           ifelse(valid_ndwi < -0.1, "Moderate Need", "Adequate"))
  
  irrigation_table <- table(irrigation_zones)
  irrigation_percent <- round(prop.table(irrigation_table) * 100, 1)
  
  print("Irrigation Needs Assessment:")
  for (zone in names(irrigation_table)) {
    print(paste(zone, ":", irrigation_table[[zone]], "pixels (", irrigation_percent[[zone]], "%)"))
  }
  
  # Visualize irrigation zones
  irrigation_raster <- ndwi
  terra::values(irrigation_raster) <- as.numeric(as.factor(irrigation_zones))
  names(irrigation_raster) <- "Irrigation_Needs"
  
  plot_raster_fast(
    irrigation_raster,
    title = "Irrigation Needs Assessment (1=Adequate, 2=High, 3=Moderate)",
    color_scheme = "water"
  )
}

# Simulate multi-year CDL data
create_rotation_analysis <- function() {
  # Create sample rotation data
  rotation_data <- data.frame(
    year = rep(2021:2023, each = 4),
    field = rep(c("Field_A", "Field_B", "Field_C", "Field_D"), 3),
    crop = c(
      "Corn", "Soybeans", "Corn", "Wheat",      # 2021
      "Soybeans", "Corn", "Soybeans", "Corn",   # 2022
      "Corn", "Soybeans", "Corn", "Soybeans"   # 2023
    ),
    yield = runif(12, 80, 200)
  )
  
  return(rotation_data)
}

rotation_analysis <- create_rotation_analysis()
print("Crop Rotation Analysis:")
print(rotation_analysis)

# Analyze rotation patterns
rotation_patterns <- list()
fields <- unique(rotation_analysis$field)

for (field in fields) {
  field_data <- rotation_analysis[rotation_analysis$field == field, ]
  rotation_patterns[[field]] <- paste(field_data$crop, collapse = " → ")
}

print("Rotation Patterns:")
for (field in names(rotation_patterns)) {
  print(paste(field, ":", rotation_patterns[[field]]))
}

# Calculate diversity index for crop rotation
calculate_crop_diversity <- function(rotation_data) {
  diversity_scores <- list()
  
  for (field in unique(rotation_data$field)) {
    field_crops <- rotation_data$crop[rotation_data$field == field]
    unique_crops <- length(unique(field_crops))
    total_years <- length(field_crops)
    
    # Simple diversity score (0-1, higher = more diverse)
    diversity_scores[[field]] <- unique_crops / total_years
  }
  
  return(diversity_scores)
}

diversity_scores <- calculate_crop_diversity(rotation_analysis)

print("Crop Diversity Scores (Higher = More Diverse):")
for (field in names(diversity_scores)) {
  print(paste(field, ":", round(diversity_scores[[field]], 2)))
}

# Sustainability recommendations
avg_diversity <- mean(unlist(diversity_scores))
print(paste("Average field diversity:", round(avg_diversity, 2)))

if (avg_diversity < 0.5) {
  print("RECOMMENDATION: Consider more diverse crop rotations")
} else if (avg_diversity > 0.8) {
  print("STATUS: Good crop diversity maintained")
}

# Create farm-ready data export
create_farm_export <- function(analysis_results, field_boundaries) {
  # Compile key metrics for farm management
  farm_data <- list(
    summary_statistics = list(),
    recommendations = list(),
    alerts = list()
  )
  
  # Extract key metrics
  if (exists("veg_stats") && length(veg_stats) > 0) {
    farm_data$summary_statistics <- veg_stats
  }
  
  # Generate management recommendations
  if (exists("alerts")) {
    if (alerts$drought_risk > 20) {
      farm_data$recommendations <- c(farm_data$recommendations, 
                                   "Increase irrigation frequency")
    }
    if (alerts$disease_risk > 15) {
      farm_data$recommendations <- c(farm_data$recommendations,
                                   "Monitor for disease symptoms")
    }
  }
  
  return(farm_data)
}

# Export field-level results
if (exists("field_analysis") && inherits(field_analysis, "sf")) {
  field_summary <- sf::st_drop_geometry(field_analysis)
  
  # Add coordinates for GPS guidance
  field_centroids <- sf::st_centroid(field_analysis)
  field_coords <- sf::st_coordinates(field_centroids)
  field_summary$centroid_lon <- field_coords[, 1]
  field_summary$centroid_lat <- field_coords[, 2]
  
  print("Field Summary for Farm Management:")
  print(field_summary)
}

# Estimate economic value based on vegetation health
calculate_economic_metrics <- function(ndvi_data, crop_prices) {
  if (!exists("ndvi_values")) return(NULL)
  
  # Simplified yield prediction based on NDVI
  # (In practice, would use crop-specific models)
  predicted_yield <- ifelse(ndvi_values > 0.7, "High",
                           ifelse(ndvi_values > 0.5, "Medium", "Low"))
  
  yield_table <- table(predicted_yield)
  
  # Economic projections (simplified)
  economic_zones <- list(
    high_yield = sum(predicted_yield == "High", na.rm = TRUE),
    medium_yield = sum(predicted_yield == "Medium", na.rm = TRUE), 
    low_yield = sum(predicted_yield == "Low", na.rm = TRUE)
  )
  
  return(economic_zones)
}

# Calculate economic zones
if (exists("ndvi_values")) {
  economic_analysis <- calculate_economic_metrics(ndvi_values, 
                                                  list(corn = 5.50, soybeans = 13.00))
  
  if (!is.null(economic_analysis)) {
    total_pixels <- sum(unlist(economic_analysis))
    
    print("Economic Potential Analysis:")
    for (zone in names(economic_analysis)) {
      pixels <- economic_analysis[[zone]]
      percentage <- round((pixels / total_pixels) * 100, 1)
      print(paste(zone, ":", pixels, "pixels (", percentage, "%)"))
    }
  }
}

# Vegetation index quality assessment
quality_check <- function(vegetation_indices) {
  qc_results <- list()
  
  if (inherits(vegetation_indices, "SpatRaster")) {
    for (i in 1:terra::nlyr(vegetation_indices)) {
      index_name <- names(vegetation_indices)[i]
      values <- terra::values(vegetation_indices[[i]], mat = FALSE)
      
      qc_results[[index_name]] <- list(
        total_pixels = length(values),
        valid_pixels = sum(!is.na(values)),
        coverage_percent = round((sum(!is.na(values)) / length(values)) * 100, 1),
        value_range = range(values, na.rm = TRUE),
        outliers = sum(values < -1 | values > 1, na.rm = TRUE)  # For normalized indices
      )
    }
  }
  
  return(qc_results)
}

# Perform quality check
if (exists("agricultural_indices")) {
  qc_results <- quality_check(agricultural_indices)
  
  print("Data Quality Assessment:")
  for (index in names(qc_results)) {
    qc <- qc_results[[index]]
    print(paste(index, "- Coverage:", qc$coverage_percent, "%, Outliers:", qc$outliers))
  }
}

# Create sampling points for field validation
create_validation_points <- function(field_boundary, n_points = 10) {
  # Generate random points within field boundary
  if (inherits(field_boundary, "sf")) {
    sample_points <- sf::st_sample(field_boundary, n_points)
    validation_sf <- sf::st_sf(
      point_id = paste0("VP_", 1:length(sample_points)),
      validation_type = "Ground_Truth",
      geometry = sample_points
    )
    return(validation_sf)
  }
  return(NULL)
}

# Create validation points for first field
if (exists("fields_sf")) {
  validation_points <- create_validation_points(fields_sf[1, ], n_points = 5)
  
  if (!is.null(validation_points)) {
    print("Validation Points Created:")
    print(sf::st_drop_geometry(validation_points))
    
    # Extract vegetation index values at validation points
    if (exists("agricultural_indices")) {
      validation_extracted <- universal_spatial_join(
        source_data = validation_points,
        target_data = agricultural_indices,
        method = "extract",
        verbose = FALSE
      )
      
      print("Extracted Values at Validation Points:")
      validation_summary <- sf::st_drop_geometry(validation_extracted)
      print(head(validation_summary))
    }
  }
}

# Complete agricultural analysis workflow
run_agricultural_workflow <- function(spectral_data, cdl_data = NULL, 
                                     region_boundary = NULL) {
  workflow_results <- list()
  
  # Step 1: Calculate vegetation indices
  message("Step 1: Calculating vegetation indices...")
  indices <- calculate_multiple_indices(
    spectral_data = spectral_data,
    indices = c("NDVI", "SAVI", "DVI"),
    output_stack = TRUE
  )
  workflow_results$vegetation_indices <- indices
  
  # Step 2: Crop classification (if CDL available)
  if (!is.null(cdl_data)) {
    message("Step 2: Analyzing crop distribution...")
    crop_mask <- create_crop_mask(cdl_data, "corn")
    workflow_results$crop_mask <- crop_mask
  }
  
  # Step 3: Stress assessment
  message("Step 3: Assessing crop stress...")
  if (inherits(indices, "SpatRaster") && "NDVI" %in% names(indices)) {
    ndvi_vals <- terra::values(indices[["NDVI"]], mat = FALSE)
    stress_percent <- sum(ndvi_vals < 0.5, na.rm = TRUE) / sum(!is.na(ndvi_vals)) * 100
    workflow_results$stress_assessment <- round(stress_percent, 1)
  }
  
  # Step 4: Generate recommendations
  message("Step 4: Generating management recommendations...")
  recommendations <- character()
  
  if (!is.null(workflow_results$stress_assessment)) {
    if (workflow_results$stress_assessment > 25) {
      recommendations <- c(recommendations, "High stress detected - investigate causes")
    }
    if (workflow_results$stress_assessment < 10) {
      recommendations <- c(recommendations, "Crop conditions appear favorable")
    }
  }
  
  workflow_results$recommendations <- recommendations
  
  return(workflow_results)
}

# Run the workflow
workflow_output <- run_agricultural_workflow(
  spectral_data = spectral_stack
)

print("Complete Agricultural Workflow Results:")
print(paste("Stress assessment:", workflow_output$stress_assessment, "% of area"))
print("Recommendations:")
for (rec in workflow_output$recommendations) {
  print(paste("-", rec))
}

# 1. Always validate your data
print("Data Validation Checklist:")
print("✓ Check coordinate reference systems")
print("✓ Verify date ranges match growing season")
print("✓ Validate vegetation index ranges")
print("✓ Confirm crop mask accuracy")

# 2. Use appropriate indices for your crop type
crop_index_recommendations <- list(
  corn = c("NDVI", "EVI2", "SAVI", "DVI"),
  soybeans = c("NDVI", "SAVI", "GNDVI", "DVI"),
  wheat = c("NDVI", "SAVI", "MSAVI"),
  general = c("NDVI", "SAVI", "DVI", "RVI")
)

print("Recommended indices by crop type:")
for (crop in names(crop_index_recommendations)) {
  indices <- crop_index_recommendations[[crop]]
  print(paste(crop, ":", paste(indices, collapse = ", ")))
}

# 3. Monitor throughout growing season
print("Seasonal Monitoring Schedule:")
print("- Planting: Establish baseline measurements")
print("- Early season: Monitor emergence and establishment") 
print("- Mid-season: Peak growth assessment and stress detection")
print("- Late season: Maturity evaluation and harvest timing")
print("- Post-harvest: Residue analysis and planning")

# Precision agriculture workflow components
precision_ag_components <- list(
  data_collection = c("Satellite imagery", "UAV surveys", "Ground sensors"),
  analysis_methods = c("Vegetation indices", "Stress detection", "Yield mapping"),
  decision_support = c("Variable rate applications", "Irrigation scheduling", "Harvest timing"),
  validation = c("Ground truth sampling", "Yield monitoring", "Economic analysis")
)

print("Precision Agriculture Integration:")
for (component in names(precision_ag_components)) {
  print(paste(component, ":", paste(precision_ag_components[[component]], collapse = ", ")))
}