Using analyze_crop_vegetation() in geospatialsuite

Complete Documentation: analyze_crop_vegetation() Output Structure

Overview

The analyze_crop_vegetation() function returns a nested list with three main components containing comprehensive crop analysis results.

Output Structure

red <- load_sample_data("sample_red.rds")
nir <- load_sample_data("sample_nir.rds")
blue <- load_sample_data("sample_blue.rds")

spectral_stack <- c(red, nir, blue)

names(spectral_stack) <- c("red", "nir", "blue")

result <- analyze_crop_vegetation(
  spectral_data = spectral_stack,
  crop_type = "corn",
  analysis_type = "comprehensive"
)

# Structure:
result$vegetation_indices      # SpatRaster with calculated indices
result$analysis_results        # Detailed analysis results
result$metadata                # Processing metadata

1. vegetation_indices (SpatRaster)

Type: terra::SpatRaster object (multi-layer)

Description: A raster stack containing all calculated vegetation indices for the specified crop type.

Contents:

Each layer is a different vegetation index (NDVI, EVI, GNDVI, etc.)
Values are the actual index calculations for each pixel
Layer names correspond to the index abbreviations

Example Access:

# View all calculated indices
names(result$vegetation_indices)
# [1] "NDVI" "EVI" "GNDVI" "DVI" "RVI" "PRI"

# Extract a specific index
ndvi <- result$vegetation_indices[["NDVI"]]

# Plot an index
plot(result$vegetation_indices[["NDVI"]], main = "NDVI")

# Get values for analysis
ndvi_values <- terra::values(result$vegetation_indices[["NDVI"]])

2. analysis_results (List)

Type: Named list with up to 5 components depending on analysis_type

2.1 stress_analysis (if analysis_type includes “stress” or “comprehensive”)

Purpose: Identifies vegetation stress based on index thresholds

Structure:

result$analysis_results$stress_analysis

Contents for each index analyzed (e.g., NDVI, EVI, SIPI, GNDVI):

Field	Type	Description	Example
`healthy_percentage`	numeric	% of pixels classified as healthy vegetation	65.3
`moderate_stress_percentage`	numeric	% of pixels showing moderate stress	25.1
`severe_stress_percentage`	numeric	% of pixels showing severe stress	9.6
`mean_value`	numeric	Mean index value across all pixels	0.72
`median_value`	numeric	Median index value	0.74
`std_dev`	numeric	Standard deviation of index values	0.15
`thresholds_used`	list	The threshold values used for classification	See below
`total_pixels_analyzed`	integer	Total number of valid pixels	125000

Threshold Structure:

result$analysis_results$stress_analysis$NDVI$thresholds_used
# $healthy: c(0.6, 1.0)         # NDVI 0.6-1.0 = healthy
# $moderate_stress: c(0.4, 0.6) # NDVI 0.4-0.6 = moderate stress
# $severe_stress: c(0.0, 0.4)   # NDVI 0.0-0.4 = severe stress

Interpretation:

Healthy: Vegetation is growing normally, adequate water/nutrients
Moderate Stress: Some stress present, may indicate water deficit or nutrient issues
Severe Stress: Significant stress, requires immediate attention

Example Usage:

# Access stress results for NDVI
ndvi_stress <- result$analysis_results$stress_analysis$NDVI

# What percentage of my field is healthy?
cat(sprintf("Healthy vegetation: %.1f%%\n", ndvi_stress$healthy_percentage))

# What's the average NDVI?
cat(sprintf("Mean NDVI: %.3f\n", ndvi_stress$mean_value))

# Check all indices analyzed
names(result$analysis_results$stress_analysis)

2.2 growth_analysis (if analysis_type includes “growth” or “comprehensive”)

Purpose: Estimates crop growth stage and provides growth statistics

Structure:

result$analysis_results$growth_analysis

Contents:

Index-specific statistics (for NDVI, EVI, GNDVI, DVI):

Field	Type	Description	Example
`mean`	numeric	Mean index value	0.68
`median`	numeric	Median index value	0.70
`std_dev`	numeric	Standard deviation	0.12
`min`	numeric	Minimum value	0.15
`max`	numeric	Maximum value	0.92
`range`	numeric	Max - Min	0.77
`percentiles`	numeric vector	10th, 25th, 75th, 90th percentiles	c(0.45, 0.58, 0.78, 0.85)
`coefficient_of_variation`	numeric	std_dev / mean (measure of variability)	0.18
`n_pixels`	integer	Number of pixels analyzed	125000

Growth stage prediction (overall):

Field	Type	Description	Example
`predicted_growth_stage`	character	Predicted crop growth stage	“reproductive”
`stage_confidence`	numeric	Confidence in prediction (0-1)	0.85
`crop_type_used`	character	Crop type used for classification	“corn”

Growth Stage Classifications:

For Corn:

“emergence”: NDVI < 0.3
“vegetative”: NDVI 0.3-0.6
“reproductive”: NDVI 0.6-0.8
“maturity”: NDVI > 0.8

For Soybeans:

“emergence”: NDVI < 0.4
“vegetative”: NDVI 0.4-0.65
“reproductive”: NDVI 0.65-0.8
“maturity”: NDVI > 0.8

For Wheat:

“tillering”: NDVI < 0.35
“stem_elongation”: NDVI 0.35-0.7
“grain_filling”: NDVI 0.7-0.8
“maturity”: NDVI > 0.8

Example Usage:

# What growth stage is the crop in?
growth <- result$analysis_results$growth_analysis
cat(sprintf("Growth stage: %s (confidence: %.2f)\n", 
            growth$predicted_growth_stage,
            growth$stage_confidence))

# Get detailed NDVI statistics
ndvi_stats <- growth$NDVI
cat(sprintf("NDVI range: %.3f - %.3f\n", ndvi_stats$min, ndvi_stats$max))
cat(sprintf("NDVI variability (CV): %.3f\n", ndvi_stats$coefficient_of_variation))

2.3 yield_analysis (if analysis_type includes “yield” or “comprehensive”)

Purpose: Estimates yield potential using multiple vegetation indices

Structure:

result$analysis_results$yield_analysis

Contents:

Field	Type	Description	Example
`composite_yield_index`	numeric	Normalized yield potential score (0-1)	0.72
`yield_potential_class`	character	Categorical yield potential	“High”
`indices_used`	character vector	Which indices contributed	c(“NDVI”, “EVI”, “GNDVI”)
`n_indices_used`	integer	Number of indices used	3
`index_contributions`	list	Individual index contributions	See below
`crop_type`	character	Crop type used	“corn”
`classification_confidence`	numeric	Confidence in classification (0-1)	0.44

Index Contributions Structure:

result$analysis_results$yield_analysis$index_contributions$NDVI
# $mean_normalized: 0.75    # Normalized contribution (0-1)
# $raw_mean: 0.68           # Raw mean NDVI value
# $raw_std: 0.12            # Raw standard deviation

Composite Yield Index Calculation:

Each index (NDVI, EVI, GNDVI, DVI, RVI) is normalized to 0-1 scale
Normalized values are averaged across all available indices
Result is a single 0-1 score where:
- 0.0 = Very low yield potential
- 0.5 = Medium yield potential
- 1.0 = Maximum yield potential

Yield Potential Classifications:

For Corn:

“Low”: composite_index < 0.3
“Medium”: composite_index 0.3-0.6
“High”: composite_index 0.6-0.8
“Very High”: composite_index > 0.8

For Soybeans:

“Low”: composite_index < 0.35
“Medium”: composite_index 0.35-0.65
“High”: composite_index 0.65-0.85
“Very High”: composite_index > 0.85

For Wheat:

“Low”: composite_index < 0.3
“Medium”: composite_index 0.3-0.6
“High”: composite_index 0.6-0.8
“Very High”: composite_index > 0.8

Interpretation:

Composite Yield Index (0.0-1.0): Higher values indicate better yield potential
Yield Potential Class: Categorical assessment for easier interpretation
Index Contributions: Shows which indices contributed and their individual scores
Classification Confidence: Higher when composite_index is far from class boundaries (e.g., 0.2 or 0.9 are more confident than 0.5)

Example Usage:

# What's the yield potential?
yield <- result$analysis_results$yield_analysis
cat(sprintf("Yield Potential: %s (score: %.2f)\n", 
            yield$yield_potential_class,
            yield$composite_yield_index))

# Which indices contributed?
cat(sprintf("Based on %d indices: %s\n",
            yield$n_indices_used,
            paste(yield$indices_used, collapse = ", ")))

# Get individual index contributions
for (idx in names(yield$index_contributions)) {
  contrib <- yield$index_contributions[[idx]]
  cat(sprintf("%s: %.3f (raw: %.3f +/- %.3f)\n",
              idx, 
              contrib$mean_normalized,
              contrib$raw_mean,
              contrib$raw_std))
}

2.4 summary_statistics

Purpose: Basic statistical summary for all calculated indices

Structure:

result$analysis_results$summary_statistics

Contents for each index:

Field	Type	Description
`mean`	numeric	Mean value across all pixels
`median`	numeric	Median value
`std_dev`	numeric	Standard deviation
`min`	numeric	Minimum value
`max`	numeric	Maximum value
`count`	integer	Number of valid pixels
`na_count`	integer	Number of NA pixels
`range`	numeric	Max - Min
`cv`	numeric	Coefficient of variation (std_dev/mean)
`percentiles`	numeric vector	5th, 25th, 75th, 95th percentiles
`coverage_percent`	numeric	% of total pixels with valid data
`histogram`	list (optional)	Histogram data if >= 100 pixels

Plus Overall Summary:

result$analysis_results$summary_statistics$summary
# $total_indices_calculated: 6
# $indices_with_valid_data: c("NDVI", "EVI", "GNDVI", ...)
# $total_indices_requested: 6
# $success_rate: 100.0

Example Usage:

# Get statistics for all indices
stats <- result$analysis_results$summary_statistics

# NDVI statistics
ndvi_stats <- stats$NDVI
cat(sprintf("NDVI: %.3f +/- %.3f (range: %.3f to %.3f)\n",
            ndvi_stats$mean,
            ndvi_stats$std_dev,
            ndvi_stats$min,
            ndvi_stats$max))

# Check data quality
cat(sprintf("Coverage: %.1f%% (%d pixels)\n",
            ndvi_stats$coverage_percent,
            ndvi_stats$count))

2.5 validation (if reference_data provided)

Purpose: Validates analysis against ground truth data

Note: This component only appears if you provide reference_data parameter

Structure: TBD (depends on reference data format)

3. metadata (List)

Purpose: Documents analysis parameters and processing information

Structure:

result$metadata

Contents:

Field	Type	Description	Example
`crop_type`	character	Crop type analyzed	“corn”
`growth_stage`	character	Growth stage specified	“mid”
`analysis_type`	character	Type of analysis performed	“comprehensive”
`indices_used`	character vector	Indices calculated	c(“NDVI”, “EVI”, …)
`processing_date`	POSIXct	When analysis was performed	2025-11-03 10:30:45
`input_bands`	integer	Number of input spectral bands	8
`spatial_resolution`	numeric vector	Resolution (x, y) in map units	c(10, 10)
`spatial_extent`	numeric vector	Extent (xmin, xmax, ymin, ymax)	c(-95.5, -95.0, 41.5, 42.0)

Example Usage:

# Check what was analyzed
meta <- result$metadata
cat(sprintf("Analyzed %s at %s growth stage\n", 
            meta$crop_type, 
            meta$growth_stage))
cat(sprintf("Used %d indices: %s\n",
            length(meta$indices_used),
            paste(meta$indices_used, collapse = ", ")))
cat(sprintf("Processed on: %s\n", meta$processing_date))

Complete Example Workflow

library(geospatialsuite)
library(terra)

# Run comprehensive crop analysis
result <- analyze_crop_vegetation(
  spectral_data = "path/to/sentinel2_data.tif",
  crop_type = "corn",
  growth_stage = "mid",
  analysis_type = "comprehensive",
  verbose = TRUE
)

# ===== 1. Check what was calculated =====
cat("Indices calculated:\n")
print(names(result$vegetation_indices))

# ===== 2. Assess crop stress =====
stress <- result$analysis_results$stress_analysis$NDVI
cat(sprintf("\nStress Assessment:\n"))
cat(sprintf("  Healthy: %.1f%%\n", stress$healthy_percentage))
cat(sprintf("  Moderate stress: %.1f%%\n", stress$moderate_stress_percentage))
cat(sprintf("  Severe stress: %.1f%%\n", stress$severe_stress_percentage))

# ===== 3. Identify growth stage =====
growth <- result$analysis_results$growth_analysis
cat(sprintf("\nGrowth Stage: %s (%.0f%% confidence)\n",
            growth$predicted_growth_stage,
            growth$stage_confidence * 100))

# ===== 4. Estimate yield potential =====
yield <- result$analysis_results$yield_analysis
cat(sprintf("\nYield Potential: %s\n", yield$yield_potential_class))
cat(sprintf("Composite Yield Index: %.3f\n", yield$composite_yield_index))
cat(sprintf("Based on %d indices: %s\n",
            yield$n_indices_used,
            paste(yield$indices_used, collapse = ", ")))

# ===== 5. Visualize results =====
# Plot stress map
plot(result$vegetation_indices[["NDVI"]], 
     main = "NDVI - Stress Detection",
     col = terrain.colors(100))

# Create stress classification map
ndvi <- result$vegetation_indices[["NDVI"]]
stress_map <- classify(ndvi, 
                      matrix(c(-Inf, 0.4, 1,    # Severe stress
                              0.4, 0.6, 2,       # Moderate stress
                              0.6, Inf, 3),      # Healthy
                            ncol = 3, byrow = TRUE))
plot(stress_map, 
     main = "Crop Stress Classification",
     col = c("red", "yellow", "green"),
     legend = FALSE)
legend("topright", 
       legend = c("Severe Stress", "Moderate Stress", "Healthy"),
       fill = c("red", "yellow", "green"))

# ===== 6. Export results =====
# Save as geotiff
writeRaster(result$vegetation_indices, 
            "crop_indices.tif", 
            overwrite = TRUE)

# Save statistics as CSV
stats_df <- data.frame(
  Index = names(result$analysis_results$summary_statistics)[-length(names(result$analysis_results$summary_statistics))],
  Mean = sapply(result$analysis_results$summary_statistics[1:(length(result$analysis_results$summary_statistics)-1)], function(x) x$mean),
  StdDev = sapply(result$analysis_results$summary_statistics[1:(length(result$analysis_results$summary_statistics)-1)], function(x) x$std_dev),
  Min = sapply(result$analysis_results$summary_statistics[1:(length(result$analysis_results$summary_statistics)-1)], function(x) x$min),
  Max = sapply(result$analysis_results$summary_statistics[1:(length(result$analysis_results$summary_statistics)-1)], function(x) x$max)
)
write.csv(stats_df, "crop_analysis_statistics.csv", row.names = FALSE)

Tips for Using Results

1. Identifying Problem Areas

# Find pixels with severe stress
ndvi <- result$vegetation_indices[["NDVI"]]
severe_stress <- ndvi < 0.4
plot(severe_stress, main = "Severe Stress Areas")

2. Comparing Multiple Fields

# Run analysis for multiple fields and compare
field1_result <- analyze_crop_vegetation(field1_data, crop_type = "corn")
field2_result <- analyze_crop_vegetation(field2_data, crop_type = "corn")

# Compare yield potential
cat(sprintf("Field 1 yield: %s (%.3f)\n", 
            field1_result$analysis_results$yield_analysis$yield_potential_class,
            field1_result$analysis_results$yield_analysis$composite_yield_index))
cat(sprintf("Field 2 yield: %s (%.3f)\n",
            field2_result$analysis_results$yield_analysis$yield_potential_class,
            field2_result$analysis_results$yield_analysis$composite_yield_index))

3. Time Series Analysis

# Analyze the same field at different dates
early_season <- analyze_crop_vegetation(june_data, growth_stage = "early")
mid_season <- analyze_crop_vegetation(july_data, growth_stage = "mid")
late_season <- analyze_crop_vegetation(august_data, growth_stage = "late")

# Track NDVI progression
ndvi_progression <- c(
  early_season$analysis_results$growth_analysis$NDVI$mean,
  mid_season$analysis_results$growth_analysis$NDVI$mean,
  late_season$analysis_results$growth_analysis$NDVI$mean
)
plot(1:3, ndvi_progression, type = "b",
     xlab = "Time Period", ylab = "Mean NDVI",
     main = "NDVI Progression Through Season")

Classification Methodology and References

analyze_crop_vegetation() classifications are based on NDVI and other indices. Each analysis component uses a different combination of indices, with thresholds grounded in peer-reviewed remote sensing literature.

Stress Detection

Thresholds are applied independently to each available index, with each index using biophysically appropriate breakpoints (Tucker, 1979; Hatfield et al., 2008):

Index	Healthy	Moderate Stress	Severe Stress
NDVI	>= 0.60	0.40 – 0.60	< 0.40
EVI	>= 0.40	0.20 – 0.40	< 0.20
GNDVI	>= 0.50	0.30 – 0.50	< 0.30
SIPI	>= 1.00	0.80 – 1.00	< 0.80

Growth Stage Classification

Growth stage is predicted from mean NDVI using crop-specific thresholds. Statistics (mean, range, percentiles, CV) are additionally reported for EVI, GNDVI, and DVI. The threshold approach is conceptually consistent with Anyamba et al. (2021), who applied explicit per-crop NDVI cutoffs for large-scale MODIS-based crop monitoring, and with the MODIS + CDL crop monitoring framework of Akanbi et al. (2024).

Corn

Stage	Mean NDVI
Emergence	< 0.30
Vegetative	0.30 – 0.60
Reproductive	0.60 – 0.80
Maturity	> 0.80

Soybean

Stage	Mean NDVI
Emergence	< 0.40
Vegetative	0.40 – 0.65
Reproductive	0.65 – 0.80
Maturity	> 0.80

Wheat

Stage	Mean NDVI
Tillering	< 0.35
Stem elongation	0.35 – 0.70
Grain filling	0.70 – 0.80
Maturity	> 0.80

These are literature-informed starting points. Exact NDVI boundaries vary by cultivar, sensor, region, and season. Calibration with local ground truth data is recommended.

Yield Potential (Composite Index)

Yield potential is assessed using a composite score that normalizes up to five indices (NDVI, EVI, GNDVI, DVI, RVI) to a 0–1 scale and averages them, reducing dependence on any single index. This multi-index approach is consistent with Bolton & Friedl (2013) and Mkhabela et al. (2011).

Composite Score	Corn / Wheat	Soybean	Class
< 0.30	Low	< 0.35	Low
0.30 – 0.60	Medium	0.35 – 0.65	Medium
0.60 – 0.80	High	0.65 – 0.85	High
> 0.80	Very High	> 0.85	Very High

The composite yield index is a vegetation-based proxy, not a direct yield prediction. Correlation with actual harvested yield varies by crop, region, year, and management practice.

References

Tucker, C.J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment, 8(2), 127–150. https://doi.org/10.1016/0034-4257(79)90013-0
Hatfield, J.L., Gitelson, A.A., Schepers, J.S. & Walthall, C.L. (2008). Application of spectral remote sensing for agronomic decisions. Agronomy Journal, 100(S3), S-117–S-131. https://doi.org/10.2134/agronj2006.0370c
Anyamba, A. et al. (2021). USA crop yield estimation with MODIS NDVI: Impact of the 2012 drought on corn production. Remote Sensing, 13(21), 4227. https://doi.org/10.3390/rs13214227
Akanbi, O.D. et al. (2024). Integrating multiscale geospatial analysis for monitoring crop growth, nutrient distribution, and hydrological dynamics in large-scale agricultural systems. Journal of Geovisualization and Spatial Analysis, 8, 9. https://doi.org/10.1007/s41651-023-00164-y
Bolton, D.K. & Friedl, M.A. (2013). Forecasting crop yield using remotely sensed vegetation indices and crop phenology metrics. Agricultural and Forest Meteorology, 173, 74–84. https://doi.org/10.1016/j.agrformet.2013.01.007
Mkhabela, M.S. et al. (2011). Crop yield forecasting on the Canadian Prairies using MODIS NDVI data. Agricultural and Forest Meteorology, 151(3), 385–393. https://doi.org/10.1016/j.agrformet.2010.11.012

Validation Notes

Important Caveats:

Threshold-based: Stress and yield classifications use literature-based thresholds that may need adjustment for your specific region/conditions
Composite Yield Index: This is a vegetation-based proxy, not a direct yield prediction. Correlation with actual yield varies by crop, region, and year
Growth Stage: Predictions are based on NDVI patterns and may not align perfectly with field observations
These are analytical tools to support decision-making, not definitive assessments

Recommended Validation:

Compare with ground truth data (yield monitors, field scouting)
Calibrate thresholds for your specific conditions
Use multiple years of data to establish local patterns
Combine with other data sources (weather, soil, management)

Acknowledgments

This work was developed by the geospatialsuite team with contributions from: Olatunde D. Akanbi, Vibha Mandayam, Yinghui Wu, Jeffrey Yarus, Erika I. Barcelos, and Roger H. French.

geospatialsuite Development Team

Complete Documentation: analyze_crop_vegetation() Output Structure

Overview

Output Structure

1. vegetation_indices (SpatRaster)

2. analysis_results (List)

2.1 stress_analysis (if analysis_type includes “stress” or “comprehensive”)

2.2 growth_analysis (if analysis_type includes “growth” or “comprehensive”)

Index-specific statistics (for NDVI, EVI, GNDVI, DVI):

Growth stage prediction (overall):

2.3 yield_analysis (if analysis_type includes “yield” or “comprehensive”)

2.4 summary_statistics

2.5 validation (if reference_data provided)

3. metadata (List)

Complete Example Workflow

Tips for Using Results

1. Identifying Problem Areas

2. Comparing Multiple Fields

3. Time Series Analysis

Classification Methodology and References

Stress Detection

Growth Stage Classification

Yield Potential (Composite Index)

References

Validation Notes

Acknowledgments