Vegetation Index Analysis with GeoSpatialSuite

GeoSpatialSuite Development Team

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:

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:

Complete Vegetation Indices Reference - All 62 Indices
Index Category Formula Range Reference
NDVI basic (NIR - Red) / (NIR + Red) [-1, 1] Rouse et al. (1974)
SAVI basic ((NIR - Red) / (NIR + Red + L)) * (1 + L), where L=0.5 [-1, 1.5] Huete (1988)
MSAVI basic 0.5 * (2NIR + 1 - sqrt((2NIR + 1)^2 - 8*(NIR - Red))) [0, 2] Qi et al. (1994)
OSAVI basic (NIR - Red) / (NIR + Red + 0.16) [-1, 1] Rondeaux et al. (1996)
EVI basic 2.5 * ((NIR - Red) / (NIR + 6Red - 7.5Blue + 1)) [-1, 3] Huete et al. (1997)
EVI2 basic 2.5 * ((NIR - Red) / (NIR + 2.4*Red + 1)) [-1, 3] Jiang et al. (2008)
DVI basic NIR - Red [-2, 2] Richardson & Wiegand (1977)
RVI basic NIR / Red [0, 30] Birth & McVey (1968)
GNDVI basic (NIR - Green) / (NIR + Green) [-1, 1] Gitelson et al. (1996)
WDVI basic NIR - 0.5 * Red [-2, 2] Clevers (1988)
ARVI enhanced (NIR - (2Red - Blue)) / (NIR + (2Red - Blue)) [-1, 1] Kaufman & Tanre (1992)
RDVI enhanced (NIR - Red) / sqrt(NIR + Red) [-2, 2] Roujean & Breon (1995)
PVI enhanced (NIR - a*Red - b) / sqrt(1 + a^2), where a=1.5, b=10 [-2, 2] Richardson & Wiegand (1977)
IPVI enhanced NIR / (NIR + Red) [0, 1] Crippen (1990)
TNDVI enhanced sqrt(((NIR - Red) / (NIR + Red)) + 0.5) [0, 1.2] Deering et al. (1975)
GEMI enhanced eta * (1 - 0.25eta) - (Red - 0.125) / (1 - Red), where eta = (2(NIR^2 - Red^2) + 1.5NIR + 0.5Red) / (NIR + Red + 0.5) [-1, 1] Pinty & Verstraete (1992)
VARI enhanced (Green - Red) / (Green + Red - Blue) [-1, 1] Gitelson et al. (2002)
TSAVI enhanced (a * (NIR - aRed - b)) / (Red + aNIR - ab + X(1 + a^2)), where a=1.5, b=10, X=0.08 [-1, 1.5] Baret & Guyot (1991)
ATSAVI enhanced (a * (NIR - aRed - b)) / (aNIR + Red - ab + X(1 + a^2)), where a=1.22, b=0.03, X=0.08 [-1, 1.5] Baret et al. (1992)
GESAVI enhanced ((NIR - Red) / (NIR + Red + Z)) * (1 + Z), where Z=0.5 [-1, 1.5] Gilabert et al. (2002)
MTVI enhanced 1.2 * (1.2(NIR - Green) - 2.5(Red - Green)) [-2, 2] Haboudane et al. (2004)
CTVI enhanced ((NDVI + 0.5) / |NDVI + 0.5|) * sqrt(|NDVI + 0.5|), where NDVI = (NIR - Red) / (NIR + Red) [-1, 1.5] Perry & Lautenschlager (1984)
NDRE advanced (NIR - RedEdge) / (NIR + RedEdge) [-1, 1] Gitelson & Merzlyak (1994)
MTCI advanced (RedEdge - Red) / (NIR - Red) [0, 8] Dash & Curran (2004)
IRECI advanced (RedEdge - Red) / (RedEdge / NIR) [-2, 5] Frampton et al. (2013)
S2REP advanced 705 + 35 * ((Red + RedEdge)/2 - Red) / (RedEdge - Red) [680, 780] Frampton et al. (2013)
PSRI advanced (Red - Green) / RedEdge [-1, 1] Merzlyak et al. (1999)
CRI1 advanced (1 / Green) - (1 / Red) [-10, 10] Gitelson et al. (2002)
CRI2 advanced (1 / Green) - (1 / RedEdge) [-10, 10] Gitelson et al. (2002)
ARI1 advanced (1 / Green) - (1 / RedEdge) [-10, 10] Gitelson et al. (2001)
ARI2 advanced NIR * ((1 / Green) - (1 / RedEdge)) [-10, 10] Gitelson et al. (2001)
MCARI advanced ((RedEdge - Red) - 0.2(RedEdge - Green)) (RedEdge / Red) [-2, 2] Daughtry et al. (2000)
PRI stress (Green - NIR) / (Green + NIR) [-1, 1] Gamon et al. (1992)
SIPI stress (NIR - Red) / (NIR - Green) [0, 2] Penuelas et al. (1995)
CCI stress (RedEdge - Red) / (RedEdge + Red) [-1, 1] Barnes et al. (2000)
NDNI stress (log(1/NIR) - log(1/SWIR1)) / (log(1/NIR) + log(1/SWIR1)) [-1, 1] Serrano et al. (2002)
CARI stress RedEdge * (Red / Green) [0, 5] Kim et al. (1994)
TCARI stress 3 * ((RedEdge - Red) - 0.2(RedEdge - Green) (RedEdge / Red)) [-2, 5] Haboudane et al. (2002)
MTVI1 stress 1.2 * (1.2(NIR - Green) - 2.5(Red - Green)) [-2, 2] Haboudane et al. (2004)
MTVI2 stress 1.5 * (1.2(NIR - Green) - 2.5(Red - Green)) / sqrt((2NIR + 1)^2 - (6NIR - 5*sqrt(Red)) - 0.5) [-2, 5] Haboudane et al. (2004)
TVI stress 0.5 * (120(NIR - Green) - 200(Red - Green)) [-100, 100] Broge & Leblanc (2001)
NPCI stress (Red - Blue) / (Red + Blue) [-1, 1] Penuelas et al. (1994)
RARS stress Red / NIR [0, 5] Chappelle et al. (1992)
NPQI stress (Red - Blue) / (Red + Blue) [-1, 1] Barnes et al. (1992)
NDWI water (Green - NIR) / (Green + NIR) [-1, 1] McFeeters (1996)
MNDWI water (Green - SWIR1) / (Green + SWIR1) [-1, 1] Xu (2006)
NDMI water (NIR - SWIR1) / (NIR + SWIR1) [-1, 1] Gao (1996)
MSI water SWIR1 / NIR [0, 5] Hunt & Rock (1989)
NDII water (NIR - SWIR1) / (NIR + SWIR1) [-1, 1] Hardisky et al. (1983)
WI water NIR / SWIR1 [0, 10] Gao (1996)
SRWI water NIR / SWIR1 [0, 10] Zarco-Tejada et al. (2003)
LSWI water (NIR - SWIR1) / (NIR + SWIR1) [-1, 1] Xiao et al. (2002)
LAI specialized 3.618 * EVI - 0.118, where EVI = 2.5 * ((NIR - Red) / (NIR + 6Red - 7.5Blue + 1)) [0, 15] Baret & Guyot (1991)
FAPAR specialized -0.161 + 1.257 * NDVI, where NDVI = (NIR - Red) / (NIR + Red) [0, 1] Myneni & Williams (1994)
FCOVER specialized -2.274 + 4.336NDVI - 1.33NDVI^2, where NDVI = (NIR - Red) / (NIR + Red) [0, 1] Baret et al. (2007)
NBR specialized (NIR - SWIR2) / (NIR + SWIR2) [-1, 1] Lopez Garcia & Caselles (1991)
BAI specialized 1 / ((0.1 - Red)^2 + (0.06 - NIR)^2) [0, 1000] Chuvieco et al. (2002)
NDSI specialized (Green - SWIR1) / (Green + SWIR1) [-1, 1] Hall et al. (1995)
GRVI specialized (Green - Red) / (Green + Red) [-1, 1] Tucker (1979)
VIG specialized (Green - Red) / (Green + Red) [-1, 1] Gitelson et al. (2002)
CI specialized (Red - Green) / Red [-1, 1] Escadafal & Huete (1991)
GBNDVI specialized (NIR - (Green + Blue)) / (NIR + Green + Blue) [-1, 1] Wang et al. (2007)

Using Index Information

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

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

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

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

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

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

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

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

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

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

# 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

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

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

# 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

# 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

# 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

# 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

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

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

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

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

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:

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

Precision Agriculture:

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

# 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

# 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

# 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

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

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

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

# 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

# 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

# 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

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

# 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

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

Additional Resources

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

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

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:

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

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

# 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

# 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

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

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

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

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

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

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