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

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

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

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

# 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

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

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

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

# 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

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

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

ndvi <- calculate_vegetation_index(red = red, nir = nir, index_type = "NDVI")

evi <- calculate_vegetation_index(red = red, nir = nir, blue = blue, index_type = "EVI")

savi <- calculate_vegetation_index(red = red, nir = nir, index_type = "SAVI")

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

crop_health <- calculate_multiple_indices(
  red = red, nir = nir, green = green,
  indices = c("NDVI", "GNDVI", "SIPI"),  # Structure Insensitive Pigment Index
  output_stack = TRUE
)

drought_indices <- calculate_multiple_indices(
  nir = nir, swir1 = swir1,
  indices = c("NDMI", "MSI"),  # Moisture indices
  output_stack = TRUE
)

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
)

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

# 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

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

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

# 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), "%"))

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

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

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

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

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

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

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

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