---
title: "Vegetation Index Analysis with GeoSpatialSuite"
author: "GeoSpatialSuite Development Team"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Vegetation Index Analysis with GeoSpatialSuite}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  warning = FALSE,
  message = FALSE
)
```

# Introduction

This vignette demonstrates comprehensive vegetation analysis using GeoSpatialSuite's 60+ vegetation indices. Learn to monitor plant health, detect stress, and analyze agricultural productivity using satellite imagery.

## Overview of Available Indices

GeoSpatialSuite supports over 60 vegetation indices across multiple categories:

```{r eval=FALSE}
library(geospatialsuite)

# See all available vegetation indices
all_indices <- list_vegetation_indices()
print(all_indices)

# View detailed information with formulas
detailed_indices <- list_vegetation_indices(detailed = TRUE)
head(detailed_indices)

# Filter by category
basic_indices <- list_vegetation_indices(category = "basic")
stress_indices <- list_vegetation_indices(category = "stress")
```

### Index Categories

**Basic Vegetation Indices (10):**
- NDVI, SAVI, MSAVI, OSAVI, EVI, EVI2, DVI, RVI, GNDVI, WDVI

**Enhanced/Improved Indices (12):**
- ARVI, RDVI, PVI, IPVI, TNDVI, GEMI, VARI, TSAVI, ATSAVI, GESAVI, MTVI, CTVI

**Red Edge and Advanced Indices (10):**
- NDRE, MTCI, IRECI, S2REP, PSRI, CRI1, CRI2, ARI1, ARI2, MCARI

**Stress and Chlorophyll Indices (12):**
- PRI, SIPI, CCI, NDNI, CARI, TCARI, MTVI1, MTVI2, TVI, NPCI, RARS, NPQI

## Single Date Analysis

### Basic NDVI Calculation

The most common vegetation analysis:

```{r eval=FALSE}
# Simple NDVI from individual bands
ndvi <- calculate_vegetation_index(
  red = "landsat_red.tif",
  nir = "landsat_nir.tif",
  index_type = "NDVI",
  verbose = TRUE
)

# View results
terra::plot(ndvi, main = "NDVI Analysis", 
            col = colorRampPalette(c("brown", "yellow", "green"))(100))

# Check value range
summary(terra::values(ndvi, mat = FALSE))
```

### Multiple Index Calculation

Calculate several indices at once for comprehensive analysis:

```{r eval=FALSE}
# Calculate multiple vegetation indices
vegetation_stack <- calculate_multiple_indices(
  red = red_band,
  nir = nir_band, 
  blue = blue_band,
  green = green_band,
  indices = c("NDVI", "EVI", "SAVI", "GNDVI", "DVI", "RVI"),
  output_stack = TRUE,
  region_boundary = "Ohio",
  verbose = TRUE
)

# Access individual indices
ndvi_layer <- vegetation_stack$NDVI
evi_layer <- vegetation_stack$EVI

# Create comparison visualization
terra::plot(vegetation_stack, main = names(vegetation_stack))
```

### Multi-Band Data Processing

Work with multi-band satellite imagery:

```{r eval=FALSE}
# From multi-band raster with auto-detection
evi <- calculate_vegetation_index(
  spectral_data = "sentinel2_image.tif",
  index_type = "EVI",
  auto_detect_bands = TRUE,
  verbose = TRUE
)

# From directory with band detection
savi <- calculate_vegetation_index(
  spectral_data = "/path/to/sentinel/bands/",
  index_type = "SAVI", 
  auto_detect_bands = TRUE
)

# Custom band names for non-standard data
ndvi_custom <- calculate_vegetation_index(
  spectral_data = custom_satellite_data,
  band_names = c("B4", "B3", "B2", "B8"),  # Red, Green, Blue, NIR
  index_type = "NDVI"
)
```

## Time Series Analysis

### Enhanced NDVI for Temporal Analysis

Use `calculate_ndvi_enhanced()` for time series:

```{r eval=FALSE}
# Time series NDVI with quality control
ndvi_series <- calculate_ndvi_enhanced(
  red_data = "/path/to/red/time_series/",
  nir_data = "/path/to/nir/time_series/",
  match_by_date = TRUE,           # Automatically match files by date
  quality_filter = TRUE,          # Remove outliers
  temporal_smoothing = TRUE,      # Smooth time series
  clamp_range = c(-0.2, 1),
  verbose = TRUE
)

# Plot time series layers
terra::plot(ndvi_series, main = paste("NDVI", names(ndvi_series)))

# Calculate temporal statistics
temporal_mean <- terra::app(ndvi_series, mean, na.rm = TRUE)
temporal_cv <- terra::app(ndvi_series, function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE))
```

### Temporal Change Analysis

Analyze vegetation changes over time:

```{r eval=FALSE}
# Temporal analysis workflow
temporal_results <- analyze_temporal_changes(
  data_list = c("ndvi_2020.tif", "ndvi_2021.tif", "ndvi_2022.tif", "ndvi_2023.tif"),
  dates = c("2020", "2021", "2022", "2023"),
  region_boundary = "Iowa",
  analysis_type = "trend",
  output_folder = "temporal_analysis/"
)

# Access trend results
slope_raster <- temporal_results$trend_rasters$slope
r_squared_raster <- temporal_results$trend_rasters$r_squared

# Interpret trends
# Positive slope = increasing vegetation
# Negative slope = declining vegetation  
# High R² = strong temporal trend
```

## Agricultural Applications

### Crop-Specific Analysis

Analyze vegetation for specific crops:

```{r eval=FALSE}
# Comprehensive crop vegetation analysis
corn_analysis <- analyze_crop_vegetation(
  spectral_data = "sentinel2_data.tif",
  crop_type = "corn",
  growth_stage = "mid",
  analysis_type = "comprehensive",
  cdl_mask = corn_mask,
  verbose = TRUE
)

# Access results
vegetation_indices <- corn_analysis$vegetation_indices
stress_analysis <- corn_analysis$analysis_results$stress_analysis
yield_analysis <- corn_analysis$analysis_results$yield_analysis

# View stress detection results
print(stress_analysis$NDVI)
```

### Growth Stage Monitoring

Monitor crop development through the season:

```{r eval=FALSE}
# Early season analysis (emergence to vegetative)
early_season <- analyze_crop_vegetation(
  spectral_data = "may_sentinel2.tif", 
  crop_type = "soybeans",
  growth_stage = "early",
  analysis_type = "growth"
)

# Mid-season analysis (reproductive stage)  
mid_season <- analyze_crop_vegetation(
  spectral_data = "july_sentinel2.tif",
  crop_type = "soybeans", 
  growth_stage = "mid",
  analysis_type = "stress"
)

# Compare growth stages
early_ndvi <- early_season$vegetation_indices$NDVI
mid_ndvi <- mid_season$vegetation_indices$NDVI
growth_change <- mid_ndvi - early_ndvi
```

## Stress Detection

### Identifying Plant Stress

Use specialized indices for stress detection:

```{r eval=FALSE}
# Calculate stress-sensitive indices
stress_indices <- calculate_multiple_indices(
  red = red_band,
  nir = nir_band,
  green = green_band,
  red_edge = red_edge_band,     # If available
  indices = c("NDVI", "PRI", "SIPI", "NDRE"),
  output_stack = TRUE
)

# Stress analysis thresholds
ndvi_values <- terra::values(stress_indices$NDVI, mat = FALSE)
healthy_threshold <- 0.6
stress_threshold <- 0.4

# Classify stress levels
stress_mask <- terra::classify(stress_indices$NDVI, 
                               matrix(c(-1, stress_threshold, 1,      # Severe stress
                                       stress_threshold, healthy_threshold, 2,  # Moderate stress  
                                       healthy_threshold, 1, 3), ncol = 3, byrow = TRUE))

# Visualization
stress_colors <- c("red", "orange", "green")
terra::plot(stress_mask, main = "Vegetation Stress Classification",
            col = stress_colors, legend = c("Severe", "Moderate", "Healthy"))
```

### Water Stress Detection

Monitor vegetation water content:

```{r eval=FALSE}
# Water content indices
water_stress <- calculate_multiple_indices(
  nir = nir_band,
  swir1 = swir1_band,
  indices = c("NDMI", "MSI", "NDII"),
  output_stack = TRUE
)

# NDMI interpretation:
# High NDMI (>0.4) = High water content
# Low NDMI (<0.1) = Water stress

# MSI interpretation (opposite of NDMI):
# Low MSI (<1.0) = High moisture
# High MSI (>1.6) = Water stress

# Create water stress map
water_stress_binary <- terra::classify(water_stress$NDMI,
                                       matrix(c(-1, 0.1, 1,    # Water stress
                                               0.1, 0.4, 2,    # Moderate 
                                               0.4, 1, 3), ncol = 3, byrow = TRUE))
```

## Quality Control and Validation

### Quality Filtering

Apply quality controls to vegetation data:

```{r eval=FALSE}
# Enhanced NDVI with quality filtering
ndvi_quality <- calculate_ndvi_enhanced(
  red_data = red_raster,
  nir_data = nir_raster,
  quality_filter = TRUE,          # Remove outliers
  clamp_range = c(-0.2, 1),      # Reasonable NDVI range
  temporal_smoothing = FALSE,     # For single date
  verbose = TRUE
)

# Custom quality filtering
vegetation_filtered <- calculate_vegetation_index(
  red = red_band,
  nir = nir_band,
  index_type = "NDVI",
  mask_invalid = TRUE,            # Remove invalid values
  clamp_range = c(-1, 1)         # Theoretical NDVI range
)
```

### Validation with Reference Data

```{r eval=FALSE}
# Compare with field measurements
field_ndvi_validation <- universal_spatial_join(
  source_data = "field_measurements.csv",  # Ground truth points
  target_data = ndvi_raster,              # Satellite NDVI
  method = "extract",
  buffer_distance = 30,                   # 30m buffer around points
  summary_function = "mean"
)

# Calculate correlation
observed <- field_ndvi_validation$field_ndvi
predicted <- field_ndvi_validation$extracted_NDVI

correlation <- cor(observed, predicted, use = "complete.obs")
rmse <- sqrt(mean((observed - predicted)^2, na.rm = TRUE))

print(paste("Validation R =", round(correlation, 3)))
print(paste("RMSE =", round(rmse, 4)))
```

## Advanced Applications

### Phenology Analysis

Track vegetation phenology over time:

```{r eval=FALSE}
# Create NDVI time series for phenology
phenology_analysis <- analyze_temporal_changes(
  data_list = list.files("/ndvi/2023/", pattern = "*.tif", full.names = TRUE),
  dates = extract_dates_universal(list.files("/ndvi/2023/", pattern = "*.tif")),
  region_boundary = "Iowa", 
  analysis_type = "seasonal",
  output_folder = "phenology_results/"
)

# Access seasonal patterns
spring_ndvi <- phenology_analysis$seasonal_rasters$Spring
summer_ndvi <- phenology_analysis$seasonal_rasters$Summer
fall_ndvi <- phenology_analysis$seasonal_rasters$Fall

# Calculate growing season metrics
peak_season <- terra::app(c(spring_ndvi, summer_ndvi), max, na.rm = TRUE)
growing_season_range <- summer_ndvi - spring_ndvi
```

### Multi-Sensor Integration

Combine data from different sensors:

```{r eval=FALSE}
# Landsat NDVI (30m resolution)
landsat_ndvi <- calculate_vegetation_index(
  red = "landsat8_red.tif",
  nir = "landsat8_nir.tif", 
  index_type = "NDVI"
)

# MODIS NDVI (250m resolution)  
modis_ndvi <- calculate_vegetation_index(
  red = "modis_red.tif",
  nir = "modis_nir.tif",
  index_type = "NDVI"
)

# Resample MODIS to Landsat resolution for comparison
modis_resampled <- universal_spatial_join(
  source_data = modis_ndvi,
  target_data = landsat_ndvi,
  method = "resample",
  summary_function = "bilinear"
)

# Compare sensors
sensor_difference <- landsat_ndvi - modis_resampled
correlation <- calculate_spatial_correlation(landsat_ndvi, modis_resampled)
```

## Specialized Indices by Application

### Forest Monitoring

```{r eval=FALSE}
# Forest-specific indices
forest_indices <- calculate_multiple_indices(
  red = red_band,
  nir = nir_band,
  swir1 = swir1_band,
  swir2 = swir2_band,
  indices = c("NDVI", "EVI", "NBR", "NDMI"),  # Normalized Burn Ratio, moisture
  region_boundary = "forest_boundary.shp",
  output_stack = TRUE
)

# Burn severity assessment using NBR
pre_fire_nbr <- forest_indices$NBR  # Before fire
post_fire_nbr <- calculate_vegetation_index(
  red = "post_fire_red.tif",
  nir = "post_fire_nir.tif", 
  swir2 = "post_fire_swir2.tif",
  index_type = "NBR"
)

# Calculate burn severity (dNBR)
burn_severity <- pre_fire_nbr - post_fire_nbr

# Classify burn severity
severity_classes <- terra::classify(burn_severity,
  matrix(c(-Inf, 0.1, 1,      # Unburned
           0.1, 0.27, 2,       # Low severity
           0.27, 0.44, 3,      # Moderate-low
           0.44, 0.66, 4,      # Moderate-high  
           0.66, Inf, 5), ncol = 3, byrow = TRUE))
```

### Grassland and Rangeland

```{r eval=FALSE}
# Grassland-specific analysis
grassland_health <- calculate_multiple_indices(
  red = red_band,
  nir = nir_band,
  green = green_band,
  indices = c("NDVI", "SAVI", "MSAVI", "GNDVI"),  # Soil-adjusted indices
  output_stack = TRUE
)

# Analyze grazing impact
grazed_area_ndvi <- terra::mask(grassland_health$NDVI, grazing_areas)
ungrazed_area_ndvi <- terra::mask(grassland_health$NDVI, ungrazed_areas)

# Compare means
grazed_mean <- terra::global(grazed_area_ndvi, "mean", na.rm = TRUE)
ungrazed_mean <- terra::global(ungrazed_area_ndvi, "mean", na.rm = TRUE)

grazing_impact <- ungrazed_mean - grazed_mean
```

### Urban Vegetation

```{r eval=FALSE}
# Urban vegetation analysis with atmospheric correction
urban_vegetation <- calculate_vegetation_index(
  red = urban_red,
  nir = urban_nir,
  blue = urban_blue,
  index_type = "ARVI",           # Atmospherically Resistant VI
  region_boundary = "city_boundary.shp"
)

# Green space analysis
green_space_threshold <- 0.3
urban_greenness <- terra::classify(urban_vegetation,
  matrix(c(-1, green_space_threshold, 0,           # Non-vegetated
           green_space_threshold, 1, 1), ncol = 3, byrow = TRUE))  # Vegetated

# Calculate green space statistics
total_urban_area <- terra::ncell(urban_greenness) * prod(terra::res(urban_greenness))
green_area <- terra::global(urban_greenness, "sum", na.rm = TRUE) * prod(terra::res(urban_greenness))
green_space_percentage <- (green_area / total_urban_area) * 100
```

## Index Selection Guide

### When to Use Each Index

**NDVI (Normalized Difference Vegetation Index):**
- **Best for:** General vegetation monitoring, biomass estimation
- **Range:** -1 to 1 (vegetation typically 0.3-0.9)
- **Limitations:** Saturates in dense vegetation

```{r eval=FALSE}
ndvi <- calculate_vegetation_index(red = red, nir = nir, index_type = "NDVI")
```

**EVI (Enhanced Vegetation Index):**
- **Best for:** Dense vegetation, reducing atmospheric effects
- **Range:** -1 to 3 (vegetation typically 0.2-1.0)
- **Advantages:** Less saturation than NDVI

```{r eval=FALSE}
evi <- calculate_vegetation_index(red = red, nir = nir, blue = blue, index_type = "EVI")
```

**SAVI (Soil Adjusted Vegetation Index):**
- **Best for:** Areas with exposed soil, early season crops
- **Range:** -1 to 1.5
- **Advantages:** Reduces soil background effects

```{r eval=FALSE}
savi <- calculate_vegetation_index(red = red, nir = nir, index_type = "SAVI")
```

**PRI (Photochemical Reflectance Index):**
- **Best for:** Plant stress detection, photosynthetic efficiency
- **Range:** -1 to 1
- **Applications:** Early stress detection before visible symptoms

```{r eval=FALSE}
# Note: PRI typically uses 531nm and 570nm bands
# Using green and NIR as approximation
pri <- calculate_vegetation_index(green = green, nir = nir, index_type = "PRI")
```

### Index Combinations for Different Applications

**Crop Health Monitoring:**
```{r eval=FALSE}
crop_health <- calculate_multiple_indices(
  red = red, nir = nir, green = green,
  indices = c("NDVI", "GNDVI", "SIPI"),  # Structure Insensitive Pigment Index
  output_stack = TRUE
)
```

**Drought Monitoring:**
```{r eval=FALSE}
drought_indices <- calculate_multiple_indices(
  nir = nir, swir1 = swir1,
  indices = c("NDMI", "MSI"),  # Moisture indices
  output_stack = TRUE
)
```

**Precision Agriculture:**
```{r eval=FALSE}
precision_ag <- calculate_multiple_indices(
  red = red, nir = nir, green = green, red_edge = red_edge,
  indices = c("NDVI", "NDRE", "GNDVI", "CCI"),  # Chlorophyll indices
  output_stack = TRUE
)
```

## Visualization and Interpretation

### Custom Visualization

```{r eval=FALSE}
# NDVI-specific visualization
create_spatial_map(
  spatial_data = ndvi_raster,
  color_scheme = "ndvi",          # NDVI-specific colors
  region_boundary = "study_area.shp",
  title = "NDVI Analysis - Growing Season Peak",
  output_file = "ndvi_analysis.png"
)

# Multi-index comparison
create_comparison_map(
  data1 = ndvi_raster,
  data2 = evi_raster, 
  comparison_type = "side_by_side",
  titles = c("NDVI", "EVI"),
  color_scheme = "viridis"
)
```

### Interactive Exploration

```{r eval=FALSE}
# Interactive vegetation analysis
interactive_veg_map <- create_interactive_map(
  spatial_data = vegetation_points,
  fill_variable = "NDVI",
  popup_vars = c("NDVI", "EVI", "collection_date"),
  basemap = "satellite",
  title = "Interactive Vegetation Analysis"
)

# Save interactive map
htmlwidgets::saveWidget(interactive_veg_map, "vegetation_interactive.html")
```

## Statistical Analysis

### Vegetation Statistics

```{r eval=FALSE}
# Calculate comprehensive statistics
vegetation_stats <- terra::global(vegetation_stack, 
                                  c("mean", "sd", "min", "max"), 
                                  na.rm = TRUE)

print(vegetation_stats)

# Zonal statistics by land cover
landcover_stats <- universal_spatial_join(
  source_data = vegetation_stack,
  target_data = "landcover_polygons.shp",
  method = "zonal", 
  summary_function = "mean"
)

# Statistics by vegetation class
healthy_vegetation <- vegetation_stack$NDVI > 0.6
moderate_vegetation <- vegetation_stack$NDVI > 0.3 & vegetation_stack$NDVI <= 0.6
poor_vegetation <- vegetation_stack$NDVI <= 0.3

# Calculate percentages
total_pixels <- terra::ncell(vegetation_stack$NDVI)
healthy_pct <- terra::global(healthy_vegetation, "sum", na.rm = TRUE) / total_pixels * 100
moderate_pct <- terra::global(moderate_vegetation, "sum", na.rm = TRUE) / total_pixels * 100
poor_pct <- terra::global(poor_vegetation, "sum", na.rm = TRUE) / total_pixels * 100
```

### Correlation Analysis

```{r eval=FALSE}
# Analyze relationships between indices
index_correlations <- analyze_variable_correlations(
  variable_list = list(
    NDVI = vegetation_stack$NDVI,
    EVI = vegetation_stack$EVI,
    SAVI = vegetation_stack$SAVI,
    GNDVI = vegetation_stack$GNDVI
  ),
  method = "pearson",
  create_plots = TRUE,
  output_folder = "correlation_analysis/"
)

# View correlation matrix
print(index_correlations$correlation_matrix)
```

## Practical Examples

### Example 1: Corn Field Monitoring

Complete workflow for monitoring corn fields:

```{r eval=FALSE}
# 1. Define study area
study_area <- get_region_boundary("Iowa:Story")  # Story County, Iowa

# 2. Create corn mask
corn_mask <- create_crop_mask(
  cdl_data = "cdl_iowa_2023.tif",
  crop_codes = get_comprehensive_cdl_codes("corn"),
  region_boundary = study_area
)

# 3. Calculate vegetation indices for corn areas
corn_vegetation <- calculate_multiple_indices(
  spectral_data = "sentinel2_iowa_july.tif",
  indices = c("NDVI", "EVI", "GNDVI", "SAVI"),
  auto_detect_bands = TRUE,
  output_stack = TRUE
)

# 4. Apply corn mask
corn_vegetation_masked <- terra::mask(corn_vegetation, corn_mask)

# 5. Analyze corn health
corn_stats <- terra::global(corn_vegetation_masked, 
                           c("mean", "sd", "min", "max"), 
                           na.rm = TRUE)

# 6. Create visualization
quick_map(corn_vegetation_masked$NDVI, 
          title = "Story County Corn NDVI - July 2023")

# 7. Save results
terra::writeRaster(corn_vegetation_masked, "story_county_corn_indices.tif")
```

### Example 2: Stress Detection Pipeline

Detect and map vegetation stress:

```{r eval=FALSE}
# 1. Load multi-temporal data
april_data <- calculate_ndvi_enhanced("april_sentinel2.tif")
july_data <- calculate_ndvi_enhanced("july_sentinel2.tif")

# 2. Calculate stress indices
stress_indices <- calculate_multiple_indices(
  spectral_data = "july_sentinel2.tif",
  indices = c("NDVI", "PRI", "SIPI", "NDMI"),
  auto_detect_bands = TRUE,
  output_stack = TRUE
)

# 3. Identify stressed areas
# NDVI decline from April to July
ndvi_change <- july_data - april_data
stress_areas <- ndvi_change < -0.2  # Significant decline

# Water stress (low NDMI)
water_stress <- stress_indices$NDMI < 0.2

# Combined stress map
combined_stress <- stress_areas | water_stress

# 4. Quantify stress
total_area <- terra::ncell(combined_stress) * prod(terra::res(combined_stress)) / 10000  # hectares
stressed_pixels <- terra::global(combined_stress, "sum", na.rm = TRUE)
stressed_area <- stressed_pixels * prod(terra::res(combined_stress)) / 10000  # hectares
stress_percentage <- (stressed_area / total_area) * 100

print(paste("Stressed area:", round(stressed_area, 1), "hectares"))
print(paste("Stress percentage:", round(stress_percentage, 1), "%"))
```

### Example 3: Regional Vegetation Assessment

Large-scale vegetation analysis:

```{r eval=FALSE}
# 1. Multi-state analysis
great_lakes_states <- c("Michigan", "Ohio", "Indiana", "Illinois", "Wisconsin")

regional_ndvi <- list()
for (state in great_lakes_states) {
  state_ndvi <- calculate_vegetation_index(
    spectral_data = paste0("/data/", tolower(state), "_modis.tif"),
    index_type = "NDVI", 
    region_boundary = state,
    auto_detect_bands = TRUE
  )
  regional_ndvi[[state]] <- state_ndvi
}

# 2. Create regional mosaic
great_lakes_mosaic <- create_raster_mosaic(
  input_data = regional_ndvi,
  method = "merge",
  region_boundary = "Great Lakes Region"
)

# 3. Regional statistics
regional_stats <- terra::global(great_lakes_mosaic, 
                               c("mean", "sd", "min", "max"), 
                               na.rm = TRUE)

# 4. State-by-state comparison
state_means <- sapply(regional_ndvi, function(x) {
  terra::global(x, "mean", na.rm = TRUE)[[1]]
})

names(state_means) <- great_lakes_states
print(sort(state_means, decreasing = TRUE))
```

## Troubleshooting Common Issues

### Band Detection Problems

```{r eval=FALSE}
# If auto-detection fails, specify band names manually
custom_ndvi <- calculate_vegetation_index(
  spectral_data = "unusual_satellite_data.tif",
  band_names = c("Band_4_Red", "Band_5_NIR"),  # Custom names
  index_type = "NDVI",
  auto_detect_bands = FALSE
)

# Or use individual band approach
manual_ndvi <- calculate_vegetation_index(
  red = satellite_data$Band_4_Red,
  nir = satellite_data$Band_5_NIR,
  index_type = "NDVI"
)
```

### CRS and Projection Issues

```{r eval=FALSE}
# Automatic CRS fixing
robust_indices <- calculate_multiple_indices(
  red = "red_utm.tif",        # UTM projection
  nir = "nir_geographic.tif", # Geographic coordinates
  indices = c("NDVI", "EVI"),
  auto_crs_fix = TRUE,        # Automatically handle CRS differences
  verbose = TRUE              # See what's being fixed
)
```

### Memory Management for Large Areas

```{r eval=FALSE}
# Process large areas efficiently
# 1. Use chunked processing
large_area_ndvi <- calculate_vegetation_index(
  spectral_data = "very_large_raster.tif",
  index_type = "NDVI",
  chunk_size = 500000,    # Smaller chunks
  auto_detect_bands = TRUE
)

# 2. Process by regions
us_states <- c("Ohio", "Michigan", "Indiana") 
state_results <- list()

for (state in us_states) {
  state_results[[state]] <- calculate_vegetation_index(
    spectral_data = "continental_satellite_data.tif",
    index_type = "NDVI",
    region_boundary = state,     # Process each state separately
    auto_detect_bands = TRUE
  )
}

# 3. Combine results
combined_results <- create_raster_mosaic(state_results, method = "merge")
```

## Performance Optimization

### Best Practices

1. **Apply region boundaries early** to reduce data size
2. **Use appropriate resolution** for your analysis scale
3. **Batch process** multiple indices together
4. **Cache intermediate results** for complex workflows

```{r eval=FALSE}
# Efficient workflow
optimized_workflow <- function(satellite_data, study_region) {
  # 1. Crop to region first (reduces data size)
  cropped_data <- universal_spatial_join(
    source_data = satellite_data,
    target_data = get_region_boundary(study_region),
    method = "crop"
  )
  
  # 2. Calculate multiple indices in one call
  vegetation_indices <- calculate_multiple_indices(
    spectral_data = cropped_data,
    indices = c("NDVI", "EVI", "SAVI", "GNDVI"),
    auto_detect_bands = TRUE,
    output_stack = TRUE,
    parallel = FALSE  # Use FALSE for stability
  )
  
  # 3. Cache results
  terra::writeRaster(vegetation_indices, "cached_vegetation_indices.tif")
  
  return(vegetation_indices)
}
```

## Interpretation Guidelines

### NDVI Interpretation by Land Cover

**Cropland:**
- 0.2-0.4: Bare soil/early growth
- 0.4-0.6: Developing vegetation  
- 0.6-0.8: Healthy, mature crops
- 0.8-0.9: Peak vegetation (dense crops)

**Forest:**
- 0.4-0.6: Sparse forest/stressed trees
- 0.6-0.8: Healthy forest
- 0.8-0.9: Dense, healthy forest

**Grassland:**
- 0.2-0.4: Sparse grass/dormant season
- 0.4-0.6: Moderate grass cover
- 0.6-0.8: Dense, healthy grassland

### Seasonal Patterns

**Temperate Crops (Corn/Soybeans):**
- **April-May:** 0.2-0.4 (emergence)
- **June-July:** 0.4-0.7 (vegetative growth)
- **July-August:** 0.6-0.9 (peak season)
- **September-October:** 0.4-0.6 (senescence)

## Advanced Topics

### Custom Index Development

Create your own vegetation indices:

```{r eval=FALSE}
# Custom ratio index
custom_ratio <- nir_band / red_band
names(custom_ratio) <- "Custom_Ratio"

# Custom normalized difference
custom_nd <- (nir_band - green_band) / (nir_band + green_band)
names(custom_nd) <- "Custom_ND"

# Combine with standard indices
all_indices <- c(
  calculate_vegetation_index(red = red, nir = nir, index_type = "NDVI"),
  custom_ratio,
  custom_nd
)
```

### Integration with External Data

```{r eval=FALSE}
# Combine vegetation indices with environmental data
environmental_analysis <- universal_spatial_join(
  source_data = vegetation_indices,
  target_data = list(
    precipitation = "annual_precip.tif",
    temperature = "mean_temp.tif", 
    elevation = "dem.tif"
  ),
  method = "extract"
)

# Analyze relationships
vegetation_env_correlation <- analyze_variable_correlations(
  variable_list = list(
    NDVI = vegetation_indices$NDVI,
    precipitation = environmental_data$precipitation,
    temperature = environmental_data$temperature
  ),
  create_plots = TRUE
)
```

## Conclusion

This vignette covered comprehensive vegetation analysis with GeoSpatialSuite:

- **60+ vegetation indices** for diverse applications
- **Time series analysis** for monitoring changes
- **Quality control** for reliable results  
- **Agricultural applications** for crop monitoring
- **Stress detection** for early warning systems
- **Multi-scale analysis** from field to regional scales


### Additional Resources

```{r eval=FALSE}
# Explore more indices
list_vegetation_indices(detailed = TRUE)

# Get help
?calculate_vegetation_index
?calculate_multiple_indices  
?analyze_crop_vegetation

# Test your installation
test_geospatialsuite_package_simple(verbose = TRUE)
```
## Acknowledgments

This work was developed by the GeoSpatialSuite team with contributions from:
Olatunde D. Akanbi, Erika I. Barcelos, and Roger H. French.