inst/doc/agricultural-analysis.R

In geospatialsuite: Comprehensive Geospatiotemporal Analysis and Multimodal Integration Toolkit

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE,
  fig.width = 10,
  fig.height = 8
)

## ----enhanced-ndvi, eval=FALSE------------------------------------------------
# # 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 = " - ")))

## ----multi-index-analysis, eval=FALSE-----------------------------------------
# # 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"
#     )
#   }
# }

## ----yield-assessment, eval=FALSE---------------------------------------------
# # 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"
# )

## ----crop-performance, eval=FALSE---------------------------------------------
# # 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 = "-")))
#   }
# }

## ----field-analysis, eval=FALSE-----------------------------------------------
# # 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))
# }

## ----variable-rate, eval=FALSE------------------------------------------------
# # 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"
#   )
# }

## ----stress-monitoring, eval=FALSE--------------------------------------------
# # 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]], "%)"))
#   }
# }

## ----early-warning, eval=FALSE------------------------------------------------
# # 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")
#   }
# }

## ----seasonal-monitoring, eval=FALSE------------------------------------------
# # 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)
# }

## ----harvest-timing, eval=FALSE-----------------------------------------------
# # 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, "%)"))
#   }
# }

## ----irrigation-assessment, eval=FALSE----------------------------------------
# # 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"
#   )
# }

## ----crop-rotation, eval=FALSE------------------------------------------------
# # 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]]))
# }

## ----sustainability-metrics, eval=FALSE---------------------------------------
# # 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")
# }

## ----farm-integration, eval=FALSE---------------------------------------------
# # 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)
# }

## ----economic-analysis, eval=FALSE--------------------------------------------
# # 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, "%)"))
#     }
#   }
# }

## ----quality-assurance, eval=FALSE--------------------------------------------
# # 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))
#   }
# }

## ----field-validation, eval=FALSE---------------------------------------------
# # 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-workflow, eval=FALSE--------------------------------------------
# # 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))
# }

## ----best-practices, eval=FALSE-----------------------------------------------
# # 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-ag-integration, eval=FALSE-------------------------------------
# # 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 = ", ")))
# }

## ----stress-example, eval=FALSE-----------------------------------------------
# result <- analyze_crop_vegetation(data, analysis_type = "stress")
# stress <- result$analysis_results$stress_analysis$NDVI
# 
# # What percentage of my field needs attention?
# cat(sprintf("%.1f%% of field shows stress\n",
#             stress$moderate_stress_percentage + stress$severe_stress_percentage))

## ----yield-example, eval=FALSE------------------------------------------------
# result <- analyze_crop_vegetation(data, crop_type = "corn", analysis_type = "yield")
# yield <- result$analysis_results$yield_analysis
# 
# cat(sprintf("Yield Potential: %s\n", yield$yield_potential_class))
# cat(sprintf("Composite Score: %.2f\n", yield$composite_yield_index))
# 
# # See which indices contributed
# for (idx in names(yield$index_contributions)) {
#   contrib <- yield$index_contributions[[idx]]
#   cat(sprintf("  %s: %.3f\n", idx, contrib$mean_normalized))
# }

## ----growth-example, eval=FALSE-----------------------------------------------
# result <- analyze_crop_vegetation(data, crop_type = "soybeans", analysis_type = "growth")
# growth <- result$analysis_results$growth_analysis
# 
# cat(sprintf("Predicted stage: %s\n", growth$predicted_growth_stage))
# cat(sprintf("Confidence: %.0f%%\n", growth$stage_confidence * 100))