---
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



## Complete Index Reference

All 62 vegetation indices with formulas, ranges, and references:

```{r complete-index-table, echo=FALSE, message=FALSE}
library(geospatialsuite)
detailed <- list_vegetation_indices(detailed = TRUE)

knitr::kable(
  detailed[, c("Index", "Category", "Formula", "Range", "Reference")],
  caption = "Complete Vegetation Indices Reference - All 62 Indices",
  format = "html"
)
```

### Using Index Information
```{r index-info-example, eval=FALSE}
# Get all indices with formulas
all_indices <- list_vegetation_indices(detailed = TRUE)
View(all_indices)

# Filter by category
stress <- list_vegetation_indices(category = "stress", detailed = TRUE)
water <- list_vegetation_indices(application = "water", detailed = TRUE)

# Get specific index information
ndvi_info <- all_indices[all_indices$Index == "NDVI", ]
cat(sprintf("NDVI Formula: %s\n", ndvi_info$Formula))
cat(sprintf("NDVI Range: %s\n", ndvi_info$Range))
cat(sprintf("Reference: %s\n", ndvi_info$Reference))
```








## 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)
```


# Band Naming Conventions for geospatialsuite

## Overview
The `calculate_vegetation_index()` and `analyze_crop_vegetation()` functions support automatic band detection from multiple satellite platforms. This document explains how band naming works and what formats are supported.

---

## How Band Detection Works

### 1. Case-Insensitive Matching
**All band names are case-insensitive.**

✅ These are all equivalent:
- `"red"`, `"RED"`, `"Red"`, `"rEd"`
- `"nir"`, `"NIR"`, `"Nir"`, `"nIR"`
- `"B4"`, `"b4"`, `"B04"`, `"b04"`

---

## Supported Band Names by Type

### Red Band
Recognized patterns:
- Generic: `red`, `RED`, `Red`
- Landsat 8/9: `B4`, `b4`, `band4`, `Band_4`, `Band4`
- Sentinel-2: `B04`, `b04`
- Generic numbered: `band_4`, `sr_band4`

**Example:**
```{r eval=FALSE}
# All of these will be recognized as the red band
calculate_vegetation_index(red = red_band, nir = nir_band, index_type = "NDVI")
calculate_vegetation_index(spectral_data = raster_with_band_named_"B4", auto_detect_bands = TRUE)
calculate_vegetation_index(spectral_data = raster_with_band_named_"RED", auto_detect_bands = TRUE)
```

### NIR (Near Infrared) Band
Recognized patterns:
- Generic: `nir`, `NIR`, `Nir`, `near_infrared`
- Landsat 8/9: `B5`, `b5`, `band5`, `Band_5`, `Band5`
- Sentinel-2: `B08`, `b08`, `B8`, `b8`, `band8`, `Band_8`, `Band8`

### Blue Band
Recognized patterns:
- Generic: `blue`, `BLUE`, `Blue`
- Landsat 8/9: `B2`, `b2`, `band2`, `Band_2`, `Band2`
- Sentinel-2: `B02`, `b02`

### Green Band
Recognized patterns:
- Generic: `green`, `GREEN`, `Green`
- Landsat 8/9: `B3`, `b3`, `band3`, `Band_3`, `Band3`
- Sentinel-2: `B03`, `b03`

### SWIR1 (Shortwave Infrared 1) Band
Recognized patterns:
- Generic: `swir1`, `SWIR1`, `Swir1`, `shortwave_infrared_1`, `SWIR_1`
- Landsat 8/9: `B6`, `b6`
- Sentinel-2: `B11`, `b11`

### SWIR2 (Shortwave Infrared 2) Band
Recognized patterns:
- Generic: `swir2`, `SWIR2`, `Swir2`, `shortwave_infrared_2`, `SWIR_2`
- Landsat 8/9: `B7`, `b7`
- Sentinel-2: `B12`, `b12`

### Red Edge Band
Recognized patterns:
- Generic: `rededge`, `RedEdge`, `red_edge`, `RED_EDGE`, `RE`, `re`
- Sentinel-2: `B05`, `b05`, `B5`, `b5` (note: conflicts with Landsat NIR)

### Coastal/Aerosol Band
Recognized patterns:
- Generic: `coastal`, `COASTAL`, `Coastal`, `aerosol`
- Landsat 8/9: `B1`, `b1`
- Sentinel-2: `B01`, `b01`

---

## Satellite-Specific Examples

### Landsat 8/9 (OLI)

**Band Mapping:**
| Band | Number | Wavelength | geospatialsuite Name |
|------|--------|------------|---------------------|
| Coastal/Aerosol | B1 | 0.43-0.45 µm | `coastal`, `B1` |
| Blue | B2 | 0.45-0.51 µm | `blue`, `B2` |
| Green | B3 | 0.53-0.59 µm | `green`, `B3` |
| Red | B4 | 0.64-0.67 µm | `red`, `B4` |
| NIR | B5 | 0.85-0.88 µm | `nir`, `B5` |
| SWIR1 | B6 | 1.57-1.65 µm | `swir1`, `B6` |
| SWIR2 | B7 | 2.11-2.29 µm | `swir2`, `B7` |

**Example - Landsat Scene:**
```{r eval=FALSE}
library(geospatialsuite)
library(terra)

# If your Landsat file has band names: B1, B2, B3, B4, B5, B6, B7
landsat <- rast("LC08_L2SP_029030_20230615_B1-7.tif")

# Auto-detection will work automatically
ndvi <- calculate_vegetation_index(
  spectral_data = landsat,
  index_type = "NDVI",
  auto_detect_bands = TRUE
)

# Or specify bands explicitly
ndvi <- calculate_vegetation_index(
  red = landsat[["B4"]],
  nir = landsat[["B5"]],
  index_type = "NDVI"
)
```

### Sentinel-2 (MSI)

**Band Mapping:**
| Band | Number | Wavelength | Resolution | geospatialsuite Name |
|------|--------|------------|-----------|---------------------|
| Coastal | B01 | 0.433-0.453 µm | 60m | `coastal`, `B01`, `B1` |
| Blue | B02 | 0.458-0.523 µm | 10m | `blue`, `B02`, `B2` |
| Green | B03 | 0.543-0.578 µm | 10m | `green`, `B03`, `B3` |
| Red | B04 | 0.650-0.680 µm | 10m | `red`, `B04`, `B4` |
| Red Edge 1 | B05 | 0.698-0.713 µm | 20m | `red_edge`, `B05`, `B5` |
| Red Edge 2 | B06 | 0.733-0.748 µm | 20m | N/A |
| Red Edge 3 | B07 | 0.765-0.785 µm | 20m | N/A |
| NIR | B08 | 0.785-0.900 µm | 10m | `nir`, `B08`, `B8` |
| SWIR1 | B11 | 1.565-1.655 µm | 20m | `swir1`, `B11` |
| SWIR2 | B12 | 2.100-2.280 µm | 20m | `swir2`, `B12` |

**Example - Sentinel-2 Scene:**
```{r eval=FALSE}
# If your Sentinel-2 file has standard band names
sentinel <- rast("S2A_MSIL2A_20230615_B02-B12.tif")

# Auto-detection works
ndvi <- calculate_vegetation_index(
  spectral_data = sentinel,
  index_type = "NDVI",
  auto_detect_bands = TRUE
)

# Calculate red edge index (requires Sentinel-2)
ndre <- calculate_vegetation_index(
  spectral_data = sentinel,
  index_type = "NDRE",
  auto_detect_bands = TRUE
)

# Or use specific band names
evi <- calculate_vegetation_index(
  red = sentinel[["B04"]],
  nir = sentinel[["B08"]],
  blue = sentinel[["B02"]],
  index_type = "EVI"
)
```

### MODIS (Terra/Aqua)

**Band Mapping:**
| Band | Number | Wavelength | geospatialsuite Name |
|------|--------|------------|---------------------|
| Red | 1 | 0.620-0.670 µm | `red`, `band1` |
| NIR | 2 | 0.841-0.876 µm | `nir`, `band2` |
| Blue | 3 | 0.459-0.479 µm | `blue`, `band3` |
| Green | 4 | 0.545-0.565 µm | `green`, `band4` |
| SWIR1 | 6 | 1.628-1.652 µm | `swir1`, `band6` |
| SWIR2 | 7 | 2.105-2.155 µm | `swir2`, `band7` |

**Example - MODIS:**
```{r eval=FALSE}
modis <- rast("MOD13Q1_NDVI_2023165.tif")

# If bands are named generically
ndvi <- calculate_vegetation_index(
  spectral_data = modis,
  index_type = "NDVI",
  auto_detect_bands = TRUE
)
```

---

## Custom Band Names

If your raster has non-standard band names, you have three options:

### Option 1: Rename Bands
```{r eval=FALSE}
# Your raster has bands: "band_A", "band_B", "band_C", "band_D"
my_raster <- rast("custom_satellite.tif")

# Rename to standard names
names(my_raster) <- c("red", "green", "blue", "nir")

# Now auto-detection works
ndvi <- calculate_vegetation_index(
  spectral_data = my_raster,
  index_type = "NDVI",
  auto_detect_bands = TRUE
)
```

### Option 2: Specify Band Names
```{r eval=FALSE}
# Tell the function which bands are which
ndvi <- calculate_vegetation_index(
  spectral_data = my_raster,
  band_names = c("band_A", "band_B", "band_C", "band_D"),  # Order: Red, Green, Blue, NIR
  index_type = "NDVI",
  auto_detect_bands = FALSE
)
```

### Option 3: Extract and Pass Bands Directly
```{r eval=FALSE}
# Extract specific layers
red_band <- my_raster[[1]]   # Assuming layer 1 is red
nir_band <- my_raster[[4]]   # Assuming layer 4 is NIR

# Pass explicitly
ndvi <- calculate_vegetation_index(
  red = red_band,
  nir = nir_band,
  index_type = "NDVI"
)
```

---

## Common Issues and Solutions

### Issue 1: "Band not found" error
**Problem:** Your band names don't match any recognized patterns

**Solution:**
```{r eval=FALSE}
# Check your band names
names(my_raster)
# [1] "sr_band_4" "sr_band_5" "sr_band_3" "sr_band_2"

# Rename them
names(my_raster) <- c("red", "nir", "green", "blue")

# Or extract and pass explicitly
calculate_vegetation_index(
  red = my_raster[[1]],
  nir = my_raster[[2]],
  index_type = "NDVI"
)
```

### Issue 2: Wrong bands detected
**Problem:** Auto-detection picked the wrong bands (e.g., Sentinel-2 B05 detected as NIR instead of Red Edge)

**Solution:**
```{r eval=FALSE}
# Don't use auto-detection, be explicit
calculate_vegetation_index(
  red = my_raster[["B04"]],
  nir = my_raster[["B08"]],
  red_edge = my_raster[["B05"]],
  index_type = "NDRE"
)
```

### Issue 3: Multiple files/directories
**Problem:** Your bands are in separate files

**Solution A - List of files:**
```{r eval=FALSE}
# Create file list
band_files <- c(
  red = "LC08_B4.tif",
  nir = "LC08_B5.tif",
  blue = "LC08_B2.tif"
)

# Let the function load them
evi <- calculate_vegetation_index(
  spectral_data = band_files,
  index_type = "EVI"
)
```

**Solution B - Directory:**
```{r eval=FALSE}
# If all bands are in one directory with recognizable names
ndvi <- calculate_vegetation_index(
  spectral_data = "/path/to/landsat/bands/",
  index_type = "NDVI",
  auto_detect_bands = TRUE
)
```

---

## Band Detection Priority

When using auto-detection, the function searches in this order:

1. **Exact matches** (case-insensitive)
   - `"red"` matches `"red"`, `"RED"`, `"Red"`
   
2. **Satellite-specific patterns**
   - `"B4"` for Landsat red
   - `"B04"` for Sentinel-2 red
   
3. **Generic patterns**
   - `"band4"`, `"Band_4"`, etc.
   
4. **Partial matches** (regex, case-insensitive)
   - Any layer with "red" in the name
   
5. **Position-based fallback** (if all else fails)
   - Band 1 → Red
   - Band 2 → Green  
   - Band 3 → Blue
   - Band 4 → NIR

---

## Testing Band Detection

To see which bands were detected:

```{r eval=FALSE}
# Turn on verbose mode
result <- calculate_vegetation_index(
  spectral_data = my_raster,
  index_type = "NDVI",
  auto_detect_bands = TRUE,
  verbose = TRUE  # This will print which bands were detected
)
```

Expected output:

```
Processing spectral_data input...
Extracting bands from multi-band raster...
Exact match detected red band: B4
Exact match detected nir band: B8
Starting NDVI calculation with comprehensive input handling...
```

---

## Best Practices

1. **Use standard naming** when possible (`red`, `nir`, `blue`, etc.)
2. **Keep satellite band numbers** if using Landsat/Sentinel-2 (B4, B05, etc.)
3. **Enable verbose mode** when debugging band detection issues
4. **Rename bands** for clarity rather than relying on position-based fallback
5. **Be explicit** when working with unusual sensors or custom band combinations

---

## Summary Table

| Your Data | Band Detection Method | Example |
|-----------|----------------------|---------|
| Landsat with standard names | Auto-detect | `auto_detect_bands = TRUE` |
| Sentinel-2 with standard names | Auto-detect | `auto_detect_bands = TRUE` |
| Generic names (red, nir, blue) | Auto-detect | `auto_detect_bands = TRUE` |
| Custom names | Rename first | `names(raster) <- c("red", "nir", ...)` |
| Separate files | Pass list | `spectral_data = c(red="file1.tif", ...)` |
| Non-standard structure | Extract & pass | `red = raster[[1]], nir = raster[[4]]` |

---











## Acknowledgments

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