Nothing
R/02-vegetation-indices.R
In geospatialsuite: Comprehensive Geospatiotemporal Analysis and Multimodal Integration Toolkit
Defines functions
get_index_references
get_index_ranges
get_index_formulas
calculate_vegetation_statistics
analyze_yield_potential
analyze_growth_stage
detect_vegetation_stress
select_crop_indices
validate_output
mask_invalid_values
apply_temporal_smoothing
apply_quality_filter
handle_index_edge_cases
calculate_index_by_type
validate_required_bands
load_and_validate_band
auto_detect_spectral_bands
extract_bands_from_raster
get_available_indices
load_and_stack_bands
load_bands_from_directory
analyze_crop_vegetation
calculate_multiple_indices
list_vegetation_indices
calculate_ndvi_enhanced
calculate_vegetation_index
Documented in
analyze_crop_vegetation
analyze_growth_stage
analyze_yield_potential
apply_quality_filter
apply_temporal_smoothing
auto_detect_spectral_bands
calculate_index_by_type
calculate_multiple_indices
calculate_ndvi_enhanced
calculate_vegetation_index
calculate_vegetation_statistics
detect_vegetation_stress
extract_bands_from_raster
get_available_indices
get_index_formulas
get_index_ranges
get_index_references
handle_index_edge_cases
list_vegetation_indices
load_and_stack_bands
load_and_validate_band
load_bands_from_directory
mask_invalid_values
select_crop_indices
validate_output
validate_required_bands
#' Calculate comprehensive vegetation indices
#'
#' @description
#' Calculate a wide range of vegetation indices from spectral bands with
#' automatic band detection, comprehensive error handling, and validation.
#' Supports 40+ different vegetation indices for various applications.
#' Accepts directories, file lists, and automatic CRS handling.
#'
#' @param spectral_data Either individual bands (red, nir, etc.), a multi-band raster,
#' directory path, or list of raster files
#' @param red Red band SpatRaster or file path
#' @param nir NIR band SpatRaster or file path
#' @param blue Optional blue band
#' @param green Optional green band
#' @param swir1 Optional SWIR1 band
#' @param swir2 Optional SWIR2 band
#' @param red_edge Optional Red Edge band
#' @param coastal Optional Coastal/Aerosol band
#' @param nir2 Optional second NIR band
#' @param index_type Vegetation index to calculate (see list_vegetation_indices())
#' @param auto_detect_bands Automatically detect bands from multi-band raster
#' @param band_names Custom band names for multi-band input
#' @param clamp_range Range to clamp output values (optional)
#' @param mask_invalid Mask invalid/extreme values
#' @param scale_factor Scaling factor if needed (default: 1)
#' @param auto_crs_fix Automatically fix CRS mismatches between bands
#' @param verbose Print progress messages
#'
#' @return SpatRaster of vegetation index
#'
#' @details
#' ## Input Format Support:
#'
#' ### Single Calculation:
#' ```r
#' # Individual band files
#' ndvi <- calculate_vegetation_index(red = "red.tif", nir = "nir.tif", index_type = "NDVI")
#'
#' # Multi-band raster
#' evi <- calculate_vegetation_index(spectral_data = "landsat.tif", index_type = "EVI",
#' auto_detect_bands = TRUE)
#' ```
#'
#' ### Directory/Multiple Files:
#' ```r
#' # Directory with band files
#' savi <- calculate_vegetation_index(spectral_data = "/path/to/bands/",
#' band_names = c("red", "nir"), index_type = "SAVI")
#'
#' # File list
#' arvi <- calculate_vegetation_index(spectral_data = c("red.tif", "nir.tif", "blue.tif"),
#' index_type = "ARVI")
#' ```
#'
#' ## Enhanced vs Basic NDVI:
#'
#' ### Basic `calculate_vegetation_index()`:
#' - Single time point calculation
#' - 40+ different indices
#' - Directory/file support
#' - Automatic CRS fixing
#' - **Use for**: Single-date analysis, comparing different indices
#'
#' ### `calculate_ndvi_enhanced()`:
#' - Time series support
#' - Quality filtering
#' - Temporal smoothing
#' - Date matching between red/NIR
#' - **Use for**: Multi-temporal analysis, time series trends
#'
#'#' ## Band Naming Conventions:
#'
#' The function supports case-insensitive band detection:
#' - **Generic names**: "red"/"RED"/"Red", "nir"/"NIR", "blue"/"BLUE", "green"/"GREEN"
#' - **Landsat 8/9**: B1-B7 (e.g., B4=Red, B5=NIR)
#' - **Sentinel-2**: B01-B12 (e.g., B04=Red, B08=NIR, B05=RedEdge)
#' - **MODIS**: band1-band7
#'
#' ### Satellite-Specific Examples:
#'
#' **Landsat 8/9:**
#' ```r
#' # Bands automatically detected
#' ndvi <- calculate_vegetation_index(
#' spectral_data = "LC08_stack.tif", # Has bands named B1-B7
#' index_type = "NDVI",
#' auto_detect_bands = TRUE
#' )
#' ```
#'
#' **Sentinel-2:**
#' ```r
#' # Red Edge indices need Sentinel-2
#' ndre <- calculate_vegetation_index(
#' spectral_data = sentinel_data, # Has B01-B12
#' index_type = "NDRE",
#' auto_detect_bands = TRUE
#' )
#' ```
#'
#' **Custom band names:**
#' ```r
#' # Rename your bands first
#' names(my_raster) <- c("red", "nir", "blue", "green")
#'
#' # Or specify explicitly
#' ndvi <- calculate_vegetation_index(
#' red = my_raster[[1]],
#' nir = my_raster[[4]],
#' index_type = "NDVI"
#' )
#' ```
#'
#' For complete band naming documentation, see:
#' \code{vignette("vegetation-indices", package = "geospatialsuite")}
#'
#' @examples
#' \dontrun{
#' # These examples require satellite imagery files (Landsat/Sentinel data etc.)
#' # Basic NDVI calculation
#' ndvi <- calculate_vegetation_index(red = red_band, nir = nir_band, index_type = "NDVI")
#'
#' # Multi-band raster with auto-detection
#' evi <- calculate_vegetation_index(spectral_data = landsat_stack,
#' index_type = "EVI", auto_detect_bands = TRUE)
#'
#' # Directory with automatic band detection
#' savi <- calculate_vegetation_index(spectral_data = "/path/to/sentinel/bands/",
#' index_type = "SAVI", auto_detect_bands = TRUE)
#'
#' # Advanced index with custom parameters
#' pri <- calculate_vegetation_index(red = red_band, nir = nir_band, green = green_band,
#' index_type = "PRI", clamp_range = c(-1, 1))
#'
#' # Custom band names for multi-band data
#' ndvi <- calculate_vegetation_index(spectral_data = sentinel_data,
#' band_names = c("B4", "B3", "B2", "B8"),
#' index_type = "NDVI")
#' }
#'
#' @export
calculate_vegetation_index <- function(spectral_data = NULL, red = NULL, nir = NULL,
blue = NULL, green = NULL, swir1 = NULL, swir2 = NULL,
red_edge = NULL, coastal = NULL, nir2 = NULL,
index_type = "NDVI", auto_detect_bands = FALSE,
band_names = NULL, clamp_range = NULL,
mask_invalid = TRUE, scale_factor = 1,
auto_crs_fix = TRUE, verbose = FALSE) {
if (verbose) message(sprintf("Starting %s calculation with comprehensive input handling...", index_type))
# Validate index type
available_indices <- get_available_indices()
if (!index_type %in% available_indices$Index) {
stop(sprintf("Index '%s' not supported. Use list_vegetation_indices() to see available indices.",
index_type), call. = FALSE)
}
# Handle input processing
if (!is.null(spectral_data)) {
if (verbose) message("Processing spectral_data input...")
if (is.character(spectral_data)) {
if (length(spectral_data) == 1) {
if (dir.exists(spectral_data)) {
# Directory provided - load and combine bands
if (verbose) message(sprintf("Loading bands from directory: %s", spectral_data))
spectral_data <- load_bands_from_directory(spectral_data, band_names, verbose)
} else if (file.exists(spectral_data)) {
# Single file provided
if (verbose) message(sprintf("Loading multi-band raster: %s", basename(spectral_data)))
spectral_data <- terra::rast(spectral_data)
} else {
stop(sprintf("Spectral data path does not exist: %s", spectral_data), call. = FALSE)
}
} else {
# Multiple files provided - stack them
if (verbose) message(sprintf("Loading and stacking %d band files", length(spectral_data)))
spectral_data <- load_and_stack_bands(spectral_data, band_names, verbose)
}
}
# Extract bands from processed spectral data
bands <- extract_bands_from_raster(spectral_data, auto_detect_bands, band_names, verbose)
red <- bands$red
nir <- bands$nir
blue <- bands$blue
green <- bands$green
swir1 <- bands$swir1
swir2 <- bands$swir2
red_edge <- bands$red_edge
coastal <- bands$coastal
nir2 <- bands$nir2
}
# Input validation and loading with CRS checking
if (verbose) message("Loading and validating input bands with CRS checking...")
red <- load_and_validate_band(red, "red", required = TRUE)
nir <- load_and_validate_band(nir, "nir", required = TRUE)
# CRS compatibility checking
if (verbose) message("Checking and fixing spatial compatibility...")
reference_crs <- terra::crs(red)
# Check and fix NIR CRS
if (!identical(terra::crs(nir), reference_crs)) {
if (auto_crs_fix) {
if (verbose) message("CRS mismatch detected for NIR band. Reprojecting...")
nir <- terra::project(nir, reference_crs)
} else {
stop("CRS mismatch between red and NIR bands. Set auto_crs_fix=TRUE to fix automatically.", call. = FALSE)
}
}
# Align geometries if needed
if (!terra::compareGeom(red, nir, stopOnError = FALSE)) {
if (verbose) message("Geometry mismatch detected. Resampling NIR to match red...")
nir <- terra::resample(nir, red)
}
# Load and align optional bands with CRS fixing
bands_to_process <- list(
blue = blue, green = green, swir1 = swir1, swir2 = swir2,
red_edge = red_edge, coastal = coastal, nir2 = nir2
)
processed_bands <- list()
for (band_name in names(bands_to_process)) {
band_data <- bands_to_process[[band_name]]
if (!is.null(band_data)) {
band_raster <- load_and_validate_band(band_data, band_name, required = FALSE)
if (!is.null(band_raster)) {
# Fix CRS if needed
if (!identical(terra::crs(band_raster), reference_crs) && auto_crs_fix) {
if (verbose) message(sprintf("Reprojecting %s band to match reference CRS", band_name))
band_raster <- terra::project(band_raster, reference_crs)
}
# Align geometry
band_raster <- check_raster_compatibility(red, band_raster, auto_align = TRUE)
processed_bands[[band_name]] <- band_raster
}
}
}
# Assign processed bands
blue <- processed_bands$blue
green <- processed_bands$green
swir1 <- processed_bands$swir1
swir2 <- processed_bands$swir2
red_edge <- processed_bands$red_edge
coastal <- processed_bands$coastal
nir2 <- processed_bands$nir2
# Validate required bands for specific indices
if (!validate_required_bands(index_type, blue, green, swir1, swir2, red_edge, coastal, nir2, verbose)) {
if (verbose) message(sprintf("Skipping %s calculation due to missing required bands", index_type))
# Return NA raster with same dimensions as red band
na_raster <- red
terra::values(na_raster) <- NA
names(na_raster) <- index_type
return(na_raster)
}
# Apply scale factor if needed
if (scale_factor != 1) {
if (verbose) message(sprintf("Applying scale factor: %g", scale_factor))
red <- red * scale_factor
nir <- nir * scale_factor
if (!is.null(blue)) blue <- blue * scale_factor
if (!is.null(green)) green <- green * scale_factor
if (!is.null(swir1)) swir1 <- swir1 * scale_factor
if (!is.null(swir2)) swir2 <- swir2 * scale_factor
if (!is.null(red_edge)) red_edge <- red_edge * scale_factor
if (!is.null(coastal)) coastal <- coastal * scale_factor
if (!is.null(nir2)) nir2 <- nir2 * scale_factor
}
# Calculate vegetation index
if (verbose) message(sprintf("Computing %s index...", index_type))
index <- calculate_index_by_type(index_type, red, nir, blue, green, swir1, swir2,
red_edge, coastal, nir2, verbose)
# Apply clamp range if specified
if (!is.null(clamp_range)) {
if (verbose) message(sprintf("Clamping values to range [%.3f, %.3f]",
clamp_range[1], clamp_range[2]))
index <- terra::clamp(index, clamp_range[1], clamp_range[2])
}
# Mask invalid values if requested
if (mask_invalid) {
index <- mask_invalid_values(index, index_type, verbose)
}
# Final validation and reporting
validation_result <- validate_output(index, index_type, verbose)
if (!validation_result) {
# If validation failed, return an NA raster but don't stop the process
if (verbose) message(sprintf("Returning NA raster for failed index: %s", index_type))
na_raster <- red # Use red band as template
terra::values(na_raster) <- NA
names(na_raster) <- index_type
return(na_raster)
}
names(index) <- index_type
return(index)
}
#' Calculate NDVI with time series options
#'
#' @description
#' NDVI calculation specifically designed for time series analysis with
#' date matching, quality filtering, temporal smoothing, and multi-temporal support.
#' **Use this for time series analysis, use calculate_vegetation_index() for single dates.**
#'
#' @param red_data Red band data (files, directory, or raster objects)
#' @param nir_data NIR band data (files, directory, or raster objects)
#' @param clamp_range Range to clamp NDVI values (default: c(-0.2, 1))
#' @param match_by_date Logical: match rasters by date using filenames
#' @param quality_filter Apply quality filtering (remove outliers)
#' @param temporal_smoothing Apply temporal smoothing for time series
#' @param verbose Print progress messages
#' @param date_patterns Custom date patterns for matching
#'
#' @return SpatRaster with NDVI layers (single or multi-layer for time series)
#'
#' @details
#' ## When to Use Enhanced vs Basic NDVI:
#'
#' ### Use `calculate_ndvi_enhanced()` for:
#' - **Time series analysis**: Multiple dates, trend analysis
#' - **Quality control**: Remove outliers, temporal smoothing
#' - **Date matching**: Automatic pairing of red/NIR by date
#' - **Multi-temporal studies**: Seasonal analysis, change detection
#'
#' ### Use `calculate_vegetation_index(index_type="NDVI")` for:
#' - **Single date analysis**: One-time calculation
#' - **Different indices**: Want to calculate EVI, SAVI, etc. too
#' - **Quick calculations**: Simple, fast NDVI
#' - **Mixed workflows**: Part of larger vegetation index analysis
#'
#' @examples
#' \dontrun{
#' # These examples require external data files not included with the package
#' # Time series NDVI with date matching
#' ndvi_series <- calculate_ndvi_enhanced(
#' red_data = "/path/to/red/time_series/",
#' nir_data = "/path/to/nir/time_series/",
#' match_by_date = TRUE,
#' quality_filter = TRUE,
#' temporal_smoothing = TRUE
#' )
#'
#' # Simple NDVI (single date with quality control)
#' ndvi_clean <- calculate_ndvi_enhanced(
#' red_data = red_raster,
#' nir_data = nir_raster,
#' quality_filter = TRUE
#' )
#' }
#'
#' @export
calculate_ndvi_enhanced <- function(red_data, nir_data, clamp_range = c(-0.2, 1),
match_by_date = FALSE, quality_filter = FALSE,
temporal_smoothing = FALSE, verbose = FALSE,
date_patterns = NULL) {
# Input validation
if (missing(red_data) || missing(nir_data)) {
stop("Both red_data and nir_data are required", call. = FALSE)
}
if (is.null(red_data) || is.null(nir_data)) {
stop("red_data and nir_data cannot be NULL", call. = FALSE)
}
if (verbose) message("Starting NDVI calculation for time series analysis...")
# Handle different input types with CRS support
if (inherits(red_data, "SpatRaster") && inherits(nir_data, "SpatRaster")) {
# Direct raster objects provided - use them directly
if (verbose) message("Using provided SpatRaster objects directly")
# CRS and geometry checking
if (!identical(terra::crs(red_data), terra::crs(nir_data))) {
if (verbose) message("CRS mismatch detected. Reprojecting NIR to match RED...")
nir_data <- terra::project(nir_data, terra::crs(red_data))
}
if (!terra::compareGeom(red_data, nir_data, stopOnError = FALSE)) {
if (verbose) message("Geometry mismatch detected. Resampling NIR to match RED...")
nir_data <- terra::resample(nir_data, red_data)
}
# Calculate NDVI directly
ndvi <- (nir_data - red_data) / (nir_data + red_data)
# Handle division by zero
zero_sum_mask <- (red_data + nir_data) == 0
if (any(terra::values(zero_sum_mask, mat = FALSE), na.rm = TRUE)) {
ndvi[zero_sum_mask] <- NA
}
# Apply clamp range
if (!is.null(clamp_range)) {
ndvi <- terra::clamp(ndvi, clamp_range[1], clamp_range[2])
}
# Quality filtering
if (quality_filter) {
ndvi <- apply_quality_filter(ndvi, verbose)
}
names(ndvi) <- "NDVI_enhanced"
if (verbose) {
message("NDVI calculation completed successfully!")
message(sprintf(" - %d layer(s) created", terra::nlyr(ndvi)))
valid_values <- terra::values(ndvi, mat = FALSE)
valid_values <- valid_values[!is.na(valid_values)]
if (length(valid_values) > 0) {
message(sprintf(" - Value range: [%.3f, %.3f]",
min(valid_values), max(valid_values)))
}
}
return(ndvi)
}
# Handle file/directory inputs with error handling
if (match_by_date && (is.character(red_data) && is.character(nir_data))) {
if (verbose) message("Matching red and NIR bands by date...")
tryCatch({
matched_df <- match_rasters_by_date(red_data, nir_data, date_patterns, verbose)
red_files <- matched_df$red
nir_files <- matched_df$nir
layer_dates <- matched_df$date
}, error = function(e) {
if (verbose) message("Date matching failed, using direct file order")
red_files <- red_data
nir_files <- nir_data
layer_dates <- paste0("Layer_", seq_along(if(length(red_data) > 1) red_data else c(red_data)))
})
} else {
red_files <- red_data
nir_files <- nir_data
layer_dates <- paste0("Layer_", seq_along(if(length(red_data) > 1) red_data else c(red_data)))
}
# Load rasters with error handling
tryCatch({
red_rasters <- load_raster_data(red_files, verbose = verbose)
nir_rasters <- load_raster_data(nir_files, verbose = verbose)
}, error = function(e) {
stop(sprintf("Failed to load raster data: %s", e$message), call. = FALSE)
})
# Ensure we have matching numbers of rasters
if (length(red_rasters) != length(nir_rasters)) {
min_length <- min(length(red_rasters), length(nir_rasters))
warning(sprintf("Mismatched number of red (%d) and NIR (%d) rasters. Using first %d pairs.",
length(red_rasters), length(nir_rasters), min_length))
red_rasters <- red_rasters[1:min_length]
nir_rasters <- nir_rasters[1:min_length]
layer_dates <- layer_dates[1:min_length]
}
# Calculate NDVI for each pair with CRS handling
ndvi_list <- vector("list", length(red_rasters))
for (i in seq_along(red_rasters)) {
if (verbose && (i %% 5 == 0 || i == length(red_rasters))) {
message(sprintf("Processing pair %d/%d (%s)...", i, length(red_rasters), layer_dates[i]))
}
tryCatch({
red <- red_rasters[[i]]
nir <- nir_rasters[[i]]
# CRS and geometry compatibility
if (!identical(terra::crs(red), terra::crs(nir))) {
if (verbose && i == 1) message("CRS mismatch detected in time series. Reprojecting...")
nir <- terra::project(nir, terra::crs(red))
}
if (!terra::compareGeom(red, nir, stopOnError = FALSE)) {
nir <- terra::resample(nir, red)
}
# Calculate NDVI
ndvi <- (nir - red) / (nir + red)
# Handle division by zero
zero_sum_mask <- (red + nir) == 0
if (any(terra::values(zero_sum_mask, mat = FALSE), na.rm = TRUE)) {
ndvi[zero_sum_mask] <- NA
}
# Apply clamp range
if (!is.null(clamp_range)) {
ndvi <- terra::clamp(ndvi, clamp_range[1], clamp_range[2])
}
# Quality filtering
if (quality_filter) {
ndvi <- apply_quality_filter(ndvi, verbose = FALSE) # Suppress per-layer messages
}
names(ndvi) <- paste0("NDVI_", gsub("[^A-Za-z0-9_]", "_", layer_dates[i]))
ndvi_list[[i]] <- ndvi
}, error = function(e) {
warning(sprintf("Failed to process pair %d: %s", i, e$message))
# Create an NA raster with same dimensions as the first successful one
if (i > 1 && !is.null(ndvi_list[[1]])) {
na_raster <- ndvi_list[[1]]
terra::values(na_raster) <- NA
names(na_raster) <- paste0("NDVI_", gsub("[^A-Za-z0-9_]", "_", layer_dates[i]))
ndvi_list[[i]] <- na_raster
} else {
ndvi_list[[i]] <- NULL
}
})
}
# Remove NULL entries
ndvi_list <- Filter(Negate(is.null), ndvi_list)
if (length(ndvi_list) == 0) {
stop("No NDVI layers could be calculated successfully", call. = FALSE)
}
# Combine results
result <- if (length(ndvi_list) == 1) {
ndvi_list[[1]]
} else {
tryCatch({
terra::rast(ndvi_list)
}, error = function(e) {
warning("Could not combine layers into stack, returning first layer only")
ndvi_list[[1]]
})
}
# Apply temporal smoothing if requested and we have multiple layers
if (temporal_smoothing && terra::nlyr(result) > 1) {
if (verbose) message("Applying temporal smoothing...")
tryCatch({
result <- apply_temporal_smoothing(result, verbose)
}, error = function(e) {
warning(sprintf("Temporal smoothing failed: %s", e$message))
})
}
if (verbose) {
message("NDVI time series calculation completed successfully!")
message(sprintf(" - %d layer(s) created", terra::nlyr(result)))
if (terra::nlyr(result) > 0) {
valid_values <- terra::values(result, mat = FALSE)
valid_values <- valid_values[!is.na(valid_values)]
if (length(valid_values) > 0) {
message(sprintf(" - Value range: [%.3f, %.3f]",
min(valid_values), max(valid_values)))
}
}
}
return(result)
}
#' Get comprehensive list of available vegetation indices
#'
#' @description
#' Returns detailed information about all 40+ available vegetation indices including
#' formulas, required bands, applications, and references.
#'
#' @param category Filter by category: "all", "basic", "enhanced", "specialized", "stress"
#' @param application Filter by application: "general", "agriculture", "forestry", "stress", "water"
#' @param detailed Return detailed information including formulas and references
#'
#' @return Data frame with vegetation index information
#'
#' @examples
#' \donttest{
#' # All available indices
#' all_indices <- list_vegetation_indices()
#'
#' # Only stress detection indices
#' stress_indices <- list_vegetation_indices(category = "stress")
#'
#' # Detailed information with formulas
#' detailed_info <- list_vegetation_indices(detailed = TRUE)
#'
#' # Agricultural applications only
#' ag_indices <- list_vegetation_indices(application = "agriculture")
#' }
#'
#' @export
list_vegetation_indices <- function(category = "all", application = "all", detailed = FALSE) {
# Comprehensive index database - 40+ indices
indices_db <- data.frame(
Index = c(
# 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",
# Water and Moisture Indices (8)
"NDWI", "MNDWI", "NDMI", "MSI", "NDII", "WI", "SRWI", "LSWI",
# Specialized Applications (10)
"LAI", "FAPAR", "FCOVER", "NBR", "BAI", "NDSI", "GRVI", "VIG", "CI", "GBNDVI"
),
Category = c(
# Basic (10)
rep("basic", 10),
# Enhanced (12)
rep("enhanced", 12),
# Red Edge (10)
rep("advanced", 10),
# Stress (12)
rep("stress", 12),
# Water (8)
rep("water", 8),
# Specialized (10)
rep("specialized", 10)
),
Application = c(
# Basic applications (10)
"general", "soil", "soil", "soil", "general", "general", "biomass", "general", "chlorophyll", "soil",
# Enhanced applications (12)
"atmospheric", "general", "general", "general", "general", "general", "general", "soil", "soil", "soil", "general", "canopy",
# Advanced applications (10)
"stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress",
# Stress applications (12)
"stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress",
# Water applications (8)
"water", "water", "water", "water", "water", "water", "water", "water",
# Specialized applications (10)
"forestry", "forestry", "forestry", "forestry", "forestry", "snow", "general", "general", "canopy", "general"
),
Required_Bands = c(
# Basic indices (10)
"Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR, Blue",
"Red, NIR", "Red, NIR", "Red, NIR", "Green, NIR", "Red, NIR",
# Enhanced indices (12)
"Red, NIR, Blue", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR",
"Red, NIR", "Red, Green, Blue", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR",
# Advanced indices (10)
"NIR, RedEdge", "RedEdge, NIR", "RedEdge, NIR", "RedEdge", "RedEdge, NIR",
"Red, Green", "RedEdge, Green", "RedEdge, Green", "RedEdge, NIR", "Red, Green",
# Stress indices (12)
"Green, NIR", "Red, NIR", "RedEdge, Green", "NIR, SWIR1", "Red, Green",
"Red, Green", "Red, NIR", "Red, NIR", "Red, NIR", "Red, Blue", "Red, NIR", "Red, Blue",
# Water indices (8)
"Green, NIR", "Green, SWIR1", "NIR, SWIR1", "NIR, SWIR1", "NIR, SWIR1",
"NIR, SWIR1", "NIR, SWIR1", "NIR, SWIR1",
# Specialized indices (10)
"Red, NIR", "Red, NIR", "Red, NIR", "NIR, SWIR2", "Red, NIR",
"Green, SWIR1", "Red, Green", "Green, NIR", "Red, Green", "Green, Blue, NIR"
),
stringsAsFactors = FALSE
)
# Add descriptions
indices_db$Description <- c(
# Basic (10)
"Normalized Difference Vegetation Index", "Soil Adjusted Vegetation Index",
"Modified Soil Adjusted Vegetation Index", "Optimized Soil Adjusted Vegetation Index",
"Enhanced Vegetation Index", "Two-band Enhanced Vegetation Index",
"Difference Vegetation Index", "Ratio Vegetation Index", "Green NDVI", "Weighted Difference Vegetation Index",
# Enhanced (12)
"Atmospherically Resistant Vegetation Index", "Renormalized Difference Vegetation Index",
"Perpendicular Vegetation Index", "Infrared Percentage Vegetation Index", "Transformed NDVI",
"Global Environment Monitoring Index", "Visible Atmospherically Resistant Index", "Transformed Soil Adjusted VI",
"Adjusted Transformed Soil Adjusted VI", "Generalized Soil Adjusted VI", "Modified Triangular VI", "Corrected Transformed VI",
# Advanced (10)
"Normalized Difference Red Edge", "MERIS Terrestrial Chlorophyll Index",
"Inverted Red-Edge Chlorophyll Index", "Sentinel-2 Red-Edge Position",
"Plant Senescence Reflectance Index", "Carotenoid Reflectance Index 1",
"Carotenoid Reflectance Index 2", "Anthocyanin Reflectance Index 1",
"Anthocyanin Reflectance Index 2", "Modified Chlorophyll Absorption Ratio Index",
# Stress (12)
"Photochemical Reflectance Index", "Structure Insensitive Pigment Index",
"Canopy Chlorophyll Index", "Normalized Difference Nitrogen Index",
"Chlorophyll Absorption Ratio Index", "Transformed Chlorophyll Absorption Ratio Index",
"Modified Triangular Vegetation Index 1", "Modified Triangular Vegetation Index 2",
"Triangular Vegetation Index", "Normalized Pigment Chlorophyll Index",
"Ratio Analysis of Reflectance Spectra", "Normalized Phaeophytinization Index",
# Water (8)
"Normalized Difference Water Index", "Modified Normalized Difference Water Index",
"Normalized Difference Moisture Index", "Moisture Stress Index",
"Normalized Difference Infrared Index", "Water Index",
"Simple Ratio Water Index", "Land Surface Water Index",
# Specialized (10)
"Leaf Area Index", "Fraction of Absorbed PAR", "Fraction of Vegetation Cover",
"Normalized Burn Ratio", "Burn Area Index", "Normalized Difference Snow Index",
"Green-Red Vegetation Index", "Vegetation Index Green", "Coloration Index",
"Green-Blue NDVI"
)
# Filter by category
if (category != "all") {
indices_db <- indices_db[indices_db$Category == category, ]
}
# Filter by application
if (application != "all") {
indices_db <- indices_db[indices_db$Application == application, ]
}
# Add detailed information if requested
if (detailed) {
indices_db$Formula <- get_index_formulas(indices_db$Index)
indices_db$Range <- get_index_ranges(indices_db$Index)
indices_db$Reference <- get_index_references(indices_db$Index)
}
return(indices_db)
}
#' Calculate multiple vegetation indices at once
#'
#' @description
#' Calculate multiple vegetation indices from the same spectral data in a single operation.
#' Efficient for comparative analysis and comprehensive vegetation assessment.
#' Supports directory input and automatic CRS handling.
#'
#' @param spectral_data Multi-band raster, directory path, or individual bands
#' @param indices Vector of index names to calculate
#' @param output_stack Return as single multi-layer raster (TRUE) or list (FALSE)
#' @param region_boundary Optional region boundary for clipping
#' @param parallel Use parallel processing for multiple indices
#' @param verbose Print progress messages
#' @param ... Additional arguments passed to calculate_vegetation_index
#'
#' @return SpatRaster stack or list of indices
#'
#' @examples
#' \dontrun{
#' # These examples require satellite imagery files (Landsat/Sentinel data etc.)
#' # Calculate multiple basic indices from directory
#' multi_indices <- calculate_multiple_indices(
#' spectral_data = "/path/to/sentinel/bands/",
#' indices = c("NDVI", "EVI", "SAVI", "MSAVI"),
#' auto_detect_bands = TRUE
#' )
#'
#' # Comprehensive vegetation analysis from individual files
#' veg_analysis <- calculate_multiple_indices(
#' red = red_band, nir = nir_band, blue = blue_band,
#' indices = c("NDVI", "EVI", "ARVI", "GNDVI", "DVI"),
#' output_stack = TRUE,
#' region_boundary = "Iowa"
#' )
#'
#' # Directory with custom band matching
#' stress_indices <- calculate_multiple_indices(
#' spectral_data = "/path/to/bands/",
#' indices = c("PRI", "SIPI", "NDRE"),
#' band_names = c("red", "green", "nir", "red_edge"),
#' output_stack = TRUE
#' )
#' }
#'
#' @export
calculate_multiple_indices <- function(spectral_data = NULL, indices = c("NDVI", "EVI", "SAVI"),
output_stack = TRUE, region_boundary = NULL,
parallel = FALSE, verbose = FALSE, ...) {
if (verbose) message(sprintf("Calculating %d vegetation indices...", length(indices)))
# Validate indices
available_indices <- get_available_indices()$Index
invalid_indices <- setdiff(indices, available_indices)
if (length(invalid_indices) > 0) {
stop(sprintf("Invalid indices: %s", paste(invalid_indices, collapse = ", ")))
}
# Calculate each index with error handling
index_results <- list()
if (parallel && requireNamespace("parallel", quietly = TRUE)) {
if (verbose) message("Using parallel processing...")
# Setup parallel processing
n_cores <- min(parallel::detectCores() - 1, length(indices))
if (verbose) message(sprintf("Using %d cores", n_cores))
index_results <- parallel::mclapply(indices, function(idx) {
tryCatch({
calculate_vegetation_index(spectral_data = spectral_data,
index_type = idx, verbose = FALSE, ...)
}, error = function(e) {
warning(sprintf("Failed to calculate %s: %s", idx, e$message))
return(NULL)
})
}, mc.cores = n_cores)
# Remove NULL results
failed_indices <- names(index_results)[sapply(index_results, is.null)]
index_results <- index_results[!sapply(index_results, is.null)]
if (verbose && length(failed_indices) > 0) {
message(sprintf("Failed indices (parallel): %s", paste(failed_indices, collapse = ", ")))
}
} else {
# Sequential processing with better error handling
failed_indices <- c()
for (i in seq_along(indices)) {
idx <- indices[i]
if (verbose) message(sprintf("Calculating %s (%d/%d)...", idx, i, length(indices)))
tryCatch({
result <- calculate_vegetation_index(
spectral_data = spectral_data,
index_type = idx, verbose = FALSE, ...
)
# Check if result is valid (not all NA)
if (all(is.na(terra::values(result, mat = FALSE)))) {
warning(sprintf("Index %s produced only NA values - skipping", idx))
failed_indices <- c(failed_indices, idx)
} else {
index_results[[idx]] <- result
}
}, error = function(e) {
warning(sprintf("Failed to calculate %s: %s", idx, e$message))
failed_indices <- c(failed_indices, idx)
})
}
if (verbose && length(failed_indices) > 0) {
message(sprintf("Failed/skipped indices: %s", paste(failed_indices, collapse = ", ")))
message(sprintf("Successfully calculated: %s", paste(names(index_results), collapse = ", ")))
}
}
# Ensure we have at least some results
if (length(index_results) == 0) {
stop("No vegetation indices could be calculated successfully", call. = FALSE)
}
names(index_results) <- names(index_results) %||% indices[1:length(index_results)]
# Get boundary
if (!is.null(region_boundary)) {
if (verbose) message("Applying region boundary...")
boundary <- get_region_boundary(region_boundary)
boundary_vect <- terra::vect(boundary)
index_results <- lapply(index_results, function(idx_raster) {
if (!is.null(idx_raster) && inherits(idx_raster, "SpatRaster")) {
cropped <- terra::crop(idx_raster, boundary_vect)
terra::mask(cropped, boundary_vect)
} else {
idx_raster
}
})
}
# Return as stack or list - FIXED VERSION
if (output_stack && length(index_results) > 0) {
if (verbose) message("Creating multi-layer stack...")
# Filter out NULL or invalid results
valid_results <- index_results[!sapply(index_results, is.null)]
valid_results <- valid_results[sapply(valid_results, function(x) inherits(x, "SpatRaster"))]
if (length(valid_results) == 0) {
stop("No valid indices to create stack", call. = FALSE)
}
tryCatch({
# Create stack from valid results
result <- terra::rast(valid_results)
# Set names properly - use the names from the list
names(result) <- names(valid_results)
if (verbose) {
message(sprintf("Successfully created stack with %d indices: %s",
terra::nlyr(result), paste(names(result), collapse = ", ")))
}
return(result)
}, error = function(e) {
warning(sprintf("Failed to create raster stack: %s. Returning list instead.", e$message))
return(valid_results)
})
} else {
# Return as list
valid_results <- index_results[!sapply(index_results, is.null)]
if (verbose && length(valid_results) > 0) {
message(sprintf("Returning list with %d indices: %s",
length(valid_results), paste(names(valid_results), collapse = ", ")))
}
return(valid_results)
}
}
#' Specialized crop vegetation analysis
#'
#' @description
#' Perform comprehensive vegetation analysis specifically designed for crop monitoring
#' including growth stage detection, stress identification, and yield prediction support.
#' Handles test scenarios properly with better input validation.
#'
#' @param spectral_data Multi-band spectral data (file, directory, or SpatRaster)
#' @param crop_type Crop type for specialized analysis ("corn", "soybeans", "wheat", "general")
#' @param growth_stage Growth stage if known ("early", "mid", "late", "harvest")
#' @param analysis_type Type of analysis: "comprehensive", "stress", "growth", "yield"
#' @param cdl_mask Optional CDL mask for crop-specific analysis
#' @param reference_data Optional reference data for validation
#' @param output_folder Optional output folder for results
#' @param verbose Print detailed progress
#'
#' @return List with comprehensive vegetation analysis results:
#' \itemize{
#' \item \code{vegetation_indices}: SpatRaster with calculated indices
#' \item \code{analysis_results}: Detailed analysis results by type
#' \item \code{metadata}: Analysis metadata and parameters
#' }
#'
#' @details
#' ## Crop-Specific Index Selection:
#' - **Corn**: NDVI, EVI, GNDVI, DVI, RVI, PRI
#' - **Soybeans**: NDVI, EVI, SAVI, GNDVI, PRI
#' - **Wheat**: NDVI, EVI, SAVI, DVI
#' - **General**: NDVI, EVI, SAVI, GNDVI, DVI, RVI
#'
#' ## Analysis Types:
#' - **comprehensive**: All analyses (stress, growth, yield)
#' - **stress**: Focus on stress detection indices
#' - **growth**: Growth stage analysis
#' - **yield**: Yield prediction support
#'
#' ## Output Structure:
#'
#' The function returns a list with three main components:
#'
#' ### 1. vegetation_indices (SpatRaster)
#' Multi-layer raster with calculated indices (NDVI, EVI, etc.)
#'
#' ### 2. analysis_results (List)
#'
#' **stress_analysis** (if requested):
#' - Percentage of pixels in each stress category
#' - Categories: healthy (NDVI 0.6-1.0), moderate stress (0.4-0.6), severe stress (0.0-0.4)
#' - Includes mean, median, std_dev, and thresholds used
#'
#' **growth_analysis** (if requested):
#' - Predicted growth stage based on NDVI patterns
#' - Stage confidence (0-1 scale)
#' - Detailed statistics for each index
#' - Growth stages: emergence, vegetative, reproductive, maturity (crop-specific)
#'
#' **yield_analysis** (if requested):
#' - **Composite Yield Index**: Normalized 0-1 score combining multiple indices
#' - 0.0 = Very low yield potential
#' - 0.5 = Medium yield potential
#' - 1.0 = Maximum yield potential
#' - **Yield Potential Class**: Categorical (Low, Medium, High, Very High)
#' - **Index Contributions**: How each index contributed to composite score
#' - Calculation: Each index (NDVI, EVI, GNDVI, DVI, RVI) is normalized to 0-1,
#' then averaged to create composite score
#'
#' **summary_statistics**:
#' - Basic stats (mean, std, min, max, percentiles) for all indices
#' - Coverage percentage and pixel counts
#'
#' ### 3. metadata (List)
#' Processing information: crop_type, indices_used, processing_date, spatial properties
#'
#' ## Example Interpretation:
#'
#' ```r
#' result <- analyze_crop_vegetation(data, crop_type = "corn")
#'
#' # Stress Assessment
#' stress <- result$analysis_results$stress_analysis$NDVI
#' cat(sprintf("Healthy: %.1f%%, Stressed: %.1f%%\n",
#' stress$healthy_percentage,
#' stress$severe_stress_percentage))
#'
#' # Growth Stage
#' stage <- result$analysis_results$growth_analysis$predicted_growth_stage
#' cat(sprintf("Growth stage: %s\n", stage))
#'
#' # Yield Potential
#' yield <- result$analysis_results$yield_analysis
#' cat(sprintf("Yield potential: %s (score: %.2f)\n",
#' yield$yield_potential_class,
#' yield$composite_yield_index))
#' ```
#'
#' ## Important Notes:
#'
#' - **Composite Yield Index** is a vegetation-based proxy, not a direct yield prediction
#' - Thresholds are based on literature and may need regional calibration
#' - Results should be validated with ground truth data
#' - For detailed output documentation, see package vignette
#'
#' ## Analysis Types:
#' - **comprehensive**: All analyses (stress, growth, yield)
#' - **stress**: Focus on stress detection indices
#' - **growth**: Growth stage analysis
#' - **yield**: Yield prediction support
#'
#' @examples
#' \dontrun{
#' # These examples require actual spectral data
#' # Comprehensive corn analysis
#' corn_analysis <- analyze_crop_vegetation(
#' spectral_data = sentinel_data,
#' crop_type = "corn",
#' analysis_type = "comprehensive",
#' cdl_mask = corn_mask
#' )
#'
#' # Access results
#' corn_analysis$vegetation_indices # SpatRaster with indices
#' corn_analysis$analysis_results$stress_analysis # Stress detection results
#' corn_analysis$metadata$indices_used # Which indices were calculated
#'
#' # Stress detection in soybeans
#' stress_analysis <- analyze_crop_vegetation(
#' spectral_data = landsat_stack,
#' crop_type = "soybeans",
#' analysis_type = "stress",
#' growth_stage = "mid"
#' )
#' }
#'
#' \donttest{
#' # Example with mock spectral data
#' # Create mock multi-band raster (simulating satellite data)
#' red_band <- terra::rast(nrows = 5, ncols = 5, crs = "EPSG:4326")
#' nir_band <- terra::rast(nrows = 5, ncols = 5, crs = "EPSG:4326")
#' terra::values(red_band) <- runif(25, 0.1, 0.3) # Typical red values
#' terra::values(nir_band) <- runif(25, 0.4, 0.8) # Typical NIR values
#' spectral_stack <- c(red_band, nir_band)
#' names(spectral_stack) <- c("red", "nir")
#'
#' # Analyze with mock data
#' result <- analyze_crop_vegetation(spectral_stack, crop_type = "general")
#' print(names(result)) # Should show analysis components
#' }
#'
#' @export
analyze_crop_vegetation <- function(spectral_data, crop_type = "general",
growth_stage = "unknown", analysis_type = "comprehensive",
cdl_mask = NULL, reference_data = NULL,
output_folder = NULL, verbose = FALSE) {
if (verbose) message(sprintf("Starting specialized %s vegetation analysis for %s...",
analysis_type, crop_type))
# Input validation and loading
if (is.character(spectral_data)) {
if (length(spectral_data) == 1) {
if (dir.exists(spectral_data)) {
if (verbose) message(sprintf("Loading spectral data from directory: %s", spectral_data))
spectral_data <- load_bands_from_directory(spectral_data, verbose = verbose)
} else if (file.exists(spectral_data)) {
if (verbose) message(sprintf("Loading spectral data from file: %s", basename(spectral_data)))
spectral_data <- terra::rast(spectral_data)
} else {
stop(sprintf("Spectral data path does not exist: %s", spectral_data), call. = FALSE)
}
} else {
if (verbose) message("Loading and stacking multiple spectral files")
spectral_data <- load_and_stack_bands(spectral_data, verbose = verbose)
}
}
# Validate spectral data
if (!inherits(spectral_data, "SpatRaster")) {
stop("spectral_data must be a SpatRaster object, file path, or directory", call. = FALSE)
}
if (terra::nlyr(spectral_data) < 2) {
stop("spectral_data must have at least 2 bands for vegetation analysis", call. = FALSE)
}
# Select indices based on crop type and analysis
selected_indices <- select_crop_indices(crop_type, analysis_type, growth_stage)
if (verbose) {
message(sprintf("Selected %d indices for analysis: %s",
length(selected_indices), paste(selected_indices, collapse = ", ")))
}
# Apply CDL mask if provided
if (!is.null(cdl_mask)) {
if (verbose) message("Applying crop mask...")
spectral_data <- terra::mask(spectral_data, cdl_mask)
}
# Calculate vegetation indices using the main function
vegetation_indices <- calculate_multiple_indices(
spectral_data = spectral_data,
indices = selected_indices,
output_stack = TRUE,
auto_detect_bands = TRUE,
verbose = verbose
)
# Perform analysis based on type
analysis_results <- list()
if (analysis_type %in% c("comprehensive", "stress")) {
# Stress detection analysis
if (verbose) message("Performing stress detection analysis...")
analysis_results$stress_analysis <- detect_vegetation_stress(
vegetation_indices, crop_type, selected_indices, verbose
)
}
if (analysis_type %in% c("comprehensive", "growth")) {
# Growth stage analysis
if (verbose) message("Performing growth stage analysis...")
analysis_results$growth_analysis <- analyze_growth_stage(
vegetation_indices, crop_type, growth_stage, verbose
)
}
if (analysis_type %in% c("comprehensive", "yield")) {
# Yield prediction support
if (verbose) message("Performing yield prediction analysis...")
analysis_results$yield_analysis <- analyze_yield_potential(
vegetation_indices, crop_type, selected_indices, verbose
)
}
# Generate summary statistics
analysis_results$summary_statistics <- calculate_vegetation_statistics(
vegetation_indices, selected_indices, verbose
)
# Validation if reference data provided
if (!is.null(reference_data)) {
if (verbose) message("Performing validation analysis...")
analysis_results$validation <- validate_vegetation_analysis(
vegetation_indices, reference_data, verbose
)
}
# Save results if output folder specified
if (!is.null(output_folder)) {
save_vegetation_analysis_results(analysis_results, vegetation_indices,
output_folder, crop_type, verbose)
}
# Combine all results
final_results <- list(
vegetation_indices = vegetation_indices,
analysis_results = analysis_results,
metadata = list(
crop_type = crop_type,
growth_stage = growth_stage,
analysis_type = analysis_type,
indices_used = selected_indices,
processing_date = Sys.time(),
input_bands = terra::nlyr(spectral_data),
spatial_resolution = terra::res(spectral_data),
spatial_extent = as.vector(terra::ext(spectral_data))
)
)
if (verbose) message("Crop vegetation analysis completed successfully!")
return(final_results)
}
# ==================== HELPER FUNCTIONS ==================== #
#' Load bands from directory
#' @keywords internal
load_bands_from_directory <- function(directory, band_names = NULL, verbose = FALSE) {
# Find all raster files
files <- list.files(directory, pattern = "\\.(tif|tiff|img)$",
full.names = TRUE, ignore.case = TRUE)
if (length(files) == 0) {
stop(sprintf("No raster files found in directory: %s", directory), call. = FALSE)
}
if (verbose) message(sprintf("Found %d raster files in directory", length(files)))
# Try to auto-detect which files are which bands
if (is.null(band_names)) {
# Use standard patterns to identify bands
band_patterns <- list(
red = "red|RED|Red|B04|b04|B4|b4",
nir = "nir|NIR|Nir|B08|b08|B8|b8|B05|b05|B5|b5",
blue = "blue|BLUE|Blue|B02|b02|B2|b2",
green = "green|GREEN|Green|B03|b03|B3|b3",
swir1 = "swir1|SWIR1|B11|b11|B6|b6",
swir2 = "swir2|SWIR2|B12|b12|B7|b7",
red_edge = "rededge|RedEdge|red_edge|B05|b05"
)
band_files <- list()
for (band_type in names(band_patterns)) {
pattern <- band_patterns[[band_type]]
matches <- files[grepl(pattern, basename(files), ignore.case = TRUE)]
if (length(matches) > 0) {
band_files[[band_type]] <- matches[1] # Take first match
if (verbose) message(sprintf(" Detected %s band: %s", band_type, basename(matches[1])))
}
}
if (length(band_files) == 0) {
stop("Could not auto-detect any bands. Please specify band_names parameter.", call. = FALSE)
}
# Stack the detected bands
band_rasters <- lapply(band_files, terra::rast)
stacked <- terra::rast(band_rasters)
names(stacked) <- names(band_files)
} else {
# Use provided band names to match files
if (verbose) message("Using provided band names for file matching")
band_files <- list()
for (i in seq_along(band_names)) {
band_name <- band_names[i]
# Look for files containing the band name
matches <- files[grepl(band_name, basename(files), ignore.case = TRUE)]
if (length(matches) > 0) {
band_files[[band_name]] <- matches[1]
if (verbose) message(sprintf(" Matched %s: %s", band_name, basename(matches[1])))
} else {
warning(sprintf("No file found for band: %s", band_name))
}
}
if (length(band_files) == 0) {
stop("Could not match any files to provided band names.", call. = FALSE)
}
# Stack the matched bands
band_rasters <- lapply(band_files, terra::rast)
stacked <- terra::rast(band_rasters)
names(stacked) <- names(band_files)
}
return(stacked)
}
#' Load and stack individual band files
#' @keywords internal
load_and_stack_bands <- function(file_list, band_names = NULL, verbose = FALSE) {
# Validate files exist
missing_files <- file_list[!file.exists(file_list)]
if (length(missing_files) > 0) {
stop(sprintf("Files not found: %s", paste(missing_files, collapse = ", ")), call. = FALSE)
}
# Load all rasters
rasters <- lapply(file_list, terra::rast)
# Check compatibility and align if needed
reference_raster <- rasters[[1]]
reference_crs <- terra::crs(reference_raster)
for (i in 2:length(rasters)) {
# Fix CRS if needed
if (!identical(terra::crs(rasters[[i]]), reference_crs)) {
if (verbose) message(sprintf("Reprojecting band %d to match reference CRS", i))
rasters[[i]] <- terra::project(rasters[[i]], reference_crs)
}
# Align geometry if needed
if (!terra::compareGeom(reference_raster, rasters[[i]], stopOnError = FALSE)) {
if (verbose) message(sprintf("Resampling band %d to match reference geometry", i))
rasters[[i]] <- terra::resample(rasters[[i]], reference_raster)
}
}
# Stack the rasters
stacked <- terra::rast(rasters)
# Name the layers
if (!is.null(band_names) && length(band_names) == length(file_list)) {
names(stacked) <- band_names
} else {
# Use filenames without extension
names(stacked) <- tools::file_path_sans_ext(basename(file_list))
}
if (verbose) message(sprintf("Successfully stacked %d bands", terra::nlyr(stacked)))
return(stacked)
}
#' Get available indices
#' @keywords internal
get_available_indices <- function() {
list_vegetation_indices(detailed = FALSE)
}
#' Extract bands from multi-band raster
#' @keywords internal
extract_bands_from_raster <- function(spectral_data, auto_detect, band_names, verbose) {
if (verbose) message("Extracting bands from multi-band raster...")
n_bands <- terra::nlyr(spectral_data)
layer_names <- names(spectral_data)
bands <- list(red = NULL, nir = NULL, blue = NULL, green = NULL,
swir1 = NULL, swir2 = NULL, red_edge = NULL, coastal = NULL, nir2 = NULL)
if (auto_detect || is.null(band_names)) {
# Auto-detect bands based on common naming conventions
bands <- auto_detect_spectral_bands(spectral_data, layer_names, verbose)
} else {
# Use provided band names
bands <- map_custom_band_names(spectral_data, band_names, verbose)
}
# FALLBACK: If auto-detection failed, use direct name matching
if (is.null(bands$red) && "red" %in% layer_names) {
bands$red <- spectral_data[["red"]]
if (verbose) message("Direct match found for red band")
}
if (is.null(bands$green) && "green" %in% layer_names) {
bands$green <- spectral_data[["green"]]
if (verbose) message("Direct match found for green band")
}
if (is.null(bands$blue) && "blue" %in% layer_names) {
bands$blue <- spectral_data[["blue"]]
if (verbose) message("Direct match found for blue band")
}
if (is.null(bands$nir) && "nir" %in% layer_names) {
bands$nir <- spectral_data[["nir"]]
if (verbose) message("Direct match found for nir band")
}
if (is.null(bands$red_edge) && "red_edge" %in% layer_names) {
bands$red_edge <- spectral_data[["red_edge"]]
if (verbose) message("Direct match found for red_edge band")
}
if (is.null(bands$swir1) && "swir1" %in% layer_names) {
bands$swir1 <- spectral_data[["swir1"]]
if (verbose) message("Direct match found for swir1 band")
}
if (is.null(bands$swir2) && "swir2" %in% layer_names) {
bands$swir2 <- spectral_data[["swir2"]]
if (verbose) message("Direct match found for swir2 band")
}
if (is.null(bands$coastal) && "coastal" %in% layer_names) {
bands$coastal <- spectral_data[["coastal"]]
if (verbose) message("Direct match found for coastal band")
}
# Final fallback: use band positions if names don't work
if (is.null(bands$red) && n_bands >= 1) {
bands$red <- spectral_data[[1]]
if (verbose) message("Using band 1 as red (fallback)")
}
if (is.null(bands$nir) && n_bands >= 4) {
bands$nir <- spectral_data[[4]]
if (verbose) message("Using band 4 as nir (fallback)")
}
if (is.null(bands$green) && n_bands >= 2) {
bands$green <- spectral_data[[2]]
if (verbose) message("Using band 2 as green (fallback)")
}
if (is.null(bands$blue) && n_bands >= 3) {
bands$blue <- spectral_data[[3]]
if (verbose) message("Using band 3 as blue (fallback)")
}
return(bands)
}
#' Auto-detect spectral bands
#' @keywords internal
auto_detect_spectral_bands <- function(spectral_data, layer_names, verbose) {
bands <- list(red = NULL, nir = NULL, blue = NULL, green = NULL,
swir1 = NULL, swir2 = NULL, red_edge = NULL, coastal = NULL, nir2 = NULL)
# Comprehensive band naming patterns
patterns <- list(
red = c("red", "RED", "Red", "B4", "b4", "band4", "Band_4", "Band4", "B04", "b04"),
nir = c("nir", "NIR", "Nir", "B8", "b8", "band8", "Band_8", "Band8", "B5", "b5", "near_infrared", "B08", "b08"),
blue = c("blue", "BLUE", "Blue", "B2", "b2", "band2", "Band_2", "Band2", "B02", "b02"),
green = c("green", "GREEN", "Green", "B3", "b3", "band3", "Band_3", "Band3", "B03", "b03"),
swir1 = c("swir1", "SWIR1", "Swir1", "B11", "b11", "B6", "b6", "shortwave_infrared_1", "SWIR_1"),
swir2 = c("swir2", "SWIR2", "Swir2", "B12", "b12", "B7", "b7", "shortwave_infrared_2", "SWIR_2"),
red_edge = c("rededge", "RedEdge", "red_edge", "RED_EDGE", "B5", "b5", "RE", "re", "B05", "b05"),
coastal = c("coastal", "COASTAL", "Coastal", "B1", "b1", "aerosol", "B01", "b01", "COASTAL_AEROSOL")
)
# Try exact matches first
for (band_type in names(patterns)) {
for (pattern in patterns[[band_type]]) {
if (pattern %in% layer_names) {
bands[[band_type]] <- spectral_data[[pattern]]
if (verbose) message(sprintf("Exact match detected %s band: %s", band_type, pattern))
break
}
}
}
# Try partial matches if exact matches failed
for (band_type in names(patterns)) {
if (is.null(bands[[band_type]])) {
for (pattern in patterns[[band_type]]) {
matches <- which(grepl(pattern, layer_names, ignore.case = TRUE))
if (length(matches) > 0) {
bands[[band_type]] <- spectral_data[[matches[1]]]
if (verbose) message(sprintf("Partial match detected %s band: %s", band_type, layer_names[matches[1]]))
break
}
}
}
}
return(bands)
}
#' Load and validate spectral band
#' @keywords internal
load_and_validate_band <- function(band_data, band_name, required = FALSE) {
if (is.null(band_data)) {
if (required) {
stop(sprintf("%s band is required but not provided", band_name), call. = FALSE)
}
return(NULL)
}
if (is.character(band_data)) {
band_raster <- terra::rast(band_data)
} else {
band_raster <- band_data
}
# Basic validation
if (terra::ncell(band_raster) == 0) {
stop(sprintf("%s band has no cells", band_name), call. = FALSE)
}
return(band_raster)
}
#' Validate required bands for specific indices
#' @keywords internal
validate_required_bands <- function(index_type, blue, green, swir1, swir2, red_edge, coastal, nir2, verbose = FALSE) {
requirements <- list(
"EVI" = list(blue = blue),
"GNDVI" = list(green = green),
"ARVI" = list(blue = blue),
"VARI" = list(blue = blue, green = green),
"NDWI" = list(green = green),
"MNDWI" = list(green = green, swir1 = swir1),
"NDMI" = list(swir1 = swir1),
"MSI" = list(swir1 = swir1),
"NDII" = list(swir1 = swir1),
"NDRE" = list(red_edge = red_edge),
"MTCI" = list(red_edge = red_edge),
"IRECI" = list(red_edge = red_edge),
"S2REP" = list(red_edge = red_edge),
"PSRI" = list(red_edge = red_edge),
"NBR" = list(swir2 = swir2),
"BAI" = list(swir1 = swir1),
"NDSI" = list(green = green, swir1 = swir1),
"CCI" = list(red_edge = red_edge),
"LAI" = list(blue = blue),
"FAPAR" = list(blue = blue),
"MTVI1" = list(green = green),
"MTVI2" = list(green = green),
"WI" = list(swir1 = swir1),
"SRWI" = list(swir1 = swir1),
"LSWI" = list(swir1 = swir1),
"CI" = list(green = green),
"GBNDVI" = list(green = green, blue = blue),
"NPCI" = list(blue = blue),
"NPQI" = list(blue = blue)
)
if (index_type %in% names(requirements)) {
required_bands <- requirements[[index_type]]
for (band_name in names(required_bands)) {
band_data <- required_bands[[band_name]]
# Check if band is NULL or invalid
if (is.null(band_data)) {
if (verbose) message(sprintf("%s index requires %s band but it's NULL", index_type, band_name))
return(FALSE)
}
# Check if it's a valid SpatRaster
if (!inherits(band_data, "SpatRaster")) {
if (verbose) message(sprintf("%s index requires %s band but it's not a SpatRaster", index_type, band_name))
return(FALSE)
}
# Check if it has valid data
if (terra::ncell(band_data) == 0) {
if (verbose) message(sprintf("%s index requires %s band but it has no cells", index_type, band_name))
return(FALSE)
}
}
}
return(TRUE)
}
#' Calculate index by type
#' @keywords internal
calculate_index_by_type <- function(index_type, red, nir, blue, green, swir1, swir2,
red_edge, coastal, nir2, verbose) {
index <- switch(index_type,
# ==================== BASIC VEGETATION INDICES ====================
"NDVI" = (nir - red) / (nir + red),
"SAVI" = {
L <- 0.5 # Standard soil brightness correction factor
((nir - red) / (nir + red + L)) * (1 + L)
},
"MSAVI" = {
# Modified SAVI - self-adjusting L parameter
0.5 * (2 * nir + 1 - sqrt((2 * nir + 1)^2 - 8 * (nir - red)))
},
"OSAVI" = {
# Optimized SAVI with L = 0.16
(nir - red) / (nir + red + 0.16)
},
"EVI" = {
# Enhanced Vegetation Index
2.5 * ((nir - red) / (nir + 6 * red - 7.5 * blue + 1))
},
"EVI2" = {
# Two-band Enhanced Vegetation Index
2.5 * ((nir - red) / (nir + 2.4 * red + 1))
},
"DVI" = nir - red,
"RVI" = nir / red,
"GNDVI" = (nir - green) / (nir + green),
"WDVI" = {
# Weighted Difference Vegetation Index
# Using a typical soil line slope of 0.5
nir - 0.5 * red
},
# ==================== ENHANCED VEGETATION INDICES ====================
"ARVI" = {
# Atmospherically Resistant Vegetation Index
rb <- red - 1 * (red - blue) # Self-correction for atmospheric scattering
(nir - rb) / (nir + rb)
},
"RDVI" = {
# Renormalized Difference Vegetation Index
(nir - red) / sqrt(nir + red)
},
"PVI" = {
# Perpendicular Vegetation Index (using standard coefficients)
# Based on soil line: NIR = a*Red + b, with typical values a=1.5, b=10
a <- 1.5; b <- 10
(nir - a * red - b) / sqrt(1 + a^2)
},
"IPVI" = {
# Infrared Percentage Vegetation Index
nir / (nir + red)
},
"TNDVI" = {
# Transformed NDVI
sqrt(((nir - red) / (nir + red)) + 0.5)
},
"GEMI" = {
# Global Environment Monitoring Index
eta <- (2 * (nir^2 - red^2) + 1.5 * nir + 0.5 * red) / (nir + red + 0.5)
eta * (1 - 0.25 * eta) - (red - 0.125) / (1 - red)
},
"VARI" = {
# Visible Atmospherically Resistant Index
(green - red) / (green + red - blue)
},
"TSAVI" = {
# Transformed Soil Adjusted Vegetation Index
# Using standard parameters: a=1.5, b=10, X=0.08
a <- 1.5; b <- 10; X <- 0.08
(a * (nir - a * red - b)) / (red + a * nir - a * b + X * (1 + a^2))
},
"ATSAVI" = {
# Adjusted Transformed Soil Adjusted Vegetation Index
a <- 1.22; b <- 0.03; X <- 0.08
(a * (nir - a * red - b)) / (a * nir + red - a * b + X * (1 + a^2))
},
"GESAVI" = {
# Generalized Soil Adjusted Vegetation Index
Z <- 0.5 # Adjustment factor
((nir - red) / (nir + red + Z)) * (1 + Z)
},
"MTVI" = {
# Modified Triangular Vegetation Index (without red edge)
1.2 * (1.2 * (nir - green) - 2.5 * (red - green))
},
"CTVI" = {
# Corrected Transformed Vegetation Index
ndvi_val <- (nir - red) / (nir + red)
((ndvi_val + 0.5) / abs(ndvi_val + 0.5)) * sqrt(abs(ndvi_val + 0.5))
},
# ==================== RED EDGE INDICES ====================
"NDRE" = (nir - red_edge) / (nir + red_edge),
"MTCI" = {
# MERIS Terrestrial Chlorophyll Index
(red_edge - red) / (nir - red)
},
"IRECI" = {
# Inverted Red-Edge Chlorophyll Index
(red_edge - red) / (red_edge / nir)
},
"S2REP" = {
# Sentinel-2 Red-Edge Position
705 + 35 * ((red + red_edge) / 2 - red) / (red_edge - red)
},
"PSRI" = {
# Plant Senescence Reflectance Index
(red - green) / red_edge
},
"CRI1" = (1 / green) - (1 / red),
"CRI2" = (1 / green) - (1 / red_edge),
"ARI1" = (1 / green) - (1 / red_edge),
"ARI2" = nir * ((1 / green) - (1 / red_edge)),
"MCARI" = {
# Modified Chlorophyll Absorption Ratio Index
((red_edge - red) - 0.2 * (red_edge - green)) * (red_edge / red)
},
# ==================== STRESS DETECTION INDICES ====================
"PRI" = {
# Photochemical Reflectance Index
# Using 531nm (green) as proxy for 531nm reference band
(green - nir) / (green + nir)
},
"SIPI" = {
# Structure Insensitive Pigment Index
(nir - red) / (nir - green)
},
"CCI" = {
# Canopy Chlorophyll Index
(red_edge - red) / (red_edge + red)
},
"NDNI" = {
# Normalized Difference Nitrogen Index
(log(1 / nir) - log(1 / swir1)) / (log(1 / nir) + log(1 / swir1))
},
"CARI" = {
# Chlorophyll Absorption Ratio Index
red_edge * (red / green)
},
"TCARI" = {
# Transformed Chlorophyll Absorption Ratio Index
3 * ((red_edge - red) - 0.2 * (red_edge - green) * (red_edge / red))
},
"MTVI1" = {
# Modified Triangular Vegetation Index 1
1.2 * (1.2 * (nir - green) - 2.5 * (red - green))
},
"MTVI2" = {
# Modified Triangular Vegetation Index 2
1.5 * (1.2 * (nir - green) - 2.5 * (red - green)) /
sqrt((2 * nir + 1)^2 - (6 * nir - 5 * sqrt(red)) - 0.5)
},
"TVI" = {
# Triangular Vegetation Index
0.5 * (120 * (nir - green) - 200 * (red - green))
},
"NPCI" = {
# Normalized Pigment Chlorophyll Index
(red - blue) / (red + blue)
},
"RARS" = {
# Ratio Analysis of Reflectance Spectra
red / nir
},
"NPQI" = {
# Normalized Phaeophytinization Index
(red - blue) / (red + blue)
},
# ==================== WATER/MOISTURE INDICES ====================
"NDWI" = {
# Original NDWI (McFeeters 1996) - Green and NIR for water bodies
(green - nir) / (green + nir)
},
"MNDWI" = {
# Modified NDWI (Xu 2006) - Green and SWIR1 for water bodies
(green - swir1) / (green + swir1)
},
"NDMI" = {
# Normalized Difference Moisture Index (Gao 1996) - NIR and SWIR1 for vegetation moisture
(nir - swir1) / (nir + swir1)
},
"MSI" = {
# Moisture Stress Index
swir1 / nir
},
"NDII" = {
# Normalized Difference Infrared Index
(nir - swir1) / (nir + swir1)
},
"WI" = {
# Water Index
nir / swir1
},
"SRWI" = {
# Simple Ratio Water Index
nir / swir1
},
"LSWI" = {
# Land Surface Water Index
(nir - swir1) / (nir + swir1)
},
# ==================== SPECIALIZED INDICES ====================
"LAI" = {
# Leaf Area Index approximation
3.618 * (2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)) - 0.118
},
"FAPAR" = {
# Fraction of Absorbed PAR
-0.161 + 1.257 * ((nir - red) / (nir + red))
},
"FCOVER" = {
# Fractional Vegetation Cover
ndvi_val <- (nir - red) / (nir + red)
-2.274 + 4.336 * ndvi_val - 1.33 * ndvi_val^2
},
"NBR" = {
# Normalized Burn Ratio
(nir - swir2) / (nir + swir2)
},
"BAI" = {
# Burn Area Index
1 / ((0.1 - red)^2 + (0.06 - nir)^2)
},
"NDSI" = {
# Normalized Difference Snow Index
(green - swir1) / (green + swir1)
},
"GRVI" = {
# Green-Red Vegetation Index
(green - red) / (green + red)
},
"VIG" = {
# Vegetation Index Green
(green - red) / (green + red)
},
"CI" = {
# Coloration Index
(red - green) / red
},
"GBNDVI" = {
# Green-Blue NDVI
(nir - (green + blue)) / (nir + green + blue)
},
# If index not found, throw error
stop(sprintf("Unsupported index type: %s", index_type), call. = FALSE)
)
# Handle special cases for division by zero and invalid operations
index <- handle_index_edge_cases(index, index_type, verbose)
return(index)
}
#' Handle edge cases for index calculations
#' @keywords internal
handle_index_edge_cases <- function(index, index_type, verbose) {
# Handle division by zero cases (ratio indices)
if (index_type %in% c("RVI", "MSI", "WI", "SRWI", "RARS")) {
# For ratio indices, handle zero denominators
is_invalid <- is.infinite(terra::values(index, mat = FALSE)) |
is.nan(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to division by zero in %s", n_invalid, index_type))
}
}
}
# Handle square root of negative numbers
if (index_type %in% c("MSAVI", "MTVI2", "RDVI", "TNDVI", "CTVI")) {
is_invalid <- is.nan(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to square root of negative in %s", n_invalid, index_type))
}
}
}
# Handle logarithm edge cases
if (index_type %in% c("NDNI")) {
is_invalid <- is.infinite(terra::values(index, mat = FALSE)) |
is.nan(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to logarithm issues in %s", n_invalid, index_type))
}
}
}
# Handle reciprocal operations (CRI indices)
if (index_type %in% c("CRI1", "CRI2", "ARI1", "ARI2")) {
is_invalid <- is.infinite(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to reciprocal operations in %s", n_invalid, index_type))
}
}
}
# Handle S2REP special case (Red Edge Position)
if (index_type == "S2REP") {
# Constrain to reasonable wavelength range
index[index < 680 | index > 780] <- NA
}
# Handle BAI special case (very large values possible)
if (index_type == "BAI") {
# Cap extremely large values
index[index > 1000] <- 1000
}
return(index)
}
#' Apply quality filter to remove outliers
#' @keywords internal
apply_quality_filter <- function(ndvi, verbose) {
# Calculate statistics
values <- terra::values(ndvi, mat = FALSE)
values <- values[!is.na(values)]
if (length(values) == 0) return(ndvi)
# Use IQR method for outlier detection
q1 <- quantile(values, 0.25, na.rm = TRUE)
q3 <- quantile(values, 0.75, na.rm = TRUE)
iqr <- q3 - q1
# Conservative outlier bounds (1.5 * IQR)
lower_bound <- q1 - 1.5 * iqr
upper_bound <- q3 + 1.5 * iqr
# For NDVI, also apply reasonable physical bounds
lower_bound <- max(lower_bound, -0.5) # NDVI rarely below -0.5 for vegetation
upper_bound <- min(upper_bound, 1.0) # NDVI cannot exceed 1.0
# Mask outliers
outlier_mask <- (ndvi < lower_bound) | (ndvi > upper_bound)
ndvi[outlier_mask] <- NA
if (verbose) {
n_outliers <- sum(terra::values(outlier_mask, mat = FALSE), na.rm = TRUE)
if (n_outliers > 0) {
message(sprintf("Quality filter removed %d outlier pixels (%.1f%%)",
n_outliers, (n_outliers/terra::ncell(ndvi))*100))
}
}
return(ndvi)
}
#' Apply temporal smoothing
#' @keywords internal
apply_temporal_smoothing <- function(ndvi_stack, verbose) {
n_layers <- terra::nlyr(ndvi_stack)
if (n_layers < 3) {
if (verbose) message("Not enough layers for temporal smoothing (need at least 3)")
return(ndvi_stack)
}
if (verbose) message(sprintf("Applying temporal smoothing to %d layers", n_layers))
# Create smoothed layers using moving average
smoothed_layers <- list()
for (i in 1:n_layers) {
if (i == 1) {
# First layer: average with next
smoothed_layers[[i]] <- (ndvi_stack[[i]] + ndvi_stack[[i + 1]]) / 2
} else if (i == n_layers) {
# Last layer: average with previous
smoothed_layers[[i]] <- (ndvi_stack[[i - 1]] + ndvi_stack[[i]]) / 2
} else {
# Middle layers: 3-point moving average
smoothed_layers[[i]] <- (ndvi_stack[[i - 1]] + ndvi_stack[[i]] + ndvi_stack[[i + 1]]) / 3
}
}
# Combine into stack
result <- terra::rast(smoothed_layers)
names(result) <- names(ndvi_stack)
if (verbose) message("Temporal smoothing completed")
return(result)
}
#' Mask invalid values based on index type
#' @keywords internal
mask_invalid_values <- function(index, index_type, verbose) {
# Define reasonable ranges for different indices based on literature
ranges <- list(
# Basic vegetation indices
"NDVI" = c(-1, 1), "SAVI" = c(-1, 1.5), "MSAVI" = c(0, 2), "OSAVI" = c(-1, 1),
"EVI" = c(-1, 3), "EVI2" = c(-1, 3), "DVI" = c(-2, 2), "RVI" = c(0, 30),
"GNDVI" = c(-1, 1), "WDVI" = c(-2, 2),
# Enhanced indices
"ARVI" = c(-1, 1), "RDVI" = c(-2, 2), "PVI" = c(-2, 2), "IPVI" = c(0, 1),
"TNDVI" = c(0, 1.2), "GEMI" = c(-1, 1), "VARI" = c(-1, 1),
"TSAVI" = c(-1, 1), "ATSAVI" = c(-1, 1), "GESAVI" = c(-1, 1.5),
# Red edge indices
"NDRE" = c(-1, 1), "MTCI" = c(0, 8), "IRECI" = c(0, 8), "S2REP" = c(680, 780),
"PSRI" = c(-1, 1), "CRI1" = c(-5, 5), "CRI2" = c(-5, 5),
"ARI1" = c(-5, 5), "ARI2" = c(-50, 50), "MCARI" = c(0, 2),
# Stress indices
"PRI" = c(-1, 1), "SIPI" = c(0, 2), "CCI" = c(-1, 1), "NDNI" = c(-1, 1),
"CARI" = c(0, 10), "TCARI" = c(0, 5), "MTVI1" = c(-2, 2), "MTVI2" = c(0, 8),
"TVI" = c(0, 150), "NPCI" = c(-1, 1), "RARS" = c(0, 5), "NPQI" = c(-1, 1),
# Water/moisture indices
"NDWI" = c(-1, 1), "MNDWI" = c(-1, 1), "NDMI" = c(-1, 1), "MSI" = c(0, 5),
"NDII" = c(-1, 1), "WI" = c(0, 5), "SRWI" = c(0, 5), "LSWI" = c(-1, 1),
# Specialized indices
"LAI" = c(0, 15), "FAPAR" = c(0, 1), "FCOVER" = c(0, 1), "NBR" = c(-1, 1),
"BAI" = c(0, 1000), "NDSI" = c(-1, 1), "GRVI" = c(-1, 1), "VIG" = c(-1, 1),
"CI" = c(-1, 1), "GBNDVI" = c(-1, 1)
)
if (index_type %in% names(ranges)) {
valid_range <- ranges[[index_type]]
invalid_mask <- (index < valid_range[1]) | (index > valid_range[2])
n_invalid <- sum(terra::values(invalid_mask, mat = FALSE), na.rm = TRUE)
if (n_invalid > 0) {
index[invalid_mask] <- NA
if (verbose) {
message(sprintf("Masked %d invalid pixels outside range [%.3f, %.3f] for %s (%.1f%%)",
n_invalid, valid_range[1], valid_range[2], index_type,
(n_invalid/terra::ncell(index))*100))
}
}
} else {
if (verbose) message(sprintf("No predefined range for %s, skipping range validation", index_type))
}
return(index)
}
#' Validate output quality
#' @keywords internal
validate_output <- function(index, index_type, verbose) {
values <- terra::values(index, mat = FALSE)
n_total <- length(values)
n_valid <- sum(!is.na(values))
# Instead of stopping, just warn and continue
if (n_valid == 0) {
warning(sprintf("All %s values are invalid/NA. This index may not be suitable for your data.", index_type),
call. = FALSE)
if (verbose) {
message(sprintf("Warning: %s produced no valid values - continuing with other indices", index_type))
}
return(invisible(FALSE)) # Return FALSE to indicate failure but don't stop
}
# Warn if very low data coverage but continue
coverage_pct <- (n_valid/n_total) * 100
if (coverage_pct < 10) {
warning(sprintf("Only %.1f%% of %s values are valid. This index may not be suitable for your data.",
coverage_pct, index_type), call. = FALSE)
} else if (coverage_pct < 50 && verbose) {
message(sprintf("Warning: Only %.1f%% of %s values are valid", coverage_pct, index_type))
}
if (verbose && n_valid > 0) {
value_range <- range(values, na.rm = TRUE)
message(sprintf("%s calculation completed:", index_type))
message(sprintf(" - Valid pixels: %d/%d (%.1f%%)", n_valid, n_total, coverage_pct))
message(sprintf(" - Value range: [%.6f, %.6f]", value_range[1], value_range[2]))
# Add statistical summary for verbose output
if (n_valid >= 10) {
values_clean <- values[!is.na(values)]
message(sprintf(" - Mean: %.6f, Std Dev: %.6f", mean(values_clean), sd(values_clean)))
message(sprintf(" - Percentiles (25th, 50th, 75th): [%.4f, %.4f, %.4f]",
quantile(values_clean, 0.25), quantile(values_clean, 0.5), quantile(values_clean, 0.75)))
}
}
return(invisible(TRUE)) # Return TRUE to indicate success
}
#' Select appropriate indices for crop analysis
#' @keywords internal
select_crop_indices <- function(crop_type, analysis_type, growth_stage) {
# Base indices that are always reliable with just Red + NIR
base_indices <- c("NDVI", "SAVI", "DVI")
# Crop-specific indices - ONLY USING RED + NIR COMPATIBLE INDICES
crop_indices <- switch(crop_type,
"corn" = c("NDVI", "EVI2", "SAVI", "MSAVI", "DVI", "RVI"),
"soybeans" = c("NDVI", "EVI2", "SAVI", "OSAVI", "DVI"),
"wheat" = c("NDVI", "SAVI", "MSAVI", "DVI"),
"cotton" = c("NDVI", "EVI2", "SAVI", "MSAVI"),
"rice" = c("NDVI", "SAVI", "DVI"),
"general" = c("NDVI", "EVI2", "SAVI", "MSAVI", "OSAVI", "DVI", "RVI"),
base_indices
)
# Analysis-specific additions - ONLY ROBUST RED+NIR INDICES
if (analysis_type %in% c("stress", "comprehensive")) {
# Only add stress indices that work with Red + NIR
# SIPI needs Green band, so we exclude it
# PRI needs different bands, so we exclude it
# Just use NDVI variations for stress detection
stress_indices <- c("RDVI") # Renormalized DVI can indicate stress
crop_indices <- unique(c(crop_indices, stress_indices))
}
if (analysis_type %in% c("yield", "comprehensive")) {
# Only add yield-relevant indices that are stable with Red + NIR
yield_indices <- c("DVI", "RVI")
crop_indices <- unique(c(crop_indices, yield_indices))
}
if (analysis_type %in% c("growth", "comprehensive")) {
# Growth-related indices that work with Red + NIR
growth_indices <- c("DVI", "RVI")
crop_indices <- unique(c(crop_indices, growth_indices))
}
# Remove any indices that definitely need other bands
problematic_indices <- c("SIPI", "PRI", "GNDVI", "EVI", "ARVI", "VARI", "NDWI", "MNDWI",
"NDMI", "MTCI", "PSRI", "CCI", "NDRE", "CARI", "TCARI")
crop_indices <- setdiff(crop_indices, problematic_indices)
# Ensure we always have at least the basic indices
if (length(crop_indices) == 0) {
crop_indices <- base_indices
}
return(unique(crop_indices))
}
#' Detect vegetation stress
#' @keywords internal
detect_vegetation_stress <- function(vegetation_indices, crop_type, indices, verbose) {
stress_results <- list()
# Only analyze stress indices that are actually present
stress_indices <- intersect(indices, c("SIPI", "NDVI", "EVI", "GNDVI"))
for (idx in stress_indices) {
if (idx %in% names(vegetation_indices)) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) == 0) {
stress_results[[idx]] <- list(error = "No valid values for stress analysis")
next
}
# Define stress thresholds based on literature and crop type
thresholds <- switch(idx,
"NDVI" = list(
healthy = c(0.6, 1.0),
moderate_stress = c(0.4, 0.6),
severe_stress = c(0.0, 0.4)
),
"EVI" = list(
healthy = c(0.4, 1.0),
moderate_stress = c(0.2, 0.4),
severe_stress = c(0.0, 0.2)
),
"SIPI" = list(
healthy = c(1.0, 2.0),
moderate_stress = c(0.8, 1.0),
severe_stress = c(0.5, 0.8)
),
"GNDVI" = list(
healthy = c(0.5, 1.0),
moderate_stress = c(0.3, 0.5),
severe_stress = c(0.0, 0.3)
),
# Default thresholds
list(healthy = c(0.5, 1.0), moderate_stress = c(0.3, 0.5), severe_stress = c(0.0, 0.3))
)
# Calculate stress statistics
healthy_pixels <- sum(values >= thresholds$healthy[1] & values <= thresholds$healthy[2])
moderate_stress_pixels <- sum(values >= thresholds$moderate_stress[1] & values <= thresholds$moderate_stress[2])
severe_stress_pixels <- sum(values >= thresholds$severe_stress[1] & values <= thresholds$severe_stress[2])
total_pixels <- length(values)
stress_results[[idx]] <- list(
healthy_percentage = (healthy_pixels / total_pixels) * 100,
moderate_stress_percentage = (moderate_stress_pixels / total_pixels) * 100,
severe_stress_percentage = (severe_stress_pixels / total_pixels) * 100,
mean_value = mean(values),
median_value = median(values),
std_dev = sd(values),
thresholds_used = thresholds,
total_pixels_analyzed = total_pixels
)
if (verbose) {
message(sprintf("Stress analysis for %s:", idx))
message(sprintf(" Healthy: %.1f%%, Moderate stress: %.1f%%, Severe stress: %.1f%%",
stress_results[[idx]]$healthy_percentage,
stress_results[[idx]]$moderate_stress_percentage,
stress_results[[idx]]$severe_stress_percentage))
}
}
}
return(stress_results)
}
#' Analyze growth stage
#' @keywords internal
analyze_growth_stage <- function(vegetation_indices, crop_type, growth_stage, verbose) {
# Extract key indices for growth analysis
growth_indices <- intersect(names(vegetation_indices), c("NDVI", "EVI", "GNDVI", "DVI"))
if (length(growth_indices) == 0) {
return(list(error = "No suitable indices available for growth stage analysis"))
}
growth_analysis <- list()
for (idx in growth_indices) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) == 0) {
growth_analysis[[idx]] <- list(error = "No valid values")
next
}
growth_analysis[[idx]] <- list(
mean = mean(values),
median = median(values),
std_dev = sd(values),
min = min(values),
max = max(values),
range = max(values) - min(values),
percentiles = quantile(values, c(0.1, 0.25, 0.75, 0.9)),
coefficient_of_variation = sd(values) / mean(values),
n_pixels = length(values)
)
}
# Add growth stage classification if NDVI available
if ("NDVI" %in% names(growth_analysis) && !is.null(growth_analysis$NDVI$mean)) {
ndvi_mean <- growth_analysis$NDVI$mean
# Crop-specific growth stage classification
predicted_stage <- switch(crop_type,
"corn" = {
if (ndvi_mean < 0.3) "emergence" else
if (ndvi_mean < 0.6) "vegetative" else
if (ndvi_mean < 0.8) "reproductive" else "maturity"
},
"soybeans" = {
if (ndvi_mean < 0.4) "emergence" else
if (ndvi_mean < 0.65) "vegetative" else
if (ndvi_mean < 0.8) "reproductive" else "maturity"
},
"wheat" = {
if (ndvi_mean < 0.35) "tillering" else
if (ndvi_mean < 0.7) "stem_elongation" else
if (ndvi_mean < 0.8) "grain_filling" else "maturity"
},
# General classification
{
if (ndvi_mean < 0.3) "early" else
if (ndvi_mean < 0.6) "vegetative" else
if (ndvi_mean < 0.8) "reproductive" else "mature"
}
)
# Calculate confidence based on how close to stage boundaries
stage_boundaries <- c(0.3, 0.6, 0.8)
distances_to_boundaries <- abs(ndvi_mean - stage_boundaries)
min_distance <- min(distances_to_boundaries)
stage_confidence <- max(0, 1 - (min_distance / 0.15)) # Confidence decreases as we approach boundaries
growth_analysis$predicted_growth_stage <- predicted_stage
growth_analysis$stage_confidence <- stage_confidence
growth_analysis$crop_type_used = crop_type
if (verbose) {
message(sprintf("Predicted growth stage: %s (confidence: %.2f)", predicted_stage, stage_confidence))
}
}
return(growth_analysis)
}
#' Analyze yield potential
#' @keywords internal
analyze_yield_potential <- function(vegetation_indices, crop_type, indices, verbose) {
# Select yield-relevant indices
yield_indices <- intersect(indices, c("NDVI", "EVI", "GNDVI", "DVI", "RVI"))
if (length(yield_indices) == 0) {
return(list(error = "No suitable indices available for yield analysis"))
}
yield_analysis <- list()
# Calculate composite yield index using available indices
composite_values <- 0
n_indices <- 0
index_contributions <- list()
for (idx in yield_indices) {
if (idx %in% names(vegetation_indices)) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) > 0) {
# Normalize values to 0-1 range for each index
if (idx %in% c("NDVI", "EVI", "GNDVI")) {
# For these indices, higher values generally indicate better yield potential
normalized <- pmax(0, pmin(1, (values - 0.1) / (0.9 - 0.1)))
} else if (idx == "DVI") {
# DVI can be negative, so handle differently
normalized <- pmax(0, pmin(1, (values + 0.5) / (2.0 + 0.5)))
} else if (idx == "RVI") {
# RVI ranges from near 0 to high values
normalized <- pmax(0, pmin(1, (values - 1) / (10 - 1)))
} else {
# Generic normalization
min_val <- quantile(values, 0.05, na.rm = TRUE)
max_val <- quantile(values, 0.95, na.rm = TRUE)
normalized <- pmax(0, pmin(1, (values - min_val) / (max_val - min_val)))
}
mean_normalized <- mean(normalized, na.rm = TRUE)
composite_values <- composite_values + mean_normalized
n_indices <- n_indices + 1
index_contributions[[idx]] <- list(
mean_normalized = mean_normalized,
raw_mean = mean(values, na.rm = TRUE),
raw_std = sd(values, na.rm = TRUE)
)
if (verbose) {
message(sprintf(" %s contribution: %.3f (raw mean: %.3f)", idx, mean_normalized, mean(values, na.rm = TRUE)))
}
}
}
}
if (n_indices > 0) {
composite_index <- composite_values / n_indices
# Crop-specific yield potential classification
yield_potential_class <- switch(crop_type,
"corn" = {
if (composite_index < 0.3) "Low" else
if (composite_index < 0.6) "Medium" else
if (composite_index < 0.8) "High" else "Very High"
},
"soybeans" = {
if (composite_index < 0.35) "Low" else
if (composite_index < 0.65) "Medium" else
if (composite_index < 0.85) "High" else "Very High"
},
"wheat" = {
if (composite_index < 0.3) "Low" else
if (composite_index < 0.6) "Medium" else
if (composite_index < 0.8) "High" else "Very High"
},
# Default classification
{
if (composite_index < 0.3) "Low" else
if (composite_index < 0.6) "Medium" else
if (composite_index < 0.8) "High" else "Very High"
}
)
yield_analysis <- list(
composite_yield_index = composite_index,
yield_potential_class = yield_potential_class,
indices_used = yield_indices,
n_indices_used = n_indices,
index_contributions = index_contributions,
crop_type = crop_type,
classification_confidence = abs(composite_index - 0.5) * 2 # Higher confidence away from middle
)
if (verbose) {
message(sprintf("Composite yield index: %.3f", composite_index))
message(sprintf("Yield potential class: %s", yield_potential_class))
}
} else {
yield_analysis <- list(error = "No valid data for yield analysis")
}
return(yield_analysis)
}
#' Calculate vegetation statistics
#' @keywords internal
calculate_vegetation_statistics <- function(vegetation_indices, indices, verbose) {
if (!inherits(vegetation_indices, "SpatRaster")) {
return(list(error = "Invalid vegetation indices data"))
}
statistics <- list()
for (idx in indices) {
if (idx %in% names(vegetation_indices)) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) > 0) {
statistics[[idx]] <- list(
count = length(values),
mean = mean(values),
median = median(values),
std_dev = sd(values),
min = min(values),
max = max(values),
range = max(values) - min(values),
cv = sd(values) / abs(mean(values)), # Coefficient of variation
percentiles = quantile(values, c(0.05, 0.25, 0.75, 0.95)),
coverage_percent = (length(values) / terra::ncell(vegetation_indices)) * 100
)
# Add histogram data for distribution analysis
if (length(values) >= 100) {
hist_data <- hist(values, breaks = 20, plot = FALSE)
statistics[[idx]]$histogram <- list(
breaks = hist_data$breaks,
counts = hist_data$counts,
density = hist_data$density
)
}
if (verbose) {
stats <- statistics[[idx]]
message(sprintf("%s statistics:", idx))
message(sprintf(" Mean: %.4f +/- %.4f", stats$mean, stats$std_dev))
message(sprintf(" Range: [%.4f, %.4f]", stats$min, stats$max))
message(sprintf(" Coverage: %.1f%% (%d pixels)", stats$coverage_percent, stats$count))
}
} else {
statistics[[idx]] <- list(error = "No valid values for statistics calculation")
if (verbose) message(sprintf("%s: No valid values", idx))
}
} else {
if (verbose) message(sprintf("%s: Index not found in data", idx))
}
}
# Add overall summary
valid_indices <- names(statistics)[!sapply(statistics, function(x) "error" %in% names(x))]
if (length(valid_indices) > 0) {
statistics$summary <- list(
total_indices_calculated = length(valid_indices),
indices_with_valid_data = valid_indices,
total_indices_requested = length(indices),
success_rate = (length(valid_indices) / length(indices)) * 100
)
}
return(statistics)
}
#' Get index formulas
#' @keywords internal
get_index_formulas <- function(indices) {
formulas <- setNames(rep("", length(indices)), indices)
formula_db <- list(
# ==================== BASIC VEGETATION INDICES (10) ====================
"NDVI" = "(NIR - Red) / (NIR + Red)",
"SAVI" = "((NIR - Red) / (NIR + Red + L)) * (1 + L), where L=0.5",
"MSAVI" = "0.5 * (2*NIR + 1 - sqrt((2*NIR + 1)^2 - 8*(NIR - Red)))",
"OSAVI" = "(NIR - Red) / (NIR + Red + 0.16)",
"EVI" = "2.5 * ((NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1))",
"EVI2" = "2.5 * ((NIR - Red) / (NIR + 2.4*Red + 1))",
"DVI" = "NIR - Red",
"RVI" = "NIR / Red",
"GNDVI" = "(NIR - Green) / (NIR + Green)",
"WDVI" = "NIR - 0.5 * Red",
# ==================== ENHANCED/IMPROVED INDICES (12) ====================
"ARVI" = "(NIR - (2*Red - Blue)) / (NIR + (2*Red - Blue))",
"RDVI" = "(NIR - Red) / sqrt(NIR + Red)",
"PVI" = "(NIR - a*Red - b) / sqrt(1 + a^2), where a=1.5, b=10",
"IPVI" = "NIR / (NIR + Red)",
"TNDVI" = "sqrt(((NIR - Red) / (NIR + Red)) + 0.5)",
"GEMI" = "eta * (1 - 0.25*eta) - (Red - 0.125) / (1 - Red), where eta = (2*(NIR^2 - Red^2) + 1.5*NIR + 0.5*Red) / (NIR + Red + 0.5)",
"VARI" = "(Green - Red) / (Green + Red - Blue)",
"TSAVI" = "(a * (NIR - a*Red - b)) / (Red + a*NIR - a*b + X*(1 + a^2)), where a=1.5, b=10, X=0.08",
"ATSAVI" = "(a * (NIR - a*Red - b)) / (a*NIR + Red - a*b + X*(1 + a^2)), where a=1.22, b=0.03, X=0.08",
"GESAVI" = "((NIR - Red) / (NIR + Red + Z)) * (1 + Z), where Z=0.5",
"MTVI" = "1.2 * (1.2*(NIR - Green) - 2.5*(Red - Green))",
"CTVI" = "((NDVI + 0.5) / |NDVI + 0.5|) * sqrt(|NDVI + 0.5|), where NDVI = (NIR - Red) / (NIR + Red)",
# ==================== RED EDGE AND ADVANCED INDICES (10) ====================
"NDRE" = "(NIR - RedEdge) / (NIR + RedEdge)",
"MTCI" = "(RedEdge - Red) / (NIR - Red)",
"IRECI" = "(RedEdge - Red) / (RedEdge / NIR)",
"S2REP" = "705 + 35 * ((Red + RedEdge)/2 - Red) / (RedEdge - Red)",
"PSRI" = "(Red - Green) / RedEdge",
"CRI1" = "(1 / Green) - (1 / Red)",
"CRI2" = "(1 / Green) - (1 / RedEdge)",
"ARI1" = "(1 / Green) - (1 / RedEdge)",
"ARI2" = "NIR * ((1 / Green) - (1 / RedEdge))",
"MCARI" = "((RedEdge - Red) - 0.2*(RedEdge - Green)) * (RedEdge / Red)",
# ==================== STRESS AND CHLOROPHYLL INDICES (12) ====================
"PRI" = "(Green - NIR) / (Green + NIR)",
"SIPI" = "(NIR - Red) / (NIR - Green)",
"CCI" = "(RedEdge - Red) / (RedEdge + Red)",
"NDNI" = "(log(1/NIR) - log(1/SWIR1)) / (log(1/NIR) + log(1/SWIR1))",
"CARI" = "RedEdge * (Red / Green)",
"TCARI" = "3 * ((RedEdge - Red) - 0.2*(RedEdge - Green) * (RedEdge / Red))",
"MTVI1" = "1.2 * (1.2*(NIR - Green) - 2.5*(Red - Green))",
"MTVI2" = "1.5 * (1.2*(NIR - Green) - 2.5*(Red - Green)) / sqrt((2*NIR + 1)^2 - (6*NIR - 5*sqrt(Red)) - 0.5)",
"TVI" = "0.5 * (120*(NIR - Green) - 200*(Red - Green))",
"NPCI" = "(Red - Blue) / (Red + Blue)",
"RARS" = "Red / NIR",
"NPQI" = "(Red - Blue) / (Red + Blue)",
# ==================== WATER AND MOISTURE INDICES (8) ====================
"NDWI" = "(Green - NIR) / (Green + NIR)",
"MNDWI" = "(Green - SWIR1) / (Green + SWIR1)",
"NDMI" = "(NIR - SWIR1) / (NIR + SWIR1)",
"MSI" = "SWIR1 / NIR",
"NDII" = "(NIR - SWIR1) / (NIR + SWIR1)",
"WI" = "NIR / SWIR1",
"SRWI" = "NIR / SWIR1",
"LSWI" = "(NIR - SWIR1) / (NIR + SWIR1)",
# ==================== SPECIALIZED APPLICATIONS (10) ====================
"LAI" = "3.618 * EVI - 0.118, where EVI = 2.5 * ((NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1))",
"FAPAR" = "-0.161 + 1.257 * NDVI, where NDVI = (NIR - Red) / (NIR + Red)",
"FCOVER" = "-2.274 + 4.336*NDVI - 1.33*NDVI^2, where NDVI = (NIR - Red) / (NIR + Red)",
"NBR" = "(NIR - SWIR2) / (NIR + SWIR2)",
"BAI" = "1 / ((0.1 - Red)^2 + (0.06 - NIR)^2)",
"NDSI" = "(Green - SWIR1) / (Green + SWIR1)",
"GRVI" = "(Green - Red) / (Green + Red)",
"VIG" = "(Green - Red) / (Green + Red)",
"CI" = "(Red - Green) / Red",
"GBNDVI" = "(NIR - (Green + Blue)) / (NIR + Green + Blue)"
)
for (idx in indices) {
if (idx %in% names(formula_db)) {
formulas[idx] <- formula_db[[idx]]
} else {
formulas[idx] <- "Formula not available"
}
}
return(formulas)
}
#' Get index typical ranges
#' @keywords internal
get_index_ranges <- function(indices) {
ranges <- setNames(rep("", length(indices)), indices)
range_db <- list(
# ==================== BASIC VEGETATION INDICES (10) ====================
"NDVI" = "[-1, 1]",
"SAVI" = "[-1, 1.5]",
"MSAVI" = "[0, 2]",
"OSAVI" = "[-1, 1]",
"EVI" = "[-1, 3]",
"EVI2" = "[-1, 3]",
"DVI" = "[-2, 2]",
"RVI" = "[0, 30]",
"GNDVI" = "[-1, 1]",
"WDVI" = "[-2, 2]",
# ==================== ENHANCED/IMPROVED INDICES (12) ====================
"ARVI" = "[-1, 1]",
"RDVI" = "[-2, 2]",
"PVI" = "[-2, 2]",
"IPVI" = "[0, 1]",
"TNDVI" = "[0, 1.2]",
"GEMI" = "[-1, 1]",
"VARI" = "[-1, 1]",
"TSAVI" = "[-1, 1.5]",
"ATSAVI" = "[-1, 1.5]",
"GESAVI" = "[-1, 1.5]",
"MTVI" = "[-2, 2]",
"CTVI" = "[-1, 1.5]",
# ==================== RED EDGE AND ADVANCED INDICES (10) ====================
"NDRE" = "[-1, 1]",
"MTCI" = "[0, 8]",
"IRECI" = "[-2, 5]",
"S2REP" = "[680, 780]", # Wavelength in nm
"PSRI" = "[-1, 1]",
"CRI1" = "[-10, 10]",
"CRI2" = "[-10, 10]",
"ARI1" = "[-10, 10]",
"ARI2" = "[-10, 10]",
"MCARI" = "[-2, 2]",
# ==================== STRESS AND CHLOROPHYLL INDICES (12) ====================
"PRI" = "[-1, 1]",
"SIPI" = "[0, 2]",
"CCI" = "[-1, 1]",
"NDNI" = "[-1, 1]",
"CARI" = "[0, 5]",
"TCARI" = "[-2, 5]",
"MTVI1" = "[-2, 2]",
"MTVI2" = "[-2, 5]",
"TVI" = "[-100, 100]",
"NPCI" = "[-1, 1]",
"RARS" = "[0, 5]",
"NPQI" = "[-1, 1]",
# ==================== WATER AND MOISTURE INDICES (8) ====================
"NDWI" = "[-1, 1]",
"MNDWI" = "[-1, 1]",
"NDMI" = "[-1, 1]",
"MSI" = "[0, 5]",
"NDII" = "[-1, 1]",
"WI" = "[0, 10]",
"SRWI" = "[0, 10]",
"LSWI" = "[-1, 1]",
# ==================== SPECIALIZED APPLICATIONS (10) ====================
"LAI" = "[0, 15]",
"FAPAR" = "[0, 1]",
"FCOVER" = "[0, 1]",
"NBR" = "[-1, 1]",
"BAI" = "[0, 1000]",
"NDSI" = "[-1, 1]",
"GRVI" = "[-1, 1]",
"VIG" = "[-1, 1]",
"CI" = "[-1, 1]",
"GBNDVI" = "[-1, 1]"
)
for (idx in indices) {
if (idx %in% names(range_db)) {
ranges[idx] <- range_db[[idx]]
} else {
ranges[idx] <- "Variable"
}
}
return(ranges)
}
#' Get index references
#' @keywords internal
get_index_references <- function(indices) {
references <- setNames(rep("", length(indices)), indices)
ref_db <- list(
# ==================== BASIC VEGETATION INDICES (10) ====================
"NDVI" = "Rouse et al. (1974)",
"SAVI" = "Huete (1988)",
"MSAVI" = "Qi et al. (1994)",
"OSAVI" = "Rondeaux et al. (1996)",
"EVI" = "Huete et al. (1997)",
"EVI2" = "Jiang et al. (2008)",
"DVI" = "Richardson & Wiegand (1977)",
"RVI" = "Birth & McVey (1968)",
"GNDVI" = "Gitelson et al. (1996)",
"WDVI" = "Clevers (1988)",
# ==================== ENHANCED/IMPROVED INDICES (12) ====================
"ARVI" = "Kaufman & Tanre (1992)",
"RDVI" = "Roujean & Breon (1995)",
"PVI" = "Richardson & Wiegand (1977)",
"IPVI" = "Crippen (1990)",
"TNDVI" = "Deering et al. (1975)",
"GEMI" = "Pinty & Verstraete (1992)",
"VARI" = "Gitelson et al. (2002)",
"TSAVI" = "Baret & Guyot (1991)",
"ATSAVI" = "Baret et al. (1992)",
"GESAVI" = "Gilabert et al. (2002)",
"MTVI" = "Haboudane et al. (2004)",
"CTVI" = "Perry & Lautenschlager (1984)",
# ==================== RED EDGE AND ADVANCED INDICES (10) ====================
"NDRE" = "Gitelson & Merzlyak (1994)",
"MTCI" = "Dash & Curran (2004)",
"IRECI" = "Frampton et al. (2013)",
"S2REP" = "Frampton et al. (2013)",
"PSRI" = "Merzlyak et al. (1999)",
"CRI1" = "Gitelson et al. (2002)",
"CRI2" = "Gitelson et al. (2002)",
"ARI1" = "Gitelson et al. (2001)",
"ARI2" = "Gitelson et al. (2001)",
"MCARI" = "Daughtry et al. (2000)",
# ==================== STRESS AND CHLOROPHYLL INDICES (12) ====================
"PRI" = "Gamon et al. (1992)",
"SIPI" = "Penuelas et al. (1995)",
"CCI" = "Barnes et al. (2000)",
"NDNI" = "Serrano et al. (2002)",
"CARI" = "Kim et al. (1994)",
"TCARI" = "Haboudane et al. (2002)",
"MTVI1" = "Haboudane et al. (2004)",
"MTVI2" = "Haboudane et al. (2004)",
"TVI" = "Broge & Leblanc (2001)",
"NPCI" = "Penuelas et al. (1994)",
"RARS" = "Chappelle et al. (1992)",
"NPQI" = "Barnes et al. (1992)",
# ==================== WATER AND MOISTURE INDICES (8) ====================
"NDWI" = "McFeeters (1996)",
"MNDWI" = "Xu (2006)",
"NDMI" = "Gao (1996)",
"MSI" = "Hunt & Rock (1989)",
"NDII" = "Hardisky et al. (1983)",
"WI" = "Gao (1996)",
"SRWI" = "Zarco-Tejada et al. (2003)",
"LSWI" = "Xiao et al. (2002)",
# ==================== SPECIALIZED APPLICATIONS (10) ====================
"LAI" = "Baret & Guyot (1991)",
"FAPAR" = "Myneni & Williams (1994)",
"FCOVER" = "Baret et al. (2007)",
"NBR" = "Lopez Garcia & Caselles (1991)",
"BAI" = "Chuvieco et al. (2002)",
"NDSI" = "Hall et al. (1995)",
"GRVI" = "Tucker (1979)",
"VIG" = "Gitelson et al. (2002)",
"CI" = "Escadafal & Huete (1991)",
"GBNDVI" = "Wang et al. (2007)"
)
for (idx in indices) {
if (idx %in% names(ref_db)) {
references[idx] <- ref_db[[idx]]
} else {
references[idx] <- "Various sources"
}
}
return(references)
}
Try the geospatialsuite package in your browser
Any scripts or data that you put into this service are public.
geospatialsuite documentation built on Nov. 6, 2025, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions get_index_references get_index_ranges get_index_formulas calculate_vegetation_statistics analyze_yield_potential analyze_growth_stage detect_vegetation_stress select_crop_indices validate_output mask_invalid_values apply_temporal_smoothing apply_quality_filter handle_index_edge_cases calculate_index_by_type validate_required_bands load_and_validate_band auto_detect_spectral_bands extract_bands_from_raster get_available_indices load_and_stack_bands load_bands_from_directory analyze_crop_vegetation calculate_multiple_indices list_vegetation_indices calculate_ndvi_enhanced calculate_vegetation_index
Documented in analyze_crop_vegetation analyze_growth_stage analyze_yield_potential apply_quality_filter apply_temporal_smoothing auto_detect_spectral_bands calculate_index_by_type calculate_multiple_indices calculate_ndvi_enhanced calculate_vegetation_index calculate_vegetation_statistics detect_vegetation_stress extract_bands_from_raster get_available_indices get_index_formulas get_index_ranges get_index_references handle_index_edge_cases list_vegetation_indices load_and_stack_bands load_and_validate_band load_bands_from_directory mask_invalid_values select_crop_indices validate_output validate_required_bands
#' Calculate comprehensive vegetation indices
#'
#' @description
#' Calculate a wide range of vegetation indices from spectral bands with
#' automatic band detection, comprehensive error handling, and validation.
#' Supports 40+ different vegetation indices for various applications.
#' Accepts directories, file lists, and automatic CRS handling.
#'
#' @param spectral_data Either individual bands (red, nir, etc.), a multi-band raster,
#' directory path, or list of raster files
#' @param red Red band SpatRaster or file path
#' @param nir NIR band SpatRaster or file path
#' @param blue Optional blue band
#' @param green Optional green band
#' @param swir1 Optional SWIR1 band
#' @param swir2 Optional SWIR2 band
#' @param red_edge Optional Red Edge band
#' @param coastal Optional Coastal/Aerosol band
#' @param nir2 Optional second NIR band
#' @param index_type Vegetation index to calculate (see list_vegetation_indices())
#' @param auto_detect_bands Automatically detect bands from multi-band raster
#' @param band_names Custom band names for multi-band input
#' @param clamp_range Range to clamp output values (optional)
#' @param mask_invalid Mask invalid/extreme values
#' @param scale_factor Scaling factor if needed (default: 1)
#' @param auto_crs_fix Automatically fix CRS mismatches between bands
#' @param verbose Print progress messages
#'
#' @return SpatRaster of vegetation index
#'
#' @details
#' ## Input Format Support:
#'
#' ### Single Calculation:
#' ```r
#' # Individual band files
#' ndvi <- calculate_vegetation_index(red = "red.tif", nir = "nir.tif", index_type = "NDVI")
#'
#' # Multi-band raster
#' evi <- calculate_vegetation_index(spectral_data = "landsat.tif", index_type = "EVI",
#' auto_detect_bands = TRUE)
#' ```
#'
#' ### Directory/Multiple Files:
#' ```r
#' # Directory with band files
#' savi <- calculate_vegetation_index(spectral_data = "/path/to/bands/",
#' band_names = c("red", "nir"), index_type = "SAVI")
#'
#' # File list
#' arvi <- calculate_vegetation_index(spectral_data = c("red.tif", "nir.tif", "blue.tif"),
#' index_type = "ARVI")
#' ```
#'
#' ## Enhanced vs Basic NDVI:
#'
#' ### Basic `calculate_vegetation_index()`:
#' - Single time point calculation
#' - 40+ different indices
#' - Directory/file support
#' - Automatic CRS fixing
#' - **Use for**: Single-date analysis, comparing different indices
#'
#' ### `calculate_ndvi_enhanced()`:
#' - Time series support
#' - Quality filtering
#' - Temporal smoothing
#' - Date matching between red/NIR
#' - **Use for**: Multi-temporal analysis, time series trends
#'
#'#' ## Band Naming Conventions:
#'
#' The function supports case-insensitive band detection:
#' - **Generic names**: "red"/"RED"/"Red", "nir"/"NIR", "blue"/"BLUE", "green"/"GREEN"
#' - **Landsat 8/9**: B1-B7 (e.g., B4=Red, B5=NIR)
#' - **Sentinel-2**: B01-B12 (e.g., B04=Red, B08=NIR, B05=RedEdge)
#' - **MODIS**: band1-band7
#'
#' ### Satellite-Specific Examples:
#'
#' **Landsat 8/9:**
#' ```r
#' # Bands automatically detected
#' ndvi <- calculate_vegetation_index(
#' spectral_data = "LC08_stack.tif", # Has bands named B1-B7
#' index_type = "NDVI",
#' auto_detect_bands = TRUE
#' )
#' ```
#'
#' **Sentinel-2:**
#' ```r
#' # Red Edge indices need Sentinel-2
#' ndre <- calculate_vegetation_index(
#' spectral_data = sentinel_data, # Has B01-B12
#' index_type = "NDRE",
#' auto_detect_bands = TRUE
#' )
#' ```
#'
#' **Custom band names:**
#' ```r
#' # Rename your bands first
#' names(my_raster) <- c("red", "nir", "blue", "green")
#'
#' # Or specify explicitly
#' ndvi <- calculate_vegetation_index(
#' red = my_raster[[1]],
#' nir = my_raster[[4]],
#' index_type = "NDVI"
#' )
#' ```
#'
#' For complete band naming documentation, see:
#' \code{vignette("vegetation-indices", package = "geospatialsuite")}
#'
#' @examples
#' \dontrun{
#' # These examples require satellite imagery files (Landsat/Sentinel data etc.)
#' # Basic NDVI calculation
#' ndvi <- calculate_vegetation_index(red = red_band, nir = nir_band, index_type = "NDVI")
#'
#' # Multi-band raster with auto-detection
#' evi <- calculate_vegetation_index(spectral_data = landsat_stack,
#' index_type = "EVI", auto_detect_bands = TRUE)
#'
#' # Directory with automatic band detection
#' savi <- calculate_vegetation_index(spectral_data = "/path/to/sentinel/bands/",
#' index_type = "SAVI", auto_detect_bands = TRUE)
#'
#' # Advanced index with custom parameters
#' pri <- calculate_vegetation_index(red = red_band, nir = nir_band, green = green_band,
#' index_type = "PRI", clamp_range = c(-1, 1))
#'
#' # Custom band names for multi-band data
#' ndvi <- calculate_vegetation_index(spectral_data = sentinel_data,
#' band_names = c("B4", "B3", "B2", "B8"),
#' index_type = "NDVI")
#' }
#'
#' @export
calculate_vegetation_index <- function(spectral_data = NULL, red = NULL, nir = NULL,
blue = NULL, green = NULL, swir1 = NULL, swir2 = NULL,
red_edge = NULL, coastal = NULL, nir2 = NULL,
index_type = "NDVI", auto_detect_bands = FALSE,
band_names = NULL, clamp_range = NULL,
mask_invalid = TRUE, scale_factor = 1,
auto_crs_fix = TRUE, verbose = FALSE) {
if (verbose) message(sprintf("Starting %s calculation with comprehensive input handling...", index_type))
# Validate index type
available_indices <- get_available_indices()
if (!index_type %in% available_indices$Index) {
stop(sprintf("Index '%s' not supported. Use list_vegetation_indices() to see available indices.",
index_type), call. = FALSE)
}
# Handle input processing
if (!is.null(spectral_data)) {
if (verbose) message("Processing spectral_data input...")
if (is.character(spectral_data)) {
if (length(spectral_data) == 1) {
if (dir.exists(spectral_data)) {
# Directory provided - load and combine bands
if (verbose) message(sprintf("Loading bands from directory: %s", spectral_data))
spectral_data <- load_bands_from_directory(spectral_data, band_names, verbose)
} else if (file.exists(spectral_data)) {
# Single file provided
if (verbose) message(sprintf("Loading multi-band raster: %s", basename(spectral_data)))
spectral_data <- terra::rast(spectral_data)
} else {
stop(sprintf("Spectral data path does not exist: %s", spectral_data), call. = FALSE)
}
} else {
# Multiple files provided - stack them
if (verbose) message(sprintf("Loading and stacking %d band files", length(spectral_data)))
spectral_data <- load_and_stack_bands(spectral_data, band_names, verbose)
}
}
# Extract bands from processed spectral data
bands <- extract_bands_from_raster(spectral_data, auto_detect_bands, band_names, verbose)
red <- bands$red
nir <- bands$nir
blue <- bands$blue
green <- bands$green
swir1 <- bands$swir1
swir2 <- bands$swir2
red_edge <- bands$red_edge
coastal <- bands$coastal
nir2 <- bands$nir2
}
# Input validation and loading with CRS checking
if (verbose) message("Loading and validating input bands with CRS checking...")
red <- load_and_validate_band(red, "red", required = TRUE)
nir <- load_and_validate_band(nir, "nir", required = TRUE)
# CRS compatibility checking
if (verbose) message("Checking and fixing spatial compatibility...")
reference_crs <- terra::crs(red)
# Check and fix NIR CRS
if (!identical(terra::crs(nir), reference_crs)) {
if (auto_crs_fix) {
if (verbose) message("CRS mismatch detected for NIR band. Reprojecting...")
nir <- terra::project(nir, reference_crs)
} else {
stop("CRS mismatch between red and NIR bands. Set auto_crs_fix=TRUE to fix automatically.", call. = FALSE)
}
}
# Align geometries if needed
if (!terra::compareGeom(red, nir, stopOnError = FALSE)) {
if (verbose) message("Geometry mismatch detected. Resampling NIR to match red...")
nir <- terra::resample(nir, red)
}
# Load and align optional bands with CRS fixing
bands_to_process <- list(
blue = blue, green = green, swir1 = swir1, swir2 = swir2,
red_edge = red_edge, coastal = coastal, nir2 = nir2
)
processed_bands <- list()
for (band_name in names(bands_to_process)) {
band_data <- bands_to_process[[band_name]]
if (!is.null(band_data)) {
band_raster <- load_and_validate_band(band_data, band_name, required = FALSE)
if (!is.null(band_raster)) {
# Fix CRS if needed
if (!identical(terra::crs(band_raster), reference_crs) && auto_crs_fix) {
if (verbose) message(sprintf("Reprojecting %s band to match reference CRS", band_name))
band_raster <- terra::project(band_raster, reference_crs)
}
# Align geometry
band_raster <- check_raster_compatibility(red, band_raster, auto_align = TRUE)
processed_bands[[band_name]] <- band_raster
}
}
}
# Assign processed bands
blue <- processed_bands$blue
green <- processed_bands$green
swir1 <- processed_bands$swir1
swir2 <- processed_bands$swir2
red_edge <- processed_bands$red_edge
coastal <- processed_bands$coastal
nir2 <- processed_bands$nir2
# Validate required bands for specific indices
if (!validate_required_bands(index_type, blue, green, swir1, swir2, red_edge, coastal, nir2, verbose)) {
if (verbose) message(sprintf("Skipping %s calculation due to missing required bands", index_type))
# Return NA raster with same dimensions as red band
na_raster <- red
terra::values(na_raster) <- NA
names(na_raster) <- index_type
return(na_raster)
}
# Apply scale factor if needed
if (scale_factor != 1) {
if (verbose) message(sprintf("Applying scale factor: %g", scale_factor))
red <- red * scale_factor
nir <- nir * scale_factor
if (!is.null(blue)) blue <- blue * scale_factor
if (!is.null(green)) green <- green * scale_factor
if (!is.null(swir1)) swir1 <- swir1 * scale_factor
if (!is.null(swir2)) swir2 <- swir2 * scale_factor
if (!is.null(red_edge)) red_edge <- red_edge * scale_factor
if (!is.null(coastal)) coastal <- coastal * scale_factor
if (!is.null(nir2)) nir2 <- nir2 * scale_factor
}
# Calculate vegetation index
if (verbose) message(sprintf("Computing %s index...", index_type))
index <- calculate_index_by_type(index_type, red, nir, blue, green, swir1, swir2,
red_edge, coastal, nir2, verbose)
# Apply clamp range if specified
if (!is.null(clamp_range)) {
if (verbose) message(sprintf("Clamping values to range [%.3f, %.3f]",
clamp_range[1], clamp_range[2]))
index <- terra::clamp(index, clamp_range[1], clamp_range[2])
}
# Mask invalid values if requested
if (mask_invalid) {
index <- mask_invalid_values(index, index_type, verbose)
}
# Final validation and reporting
validation_result <- validate_output(index, index_type, verbose)
if (!validation_result) {
# If validation failed, return an NA raster but don't stop the process
if (verbose) message(sprintf("Returning NA raster for failed index: %s", index_type))
na_raster <- red # Use red band as template
terra::values(na_raster) <- NA
names(na_raster) <- index_type
return(na_raster)
}
names(index) <- index_type
return(index)
}
#' Calculate NDVI with time series options
#'
#' @description
#' NDVI calculation specifically designed for time series analysis with
#' date matching, quality filtering, temporal smoothing, and multi-temporal support.
#' **Use this for time series analysis, use calculate_vegetation_index() for single dates.**
#'
#' @param red_data Red band data (files, directory, or raster objects)
#' @param nir_data NIR band data (files, directory, or raster objects)
#' @param clamp_range Range to clamp NDVI values (default: c(-0.2, 1))
#' @param match_by_date Logical: match rasters by date using filenames
#' @param quality_filter Apply quality filtering (remove outliers)
#' @param temporal_smoothing Apply temporal smoothing for time series
#' @param verbose Print progress messages
#' @param date_patterns Custom date patterns for matching
#'
#' @return SpatRaster with NDVI layers (single or multi-layer for time series)
#'
#' @details
#' ## When to Use Enhanced vs Basic NDVI:
#'
#' ### Use `calculate_ndvi_enhanced()` for:
#' - **Time series analysis**: Multiple dates, trend analysis
#' - **Quality control**: Remove outliers, temporal smoothing
#' - **Date matching**: Automatic pairing of red/NIR by date
#' - **Multi-temporal studies**: Seasonal analysis, change detection
#'
#' ### Use `calculate_vegetation_index(index_type="NDVI")` for:
#' - **Single date analysis**: One-time calculation
#' - **Different indices**: Want to calculate EVI, SAVI, etc. too
#' - **Quick calculations**: Simple, fast NDVI
#' - **Mixed workflows**: Part of larger vegetation index analysis
#'
#' @examples
#' \dontrun{
#' # These examples require external data files not included with the package
#' # Time series NDVI with date matching
#' ndvi_series <- calculate_ndvi_enhanced(
#' red_data = "/path/to/red/time_series/",
#' nir_data = "/path/to/nir/time_series/",
#' match_by_date = TRUE,
#' quality_filter = TRUE,
#' temporal_smoothing = TRUE
#' )
#'
#' # Simple NDVI (single date with quality control)
#' ndvi_clean <- calculate_ndvi_enhanced(
#' red_data = red_raster,
#' nir_data = nir_raster,
#' quality_filter = TRUE
#' )
#' }
#'
#' @export
calculate_ndvi_enhanced <- function(red_data, nir_data, clamp_range = c(-0.2, 1),
match_by_date = FALSE, quality_filter = FALSE,
temporal_smoothing = FALSE, verbose = FALSE,
date_patterns = NULL) {
# Input validation
if (missing(red_data) || missing(nir_data)) {
stop("Both red_data and nir_data are required", call. = FALSE)
}
if (is.null(red_data) || is.null(nir_data)) {
stop("red_data and nir_data cannot be NULL", call. = FALSE)
}
if (verbose) message("Starting NDVI calculation for time series analysis...")
# Handle different input types with CRS support
if (inherits(red_data, "SpatRaster") && inherits(nir_data, "SpatRaster")) {
# Direct raster objects provided - use them directly
if (verbose) message("Using provided SpatRaster objects directly")
# CRS and geometry checking
if (!identical(terra::crs(red_data), terra::crs(nir_data))) {
if (verbose) message("CRS mismatch detected. Reprojecting NIR to match RED...")
nir_data <- terra::project(nir_data, terra::crs(red_data))
}
if (!terra::compareGeom(red_data, nir_data, stopOnError = FALSE)) {
if (verbose) message("Geometry mismatch detected. Resampling NIR to match RED...")
nir_data <- terra::resample(nir_data, red_data)
}
# Calculate NDVI directly
ndvi <- (nir_data - red_data) / (nir_data + red_data)
# Handle division by zero
zero_sum_mask <- (red_data + nir_data) == 0
if (any(terra::values(zero_sum_mask, mat = FALSE), na.rm = TRUE)) {
ndvi[zero_sum_mask] <- NA
}
# Apply clamp range
if (!is.null(clamp_range)) {
ndvi <- terra::clamp(ndvi, clamp_range[1], clamp_range[2])
}
# Quality filtering
if (quality_filter) {
ndvi <- apply_quality_filter(ndvi, verbose)
}
names(ndvi) <- "NDVI_enhanced"
if (verbose) {
message("NDVI calculation completed successfully!")
message(sprintf(" - %d layer(s) created", terra::nlyr(ndvi)))
valid_values <- terra::values(ndvi, mat = FALSE)
valid_values <- valid_values[!is.na(valid_values)]
if (length(valid_values) > 0) {
message(sprintf(" - Value range: [%.3f, %.3f]",
min(valid_values), max(valid_values)))
}
}
return(ndvi)
}
# Handle file/directory inputs with error handling
if (match_by_date && (is.character(red_data) && is.character(nir_data))) {
if (verbose) message("Matching red and NIR bands by date...")
tryCatch({
matched_df <- match_rasters_by_date(red_data, nir_data, date_patterns, verbose)
red_files <- matched_df$red
nir_files <- matched_df$nir
layer_dates <- matched_df$date
}, error = function(e) {
if (verbose) message("Date matching failed, using direct file order")
red_files <- red_data
nir_files <- nir_data
layer_dates <- paste0("Layer_", seq_along(if(length(red_data) > 1) red_data else c(red_data)))
})
} else {
red_files <- red_data
nir_files <- nir_data
layer_dates <- paste0("Layer_", seq_along(if(length(red_data) > 1) red_data else c(red_data)))
}
# Load rasters with error handling
tryCatch({
red_rasters <- load_raster_data(red_files, verbose = verbose)
nir_rasters <- load_raster_data(nir_files, verbose = verbose)
}, error = function(e) {
stop(sprintf("Failed to load raster data: %s", e$message), call. = FALSE)
})
# Ensure we have matching numbers of rasters
if (length(red_rasters) != length(nir_rasters)) {
min_length <- min(length(red_rasters), length(nir_rasters))
warning(sprintf("Mismatched number of red (%d) and NIR (%d) rasters. Using first %d pairs.",
length(red_rasters), length(nir_rasters), min_length))
red_rasters <- red_rasters[1:min_length]
nir_rasters <- nir_rasters[1:min_length]
layer_dates <- layer_dates[1:min_length]
}
# Calculate NDVI for each pair with CRS handling
ndvi_list <- vector("list", length(red_rasters))
for (i in seq_along(red_rasters)) {
if (verbose && (i %% 5 == 0 || i == length(red_rasters))) {
message(sprintf("Processing pair %d/%d (%s)...", i, length(red_rasters), layer_dates[i]))
}
tryCatch({
red <- red_rasters[[i]]
nir <- nir_rasters[[i]]
# CRS and geometry compatibility
if (!identical(terra::crs(red), terra::crs(nir))) {
if (verbose && i == 1) message("CRS mismatch detected in time series. Reprojecting...")
nir <- terra::project(nir, terra::crs(red))
}
if (!terra::compareGeom(red, nir, stopOnError = FALSE)) {
nir <- terra::resample(nir, red)
}
# Calculate NDVI
ndvi <- (nir - red) / (nir + red)
# Handle division by zero
zero_sum_mask <- (red + nir) == 0
if (any(terra::values(zero_sum_mask, mat = FALSE), na.rm = TRUE)) {
ndvi[zero_sum_mask] <- NA
}
# Apply clamp range
if (!is.null(clamp_range)) {
ndvi <- terra::clamp(ndvi, clamp_range[1], clamp_range[2])
}
# Quality filtering
if (quality_filter) {
ndvi <- apply_quality_filter(ndvi, verbose = FALSE) # Suppress per-layer messages
}
names(ndvi) <- paste0("NDVI_", gsub("[^A-Za-z0-9_]", "_", layer_dates[i]))
ndvi_list[[i]] <- ndvi
}, error = function(e) {
warning(sprintf("Failed to process pair %d: %s", i, e$message))
# Create an NA raster with same dimensions as the first successful one
if (i > 1 && !is.null(ndvi_list[[1]])) {
na_raster <- ndvi_list[[1]]
terra::values(na_raster) <- NA
names(na_raster) <- paste0("NDVI_", gsub("[^A-Za-z0-9_]", "_", layer_dates[i]))
ndvi_list[[i]] <- na_raster
} else {
ndvi_list[[i]] <- NULL
}
})
}
# Remove NULL entries
ndvi_list <- Filter(Negate(is.null), ndvi_list)
if (length(ndvi_list) == 0) {
stop("No NDVI layers could be calculated successfully", call. = FALSE)
}
# Combine results
result <- if (length(ndvi_list) == 1) {
ndvi_list[[1]]
} else {
tryCatch({
terra::rast(ndvi_list)
}, error = function(e) {
warning("Could not combine layers into stack, returning first layer only")
ndvi_list[[1]]
})
}
# Apply temporal smoothing if requested and we have multiple layers
if (temporal_smoothing && terra::nlyr(result) > 1) {
if (verbose) message("Applying temporal smoothing...")
tryCatch({
result <- apply_temporal_smoothing(result, verbose)
}, error = function(e) {
warning(sprintf("Temporal smoothing failed: %s", e$message))
})
}
if (verbose) {
message("NDVI time series calculation completed successfully!")
message(sprintf(" - %d layer(s) created", terra::nlyr(result)))
if (terra::nlyr(result) > 0) {
valid_values <- terra::values(result, mat = FALSE)
valid_values <- valid_values[!is.na(valid_values)]
if (length(valid_values) > 0) {
message(sprintf(" - Value range: [%.3f, %.3f]",
min(valid_values), max(valid_values)))
}
}
}
return(result)
}
#' Get comprehensive list of available vegetation indices
#'
#' @description
#' Returns detailed information about all 40+ available vegetation indices including
#' formulas, required bands, applications, and references.
#'
#' @param category Filter by category: "all", "basic", "enhanced", "specialized", "stress"
#' @param application Filter by application: "general", "agriculture", "forestry", "stress", "water"
#' @param detailed Return detailed information including formulas and references
#'
#' @return Data frame with vegetation index information
#'
#' @examples
#' \donttest{
#' # All available indices
#' all_indices <- list_vegetation_indices()
#'
#' # Only stress detection indices
#' stress_indices <- list_vegetation_indices(category = "stress")
#'
#' # Detailed information with formulas
#' detailed_info <- list_vegetation_indices(detailed = TRUE)
#'
#' # Agricultural applications only
#' ag_indices <- list_vegetation_indices(application = "agriculture")
#' }
#'
#' @export
list_vegetation_indices <- function(category = "all", application = "all", detailed = FALSE) {
# Comprehensive index database - 40+ indices
indices_db <- data.frame(
Index = c(
# 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",
# Water and Moisture Indices (8)
"NDWI", "MNDWI", "NDMI", "MSI", "NDII", "WI", "SRWI", "LSWI",
# Specialized Applications (10)
"LAI", "FAPAR", "FCOVER", "NBR", "BAI", "NDSI", "GRVI", "VIG", "CI", "GBNDVI"
),
Category = c(
# Basic (10)
rep("basic", 10),
# Enhanced (12)
rep("enhanced", 12),
# Red Edge (10)
rep("advanced", 10),
# Stress (12)
rep("stress", 12),
# Water (8)
rep("water", 8),
# Specialized (10)
rep("specialized", 10)
),
Application = c(
# Basic applications (10)
"general", "soil", "soil", "soil", "general", "general", "biomass", "general", "chlorophyll", "soil",
# Enhanced applications (12)
"atmospheric", "general", "general", "general", "general", "general", "general", "soil", "soil", "soil", "general", "canopy",
# Advanced applications (10)
"stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress",
# Stress applications (12)
"stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress", "stress",
# Water applications (8)
"water", "water", "water", "water", "water", "water", "water", "water",
# Specialized applications (10)
"forestry", "forestry", "forestry", "forestry", "forestry", "snow", "general", "general", "canopy", "general"
),
Required_Bands = c(
# Basic indices (10)
"Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR, Blue",
"Red, NIR", "Red, NIR", "Red, NIR", "Green, NIR", "Red, NIR",
# Enhanced indices (12)
"Red, NIR, Blue", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR",
"Red, NIR", "Red, Green, Blue", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR", "Red, NIR",
# Advanced indices (10)
"NIR, RedEdge", "RedEdge, NIR", "RedEdge, NIR", "RedEdge", "RedEdge, NIR",
"Red, Green", "RedEdge, Green", "RedEdge, Green", "RedEdge, NIR", "Red, Green",
# Stress indices (12)
"Green, NIR", "Red, NIR", "RedEdge, Green", "NIR, SWIR1", "Red, Green",
"Red, Green", "Red, NIR", "Red, NIR", "Red, NIR", "Red, Blue", "Red, NIR", "Red, Blue",
# Water indices (8)
"Green, NIR", "Green, SWIR1", "NIR, SWIR1", "NIR, SWIR1", "NIR, SWIR1",
"NIR, SWIR1", "NIR, SWIR1", "NIR, SWIR1",
# Specialized indices (10)
"Red, NIR", "Red, NIR", "Red, NIR", "NIR, SWIR2", "Red, NIR",
"Green, SWIR1", "Red, Green", "Green, NIR", "Red, Green", "Green, Blue, NIR"
),
stringsAsFactors = FALSE
)
# Add descriptions
indices_db$Description <- c(
# Basic (10)
"Normalized Difference Vegetation Index", "Soil Adjusted Vegetation Index",
"Modified Soil Adjusted Vegetation Index", "Optimized Soil Adjusted Vegetation Index",
"Enhanced Vegetation Index", "Two-band Enhanced Vegetation Index",
"Difference Vegetation Index", "Ratio Vegetation Index", "Green NDVI", "Weighted Difference Vegetation Index",
# Enhanced (12)
"Atmospherically Resistant Vegetation Index", "Renormalized Difference Vegetation Index",
"Perpendicular Vegetation Index", "Infrared Percentage Vegetation Index", "Transformed NDVI",
"Global Environment Monitoring Index", "Visible Atmospherically Resistant Index", "Transformed Soil Adjusted VI",
"Adjusted Transformed Soil Adjusted VI", "Generalized Soil Adjusted VI", "Modified Triangular VI", "Corrected Transformed VI",
# Advanced (10)
"Normalized Difference Red Edge", "MERIS Terrestrial Chlorophyll Index",
"Inverted Red-Edge Chlorophyll Index", "Sentinel-2 Red-Edge Position",
"Plant Senescence Reflectance Index", "Carotenoid Reflectance Index 1",
"Carotenoid Reflectance Index 2", "Anthocyanin Reflectance Index 1",
"Anthocyanin Reflectance Index 2", "Modified Chlorophyll Absorption Ratio Index",
# Stress (12)
"Photochemical Reflectance Index", "Structure Insensitive Pigment Index",
"Canopy Chlorophyll Index", "Normalized Difference Nitrogen Index",
"Chlorophyll Absorption Ratio Index", "Transformed Chlorophyll Absorption Ratio Index",
"Modified Triangular Vegetation Index 1", "Modified Triangular Vegetation Index 2",
"Triangular Vegetation Index", "Normalized Pigment Chlorophyll Index",
"Ratio Analysis of Reflectance Spectra", "Normalized Phaeophytinization Index",
# Water (8)
"Normalized Difference Water Index", "Modified Normalized Difference Water Index",
"Normalized Difference Moisture Index", "Moisture Stress Index",
"Normalized Difference Infrared Index", "Water Index",
"Simple Ratio Water Index", "Land Surface Water Index",
# Specialized (10)
"Leaf Area Index", "Fraction of Absorbed PAR", "Fraction of Vegetation Cover",
"Normalized Burn Ratio", "Burn Area Index", "Normalized Difference Snow Index",
"Green-Red Vegetation Index", "Vegetation Index Green", "Coloration Index",
"Green-Blue NDVI"
)
# Filter by category
if (category != "all") {
indices_db <- indices_db[indices_db$Category == category, ]
}
# Filter by application
if (application != "all") {
indices_db <- indices_db[indices_db$Application == application, ]
}
# Add detailed information if requested
if (detailed) {
indices_db$Formula <- get_index_formulas(indices_db$Index)
indices_db$Range <- get_index_ranges(indices_db$Index)
indices_db$Reference <- get_index_references(indices_db$Index)
}
return(indices_db)
}
#' Calculate multiple vegetation indices at once
#'
#' @description
#' Calculate multiple vegetation indices from the same spectral data in a single operation.
#' Efficient for comparative analysis and comprehensive vegetation assessment.
#' Supports directory input and automatic CRS handling.
#'
#' @param spectral_data Multi-band raster, directory path, or individual bands
#' @param indices Vector of index names to calculate
#' @param output_stack Return as single multi-layer raster (TRUE) or list (FALSE)
#' @param region_boundary Optional region boundary for clipping
#' @param parallel Use parallel processing for multiple indices
#' @param verbose Print progress messages
#' @param ... Additional arguments passed to calculate_vegetation_index
#'
#' @return SpatRaster stack or list of indices
#'
#' @examples
#' \dontrun{
#' # These examples require satellite imagery files (Landsat/Sentinel data etc.)
#' # Calculate multiple basic indices from directory
#' multi_indices <- calculate_multiple_indices(
#' spectral_data = "/path/to/sentinel/bands/",
#' indices = c("NDVI", "EVI", "SAVI", "MSAVI"),
#' auto_detect_bands = TRUE
#' )
#'
#' # Comprehensive vegetation analysis from individual files
#' veg_analysis <- calculate_multiple_indices(
#' red = red_band, nir = nir_band, blue = blue_band,
#' indices = c("NDVI", "EVI", "ARVI", "GNDVI", "DVI"),
#' output_stack = TRUE,
#' region_boundary = "Iowa"
#' )
#'
#' # Directory with custom band matching
#' stress_indices <- calculate_multiple_indices(
#' spectral_data = "/path/to/bands/",
#' indices = c("PRI", "SIPI", "NDRE"),
#' band_names = c("red", "green", "nir", "red_edge"),
#' output_stack = TRUE
#' )
#' }
#'
#' @export
calculate_multiple_indices <- function(spectral_data = NULL, indices = c("NDVI", "EVI", "SAVI"),
output_stack = TRUE, region_boundary = NULL,
parallel = FALSE, verbose = FALSE, ...) {
if (verbose) message(sprintf("Calculating %d vegetation indices...", length(indices)))
# Validate indices
available_indices <- get_available_indices()$Index
invalid_indices <- setdiff(indices, available_indices)
if (length(invalid_indices) > 0) {
stop(sprintf("Invalid indices: %s", paste(invalid_indices, collapse = ", ")))
}
# Calculate each index with error handling
index_results <- list()
if (parallel && requireNamespace("parallel", quietly = TRUE)) {
if (verbose) message("Using parallel processing...")
# Setup parallel processing
n_cores <- min(parallel::detectCores() - 1, length(indices))
if (verbose) message(sprintf("Using %d cores", n_cores))
index_results <- parallel::mclapply(indices, function(idx) {
tryCatch({
calculate_vegetation_index(spectral_data = spectral_data,
index_type = idx, verbose = FALSE, ...)
}, error = function(e) {
warning(sprintf("Failed to calculate %s: %s", idx, e$message))
return(NULL)
})
}, mc.cores = n_cores)
# Remove NULL results
failed_indices <- names(index_results)[sapply(index_results, is.null)]
index_results <- index_results[!sapply(index_results, is.null)]
if (verbose && length(failed_indices) > 0) {
message(sprintf("Failed indices (parallel): %s", paste(failed_indices, collapse = ", ")))
}
} else {
# Sequential processing with better error handling
failed_indices <- c()
for (i in seq_along(indices)) {
idx <- indices[i]
if (verbose) message(sprintf("Calculating %s (%d/%d)...", idx, i, length(indices)))
tryCatch({
result <- calculate_vegetation_index(
spectral_data = spectral_data,
index_type = idx, verbose = FALSE, ...
)
# Check if result is valid (not all NA)
if (all(is.na(terra::values(result, mat = FALSE)))) {
warning(sprintf("Index %s produced only NA values - skipping", idx))
failed_indices <- c(failed_indices, idx)
} else {
index_results[[idx]] <- result
}
}, error = function(e) {
warning(sprintf("Failed to calculate %s: %s", idx, e$message))
failed_indices <- c(failed_indices, idx)
})
}
if (verbose && length(failed_indices) > 0) {
message(sprintf("Failed/skipped indices: %s", paste(failed_indices, collapse = ", ")))
message(sprintf("Successfully calculated: %s", paste(names(index_results), collapse = ", ")))
}
}
# Ensure we have at least some results
if (length(index_results) == 0) {
stop("No vegetation indices could be calculated successfully", call. = FALSE)
}
names(index_results) <- names(index_results) %||% indices[1:length(index_results)]
# Get boundary
if (!is.null(region_boundary)) {
if (verbose) message("Applying region boundary...")
boundary <- get_region_boundary(region_boundary)
boundary_vect <- terra::vect(boundary)
index_results <- lapply(index_results, function(idx_raster) {
if (!is.null(idx_raster) && inherits(idx_raster, "SpatRaster")) {
cropped <- terra::crop(idx_raster, boundary_vect)
terra::mask(cropped, boundary_vect)
} else {
idx_raster
}
})
}
# Return as stack or list - FIXED VERSION
if (output_stack && length(index_results) > 0) {
if (verbose) message("Creating multi-layer stack...")
# Filter out NULL or invalid results
valid_results <- index_results[!sapply(index_results, is.null)]
valid_results <- valid_results[sapply(valid_results, function(x) inherits(x, "SpatRaster"))]
if (length(valid_results) == 0) {
stop("No valid indices to create stack", call. = FALSE)
}
tryCatch({
# Create stack from valid results
result <- terra::rast(valid_results)
# Set names properly - use the names from the list
names(result) <- names(valid_results)
if (verbose) {
message(sprintf("Successfully created stack with %d indices: %s",
terra::nlyr(result), paste(names(result), collapse = ", ")))
}
return(result)
}, error = function(e) {
warning(sprintf("Failed to create raster stack: %s. Returning list instead.", e$message))
return(valid_results)
})
} else {
# Return as list
valid_results <- index_results[!sapply(index_results, is.null)]
if (verbose && length(valid_results) > 0) {
message(sprintf("Returning list with %d indices: %s",
length(valid_results), paste(names(valid_results), collapse = ", ")))
}
return(valid_results)
}
}
#' Specialized crop vegetation analysis
#'
#' @description
#' Perform comprehensive vegetation analysis specifically designed for crop monitoring
#' including growth stage detection, stress identification, and yield prediction support.
#' Handles test scenarios properly with better input validation.
#'
#' @param spectral_data Multi-band spectral data (file, directory, or SpatRaster)
#' @param crop_type Crop type for specialized analysis ("corn", "soybeans", "wheat", "general")
#' @param growth_stage Growth stage if known ("early", "mid", "late", "harvest")
#' @param analysis_type Type of analysis: "comprehensive", "stress", "growth", "yield"
#' @param cdl_mask Optional CDL mask for crop-specific analysis
#' @param reference_data Optional reference data for validation
#' @param output_folder Optional output folder for results
#' @param verbose Print detailed progress
#'
#' @return List with comprehensive vegetation analysis results:
#' \itemize{
#' \item \code{vegetation_indices}: SpatRaster with calculated indices
#' \item \code{analysis_results}: Detailed analysis results by type
#' \item \code{metadata}: Analysis metadata and parameters
#' }
#'
#' @details
#' ## Crop-Specific Index Selection:
#' - **Corn**: NDVI, EVI, GNDVI, DVI, RVI, PRI
#' - **Soybeans**: NDVI, EVI, SAVI, GNDVI, PRI
#' - **Wheat**: NDVI, EVI, SAVI, DVI
#' - **General**: NDVI, EVI, SAVI, GNDVI, DVI, RVI
#'
#' ## Analysis Types:
#' - **comprehensive**: All analyses (stress, growth, yield)
#' - **stress**: Focus on stress detection indices
#' - **growth**: Growth stage analysis
#' - **yield**: Yield prediction support
#'
#' ## Output Structure:
#'
#' The function returns a list with three main components:
#'
#' ### 1. vegetation_indices (SpatRaster)
#' Multi-layer raster with calculated indices (NDVI, EVI, etc.)
#'
#' ### 2. analysis_results (List)
#'
#' **stress_analysis** (if requested):
#' - Percentage of pixels in each stress category
#' - Categories: healthy (NDVI 0.6-1.0), moderate stress (0.4-0.6), severe stress (0.0-0.4)
#' - Includes mean, median, std_dev, and thresholds used
#'
#' **growth_analysis** (if requested):
#' - Predicted growth stage based on NDVI patterns
#' - Stage confidence (0-1 scale)
#' - Detailed statistics for each index
#' - Growth stages: emergence, vegetative, reproductive, maturity (crop-specific)
#'
#' **yield_analysis** (if requested):
#' - **Composite Yield Index**: Normalized 0-1 score combining multiple indices
#' - 0.0 = Very low yield potential
#' - 0.5 = Medium yield potential
#' - 1.0 = Maximum yield potential
#' - **Yield Potential Class**: Categorical (Low, Medium, High, Very High)
#' - **Index Contributions**: How each index contributed to composite score
#' - Calculation: Each index (NDVI, EVI, GNDVI, DVI, RVI) is normalized to 0-1,
#' then averaged to create composite score
#'
#' **summary_statistics**:
#' - Basic stats (mean, std, min, max, percentiles) for all indices
#' - Coverage percentage and pixel counts
#'
#' ### 3. metadata (List)
#' Processing information: crop_type, indices_used, processing_date, spatial properties
#'
#' ## Example Interpretation:
#'
#' ```r
#' result <- analyze_crop_vegetation(data, crop_type = "corn")
#'
#' # Stress Assessment
#' stress <- result$analysis_results$stress_analysis$NDVI
#' cat(sprintf("Healthy: %.1f%%, Stressed: %.1f%%\n",
#' stress$healthy_percentage,
#' stress$severe_stress_percentage))
#'
#' # Growth Stage
#' stage <- result$analysis_results$growth_analysis$predicted_growth_stage
#' cat(sprintf("Growth stage: %s\n", stage))
#'
#' # Yield Potential
#' yield <- result$analysis_results$yield_analysis
#' cat(sprintf("Yield potential: %s (score: %.2f)\n",
#' yield$yield_potential_class,
#' yield$composite_yield_index))
#' ```
#'
#' ## Important Notes:
#'
#' - **Composite Yield Index** is a vegetation-based proxy, not a direct yield prediction
#' - Thresholds are based on literature and may need regional calibration
#' - Results should be validated with ground truth data
#' - For detailed output documentation, see package vignette
#'
#' ## Analysis Types:
#' - **comprehensive**: All analyses (stress, growth, yield)
#' - **stress**: Focus on stress detection indices
#' - **growth**: Growth stage analysis
#' - **yield**: Yield prediction support
#'
#' @examples
#' \dontrun{
#' # These examples require actual spectral data
#' # Comprehensive corn analysis
#' corn_analysis <- analyze_crop_vegetation(
#' spectral_data = sentinel_data,
#' crop_type = "corn",
#' analysis_type = "comprehensive",
#' cdl_mask = corn_mask
#' )
#'
#' # Access results
#' corn_analysis$vegetation_indices # SpatRaster with indices
#' corn_analysis$analysis_results$stress_analysis # Stress detection results
#' corn_analysis$metadata$indices_used # Which indices were calculated
#'
#' # Stress detection in soybeans
#' stress_analysis <- analyze_crop_vegetation(
#' spectral_data = landsat_stack,
#' crop_type = "soybeans",
#' analysis_type = "stress",
#' growth_stage = "mid"
#' )
#' }
#'
#' \donttest{
#' # Example with mock spectral data
#' # Create mock multi-band raster (simulating satellite data)
#' red_band <- terra::rast(nrows = 5, ncols = 5, crs = "EPSG:4326")
#' nir_band <- terra::rast(nrows = 5, ncols = 5, crs = "EPSG:4326")
#' terra::values(red_band) <- runif(25, 0.1, 0.3) # Typical red values
#' terra::values(nir_band) <- runif(25, 0.4, 0.8) # Typical NIR values
#' spectral_stack <- c(red_band, nir_band)
#' names(spectral_stack) <- c("red", "nir")
#'
#' # Analyze with mock data
#' result <- analyze_crop_vegetation(spectral_stack, crop_type = "general")
#' print(names(result)) # Should show analysis components
#' }
#'
#' @export
analyze_crop_vegetation <- function(spectral_data, crop_type = "general",
growth_stage = "unknown", analysis_type = "comprehensive",
cdl_mask = NULL, reference_data = NULL,
output_folder = NULL, verbose = FALSE) {
if (verbose) message(sprintf("Starting specialized %s vegetation analysis for %s...",
analysis_type, crop_type))
# Input validation and loading
if (is.character(spectral_data)) {
if (length(spectral_data) == 1) {
if (dir.exists(spectral_data)) {
if (verbose) message(sprintf("Loading spectral data from directory: %s", spectral_data))
spectral_data <- load_bands_from_directory(spectral_data, verbose = verbose)
} else if (file.exists(spectral_data)) {
if (verbose) message(sprintf("Loading spectral data from file: %s", basename(spectral_data)))
spectral_data <- terra::rast(spectral_data)
} else {
stop(sprintf("Spectral data path does not exist: %s", spectral_data), call. = FALSE)
}
} else {
if (verbose) message("Loading and stacking multiple spectral files")
spectral_data <- load_and_stack_bands(spectral_data, verbose = verbose)
}
}
# Validate spectral data
if (!inherits(spectral_data, "SpatRaster")) {
stop("spectral_data must be a SpatRaster object, file path, or directory", call. = FALSE)
}
if (terra::nlyr(spectral_data) < 2) {
stop("spectral_data must have at least 2 bands for vegetation analysis", call. = FALSE)
}
# Select indices based on crop type and analysis
selected_indices <- select_crop_indices(crop_type, analysis_type, growth_stage)
if (verbose) {
message(sprintf("Selected %d indices for analysis: %s",
length(selected_indices), paste(selected_indices, collapse = ", ")))
}
# Apply CDL mask if provided
if (!is.null(cdl_mask)) {
if (verbose) message("Applying crop mask...")
spectral_data <- terra::mask(spectral_data, cdl_mask)
}
# Calculate vegetation indices using the main function
vegetation_indices <- calculate_multiple_indices(
spectral_data = spectral_data,
indices = selected_indices,
output_stack = TRUE,
auto_detect_bands = TRUE,
verbose = verbose
)
# Perform analysis based on type
analysis_results <- list()
if (analysis_type %in% c("comprehensive", "stress")) {
# Stress detection analysis
if (verbose) message("Performing stress detection analysis...")
analysis_results$stress_analysis <- detect_vegetation_stress(
vegetation_indices, crop_type, selected_indices, verbose
)
}
if (analysis_type %in% c("comprehensive", "growth")) {
# Growth stage analysis
if (verbose) message("Performing growth stage analysis...")
analysis_results$growth_analysis <- analyze_growth_stage(
vegetation_indices, crop_type, growth_stage, verbose
)
}
if (analysis_type %in% c("comprehensive", "yield")) {
# Yield prediction support
if (verbose) message("Performing yield prediction analysis...")
analysis_results$yield_analysis <- analyze_yield_potential(
vegetation_indices, crop_type, selected_indices, verbose
)
}
# Generate summary statistics
analysis_results$summary_statistics <- calculate_vegetation_statistics(
vegetation_indices, selected_indices, verbose
)
# Validation if reference data provided
if (!is.null(reference_data)) {
if (verbose) message("Performing validation analysis...")
analysis_results$validation <- validate_vegetation_analysis(
vegetation_indices, reference_data, verbose
)
}
# Save results if output folder specified
if (!is.null(output_folder)) {
save_vegetation_analysis_results(analysis_results, vegetation_indices,
output_folder, crop_type, verbose)
}
# Combine all results
final_results <- list(
vegetation_indices = vegetation_indices,
analysis_results = analysis_results,
metadata = list(
crop_type = crop_type,
growth_stage = growth_stage,
analysis_type = analysis_type,
indices_used = selected_indices,
processing_date = Sys.time(),
input_bands = terra::nlyr(spectral_data),
spatial_resolution = terra::res(spectral_data),
spatial_extent = as.vector(terra::ext(spectral_data))
)
)
if (verbose) message("Crop vegetation analysis completed successfully!")
return(final_results)
}
# ==================== HELPER FUNCTIONS ==================== #
#' Load bands from directory
#' @keywords internal
load_bands_from_directory <- function(directory, band_names = NULL, verbose = FALSE) {
# Find all raster files
files <- list.files(directory, pattern = "\\.(tif|tiff|img)$",
full.names = TRUE, ignore.case = TRUE)
if (length(files) == 0) {
stop(sprintf("No raster files found in directory: %s", directory), call. = FALSE)
}
if (verbose) message(sprintf("Found %d raster files in directory", length(files)))
# Try to auto-detect which files are which bands
if (is.null(band_names)) {
# Use standard patterns to identify bands
band_patterns <- list(
red = "red|RED|Red|B04|b04|B4|b4",
nir = "nir|NIR|Nir|B08|b08|B8|b8|B05|b05|B5|b5",
blue = "blue|BLUE|Blue|B02|b02|B2|b2",
green = "green|GREEN|Green|B03|b03|B3|b3",
swir1 = "swir1|SWIR1|B11|b11|B6|b6",
swir2 = "swir2|SWIR2|B12|b12|B7|b7",
red_edge = "rededge|RedEdge|red_edge|B05|b05"
)
band_files <- list()
for (band_type in names(band_patterns)) {
pattern <- band_patterns[[band_type]]
matches <- files[grepl(pattern, basename(files), ignore.case = TRUE)]
if (length(matches) > 0) {
band_files[[band_type]] <- matches[1] # Take first match
if (verbose) message(sprintf(" Detected %s band: %s", band_type, basename(matches[1])))
}
}
if (length(band_files) == 0) {
stop("Could not auto-detect any bands. Please specify band_names parameter.", call. = FALSE)
}
# Stack the detected bands
band_rasters <- lapply(band_files, terra::rast)
stacked <- terra::rast(band_rasters)
names(stacked) <- names(band_files)
} else {
# Use provided band names to match files
if (verbose) message("Using provided band names for file matching")
band_files <- list()
for (i in seq_along(band_names)) {
band_name <- band_names[i]
# Look for files containing the band name
matches <- files[grepl(band_name, basename(files), ignore.case = TRUE)]
if (length(matches) > 0) {
band_files[[band_name]] <- matches[1]
if (verbose) message(sprintf(" Matched %s: %s", band_name, basename(matches[1])))
} else {
warning(sprintf("No file found for band: %s", band_name))
}
}
if (length(band_files) == 0) {
stop("Could not match any files to provided band names.", call. = FALSE)
}
# Stack the matched bands
band_rasters <- lapply(band_files, terra::rast)
stacked <- terra::rast(band_rasters)
names(stacked) <- names(band_files)
}
return(stacked)
}
#' Load and stack individual band files
#' @keywords internal
load_and_stack_bands <- function(file_list, band_names = NULL, verbose = FALSE) {
# Validate files exist
missing_files <- file_list[!file.exists(file_list)]
if (length(missing_files) > 0) {
stop(sprintf("Files not found: %s", paste(missing_files, collapse = ", ")), call. = FALSE)
}
# Load all rasters
rasters <- lapply(file_list, terra::rast)
# Check compatibility and align if needed
reference_raster <- rasters[[1]]
reference_crs <- terra::crs(reference_raster)
for (i in 2:length(rasters)) {
# Fix CRS if needed
if (!identical(terra::crs(rasters[[i]]), reference_crs)) {
if (verbose) message(sprintf("Reprojecting band %d to match reference CRS", i))
rasters[[i]] <- terra::project(rasters[[i]], reference_crs)
}
# Align geometry if needed
if (!terra::compareGeom(reference_raster, rasters[[i]], stopOnError = FALSE)) {
if (verbose) message(sprintf("Resampling band %d to match reference geometry", i))
rasters[[i]] <- terra::resample(rasters[[i]], reference_raster)
}
}
# Stack the rasters
stacked <- terra::rast(rasters)
# Name the layers
if (!is.null(band_names) && length(band_names) == length(file_list)) {
names(stacked) <- band_names
} else {
# Use filenames without extension
names(stacked) <- tools::file_path_sans_ext(basename(file_list))
}
if (verbose) message(sprintf("Successfully stacked %d bands", terra::nlyr(stacked)))
return(stacked)
}
#' Get available indices
#' @keywords internal
get_available_indices <- function() {
list_vegetation_indices(detailed = FALSE)
}
#' Extract bands from multi-band raster
#' @keywords internal
extract_bands_from_raster <- function(spectral_data, auto_detect, band_names, verbose) {
if (verbose) message("Extracting bands from multi-band raster...")
n_bands <- terra::nlyr(spectral_data)
layer_names <- names(spectral_data)
bands <- list(red = NULL, nir = NULL, blue = NULL, green = NULL,
swir1 = NULL, swir2 = NULL, red_edge = NULL, coastal = NULL, nir2 = NULL)
if (auto_detect || is.null(band_names)) {
# Auto-detect bands based on common naming conventions
bands <- auto_detect_spectral_bands(spectral_data, layer_names, verbose)
} else {
# Use provided band names
bands <- map_custom_band_names(spectral_data, band_names, verbose)
}
# FALLBACK: If auto-detection failed, use direct name matching
if (is.null(bands$red) && "red" %in% layer_names) {
bands$red <- spectral_data[["red"]]
if (verbose) message("Direct match found for red band")
}
if (is.null(bands$green) && "green" %in% layer_names) {
bands$green <- spectral_data[["green"]]
if (verbose) message("Direct match found for green band")
}
if (is.null(bands$blue) && "blue" %in% layer_names) {
bands$blue <- spectral_data[["blue"]]
if (verbose) message("Direct match found for blue band")
}
if (is.null(bands$nir) && "nir" %in% layer_names) {
bands$nir <- spectral_data[["nir"]]
if (verbose) message("Direct match found for nir band")
}
if (is.null(bands$red_edge) && "red_edge" %in% layer_names) {
bands$red_edge <- spectral_data[["red_edge"]]
if (verbose) message("Direct match found for red_edge band")
}
if (is.null(bands$swir1) && "swir1" %in% layer_names) {
bands$swir1 <- spectral_data[["swir1"]]
if (verbose) message("Direct match found for swir1 band")
}
if (is.null(bands$swir2) && "swir2" %in% layer_names) {
bands$swir2 <- spectral_data[["swir2"]]
if (verbose) message("Direct match found for swir2 band")
}
if (is.null(bands$coastal) && "coastal" %in% layer_names) {
bands$coastal <- spectral_data[["coastal"]]
if (verbose) message("Direct match found for coastal band")
}
# Final fallback: use band positions if names don't work
if (is.null(bands$red) && n_bands >= 1) {
bands$red <- spectral_data[[1]]
if (verbose) message("Using band 1 as red (fallback)")
}
if (is.null(bands$nir) && n_bands >= 4) {
bands$nir <- spectral_data[[4]]
if (verbose) message("Using band 4 as nir (fallback)")
}
if (is.null(bands$green) && n_bands >= 2) {
bands$green <- spectral_data[[2]]
if (verbose) message("Using band 2 as green (fallback)")
}
if (is.null(bands$blue) && n_bands >= 3) {
bands$blue <- spectral_data[[3]]
if (verbose) message("Using band 3 as blue (fallback)")
}
return(bands)
}
#' Auto-detect spectral bands
#' @keywords internal
auto_detect_spectral_bands <- function(spectral_data, layer_names, verbose) {
bands <- list(red = NULL, nir = NULL, blue = NULL, green = NULL,
swir1 = NULL, swir2 = NULL, red_edge = NULL, coastal = NULL, nir2 = NULL)
# Comprehensive band naming patterns
patterns <- list(
red = c("red", "RED", "Red", "B4", "b4", "band4", "Band_4", "Band4", "B04", "b04"),
nir = c("nir", "NIR", "Nir", "B8", "b8", "band8", "Band_8", "Band8", "B5", "b5", "near_infrared", "B08", "b08"),
blue = c("blue", "BLUE", "Blue", "B2", "b2", "band2", "Band_2", "Band2", "B02", "b02"),
green = c("green", "GREEN", "Green", "B3", "b3", "band3", "Band_3", "Band3", "B03", "b03"),
swir1 = c("swir1", "SWIR1", "Swir1", "B11", "b11", "B6", "b6", "shortwave_infrared_1", "SWIR_1"),
swir2 = c("swir2", "SWIR2", "Swir2", "B12", "b12", "B7", "b7", "shortwave_infrared_2", "SWIR_2"),
red_edge = c("rededge", "RedEdge", "red_edge", "RED_EDGE", "B5", "b5", "RE", "re", "B05", "b05"),
coastal = c("coastal", "COASTAL", "Coastal", "B1", "b1", "aerosol", "B01", "b01", "COASTAL_AEROSOL")
)
# Try exact matches first
for (band_type in names(patterns)) {
for (pattern in patterns[[band_type]]) {
if (pattern %in% layer_names) {
bands[[band_type]] <- spectral_data[[pattern]]
if (verbose) message(sprintf("Exact match detected %s band: %s", band_type, pattern))
break
}
}
}
# Try partial matches if exact matches failed
for (band_type in names(patterns)) {
if (is.null(bands[[band_type]])) {
for (pattern in patterns[[band_type]]) {
matches <- which(grepl(pattern, layer_names, ignore.case = TRUE))
if (length(matches) > 0) {
bands[[band_type]] <- spectral_data[[matches[1]]]
if (verbose) message(sprintf("Partial match detected %s band: %s", band_type, layer_names[matches[1]]))
break
}
}
}
}
return(bands)
}
#' Load and validate spectral band
#' @keywords internal
load_and_validate_band <- function(band_data, band_name, required = FALSE) {
if (is.null(band_data)) {
if (required) {
stop(sprintf("%s band is required but not provided", band_name), call. = FALSE)
}
return(NULL)
}
if (is.character(band_data)) {
band_raster <- terra::rast(band_data)
} else {
band_raster <- band_data
}
# Basic validation
if (terra::ncell(band_raster) == 0) {
stop(sprintf("%s band has no cells", band_name), call. = FALSE)
}
return(band_raster)
}
#' Validate required bands for specific indices
#' @keywords internal
validate_required_bands <- function(index_type, blue, green, swir1, swir2, red_edge, coastal, nir2, verbose = FALSE) {
requirements <- list(
"EVI" = list(blue = blue),
"GNDVI" = list(green = green),
"ARVI" = list(blue = blue),
"VARI" = list(blue = blue, green = green),
"NDWI" = list(green = green),
"MNDWI" = list(green = green, swir1 = swir1),
"NDMI" = list(swir1 = swir1),
"MSI" = list(swir1 = swir1),
"NDII" = list(swir1 = swir1),
"NDRE" = list(red_edge = red_edge),
"MTCI" = list(red_edge = red_edge),
"IRECI" = list(red_edge = red_edge),
"S2REP" = list(red_edge = red_edge),
"PSRI" = list(red_edge = red_edge),
"NBR" = list(swir2 = swir2),
"BAI" = list(swir1 = swir1),
"NDSI" = list(green = green, swir1 = swir1),
"CCI" = list(red_edge = red_edge),
"LAI" = list(blue = blue),
"FAPAR" = list(blue = blue),
"MTVI1" = list(green = green),
"MTVI2" = list(green = green),
"WI" = list(swir1 = swir1),
"SRWI" = list(swir1 = swir1),
"LSWI" = list(swir1 = swir1),
"CI" = list(green = green),
"GBNDVI" = list(green = green, blue = blue),
"NPCI" = list(blue = blue),
"NPQI" = list(blue = blue)
)
if (index_type %in% names(requirements)) {
required_bands <- requirements[[index_type]]
for (band_name in names(required_bands)) {
band_data <- required_bands[[band_name]]
# Check if band is NULL or invalid
if (is.null(band_data)) {
if (verbose) message(sprintf("%s index requires %s band but it's NULL", index_type, band_name))
return(FALSE)
}
# Check if it's a valid SpatRaster
if (!inherits(band_data, "SpatRaster")) {
if (verbose) message(sprintf("%s index requires %s band but it's not a SpatRaster", index_type, band_name))
return(FALSE)
}
# Check if it has valid data
if (terra::ncell(band_data) == 0) {
if (verbose) message(sprintf("%s index requires %s band but it has no cells", index_type, band_name))
return(FALSE)
}
}
}
return(TRUE)
}
#' Calculate index by type
#' @keywords internal
calculate_index_by_type <- function(index_type, red, nir, blue, green, swir1, swir2,
red_edge, coastal, nir2, verbose) {
index <- switch(index_type,
# ==================== BASIC VEGETATION INDICES ====================
"NDVI" = (nir - red) / (nir + red),
"SAVI" = {
L <- 0.5 # Standard soil brightness correction factor
((nir - red) / (nir + red + L)) * (1 + L)
},
"MSAVI" = {
# Modified SAVI - self-adjusting L parameter
0.5 * (2 * nir + 1 - sqrt((2 * nir + 1)^2 - 8 * (nir - red)))
},
"OSAVI" = {
# Optimized SAVI with L = 0.16
(nir - red) / (nir + red + 0.16)
},
"EVI" = {
# Enhanced Vegetation Index
2.5 * ((nir - red) / (nir + 6 * red - 7.5 * blue + 1))
},
"EVI2" = {
# Two-band Enhanced Vegetation Index
2.5 * ((nir - red) / (nir + 2.4 * red + 1))
},
"DVI" = nir - red,
"RVI" = nir / red,
"GNDVI" = (nir - green) / (nir + green),
"WDVI" = {
# Weighted Difference Vegetation Index
# Using a typical soil line slope of 0.5
nir - 0.5 * red
},
# ==================== ENHANCED VEGETATION INDICES ====================
"ARVI" = {
# Atmospherically Resistant Vegetation Index
rb <- red - 1 * (red - blue) # Self-correction for atmospheric scattering
(nir - rb) / (nir + rb)
},
"RDVI" = {
# Renormalized Difference Vegetation Index
(nir - red) / sqrt(nir + red)
},
"PVI" = {
# Perpendicular Vegetation Index (using standard coefficients)
# Based on soil line: NIR = a*Red + b, with typical values a=1.5, b=10
a <- 1.5; b <- 10
(nir - a * red - b) / sqrt(1 + a^2)
},
"IPVI" = {
# Infrared Percentage Vegetation Index
nir / (nir + red)
},
"TNDVI" = {
# Transformed NDVI
sqrt(((nir - red) / (nir + red)) + 0.5)
},
"GEMI" = {
# Global Environment Monitoring Index
eta <- (2 * (nir^2 - red^2) + 1.5 * nir + 0.5 * red) / (nir + red + 0.5)
eta * (1 - 0.25 * eta) - (red - 0.125) / (1 - red)
},
"VARI" = {
# Visible Atmospherically Resistant Index
(green - red) / (green + red - blue)
},
"TSAVI" = {
# Transformed Soil Adjusted Vegetation Index
# Using standard parameters: a=1.5, b=10, X=0.08
a <- 1.5; b <- 10; X <- 0.08
(a * (nir - a * red - b)) / (red + a * nir - a * b + X * (1 + a^2))
},
"ATSAVI" = {
# Adjusted Transformed Soil Adjusted Vegetation Index
a <- 1.22; b <- 0.03; X <- 0.08
(a * (nir - a * red - b)) / (a * nir + red - a * b + X * (1 + a^2))
},
"GESAVI" = {
# Generalized Soil Adjusted Vegetation Index
Z <- 0.5 # Adjustment factor
((nir - red) / (nir + red + Z)) * (1 + Z)
},
"MTVI" = {
# Modified Triangular Vegetation Index (without red edge)
1.2 * (1.2 * (nir - green) - 2.5 * (red - green))
},
"CTVI" = {
# Corrected Transformed Vegetation Index
ndvi_val <- (nir - red) / (nir + red)
((ndvi_val + 0.5) / abs(ndvi_val + 0.5)) * sqrt(abs(ndvi_val + 0.5))
},
# ==================== RED EDGE INDICES ====================
"NDRE" = (nir - red_edge) / (nir + red_edge),
"MTCI" = {
# MERIS Terrestrial Chlorophyll Index
(red_edge - red) / (nir - red)
},
"IRECI" = {
# Inverted Red-Edge Chlorophyll Index
(red_edge - red) / (red_edge / nir)
},
"S2REP" = {
# Sentinel-2 Red-Edge Position
705 + 35 * ((red + red_edge) / 2 - red) / (red_edge - red)
},
"PSRI" = {
# Plant Senescence Reflectance Index
(red - green) / red_edge
},
"CRI1" = (1 / green) - (1 / red),
"CRI2" = (1 / green) - (1 / red_edge),
"ARI1" = (1 / green) - (1 / red_edge),
"ARI2" = nir * ((1 / green) - (1 / red_edge)),
"MCARI" = {
# Modified Chlorophyll Absorption Ratio Index
((red_edge - red) - 0.2 * (red_edge - green)) * (red_edge / red)
},
# ==================== STRESS DETECTION INDICES ====================
"PRI" = {
# Photochemical Reflectance Index
# Using 531nm (green) as proxy for 531nm reference band
(green - nir) / (green + nir)
},
"SIPI" = {
# Structure Insensitive Pigment Index
(nir - red) / (nir - green)
},
"CCI" = {
# Canopy Chlorophyll Index
(red_edge - red) / (red_edge + red)
},
"NDNI" = {
# Normalized Difference Nitrogen Index
(log(1 / nir) - log(1 / swir1)) / (log(1 / nir) + log(1 / swir1))
},
"CARI" = {
# Chlorophyll Absorption Ratio Index
red_edge * (red / green)
},
"TCARI" = {
# Transformed Chlorophyll Absorption Ratio Index
3 * ((red_edge - red) - 0.2 * (red_edge - green) * (red_edge / red))
},
"MTVI1" = {
# Modified Triangular Vegetation Index 1
1.2 * (1.2 * (nir - green) - 2.5 * (red - green))
},
"MTVI2" = {
# Modified Triangular Vegetation Index 2
1.5 * (1.2 * (nir - green) - 2.5 * (red - green)) /
sqrt((2 * nir + 1)^2 - (6 * nir - 5 * sqrt(red)) - 0.5)
},
"TVI" = {
# Triangular Vegetation Index
0.5 * (120 * (nir - green) - 200 * (red - green))
},
"NPCI" = {
# Normalized Pigment Chlorophyll Index
(red - blue) / (red + blue)
},
"RARS" = {
# Ratio Analysis of Reflectance Spectra
red / nir
},
"NPQI" = {
# Normalized Phaeophytinization Index
(red - blue) / (red + blue)
},
# ==================== WATER/MOISTURE INDICES ====================
"NDWI" = {
# Original NDWI (McFeeters 1996) - Green and NIR for water bodies
(green - nir) / (green + nir)
},
"MNDWI" = {
# Modified NDWI (Xu 2006) - Green and SWIR1 for water bodies
(green - swir1) / (green + swir1)
},
"NDMI" = {
# Normalized Difference Moisture Index (Gao 1996) - NIR and SWIR1 for vegetation moisture
(nir - swir1) / (nir + swir1)
},
"MSI" = {
# Moisture Stress Index
swir1 / nir
},
"NDII" = {
# Normalized Difference Infrared Index
(nir - swir1) / (nir + swir1)
},
"WI" = {
# Water Index
nir / swir1
},
"SRWI" = {
# Simple Ratio Water Index
nir / swir1
},
"LSWI" = {
# Land Surface Water Index
(nir - swir1) / (nir + swir1)
},
# ==================== SPECIALIZED INDICES ====================
"LAI" = {
# Leaf Area Index approximation
3.618 * (2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)) - 0.118
},
"FAPAR" = {
# Fraction of Absorbed PAR
-0.161 + 1.257 * ((nir - red) / (nir + red))
},
"FCOVER" = {
# Fractional Vegetation Cover
ndvi_val <- (nir - red) / (nir + red)
-2.274 + 4.336 * ndvi_val - 1.33 * ndvi_val^2
},
"NBR" = {
# Normalized Burn Ratio
(nir - swir2) / (nir + swir2)
},
"BAI" = {
# Burn Area Index
1 / ((0.1 - red)^2 + (0.06 - nir)^2)
},
"NDSI" = {
# Normalized Difference Snow Index
(green - swir1) / (green + swir1)
},
"GRVI" = {
# Green-Red Vegetation Index
(green - red) / (green + red)
},
"VIG" = {
# Vegetation Index Green
(green - red) / (green + red)
},
"CI" = {
# Coloration Index
(red - green) / red
},
"GBNDVI" = {
# Green-Blue NDVI
(nir - (green + blue)) / (nir + green + blue)
},
# If index not found, throw error
stop(sprintf("Unsupported index type: %s", index_type), call. = FALSE)
)
# Handle special cases for division by zero and invalid operations
index <- handle_index_edge_cases(index, index_type, verbose)
return(index)
}
#' Handle edge cases for index calculations
#' @keywords internal
handle_index_edge_cases <- function(index, index_type, verbose) {
# Handle division by zero cases (ratio indices)
if (index_type %in% c("RVI", "MSI", "WI", "SRWI", "RARS")) {
# For ratio indices, handle zero denominators
is_invalid <- is.infinite(terra::values(index, mat = FALSE)) |
is.nan(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to division by zero in %s", n_invalid, index_type))
}
}
}
# Handle square root of negative numbers
if (index_type %in% c("MSAVI", "MTVI2", "RDVI", "TNDVI", "CTVI")) {
is_invalid <- is.nan(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to square root of negative in %s", n_invalid, index_type))
}
}
}
# Handle logarithm edge cases
if (index_type %in% c("NDNI")) {
is_invalid <- is.infinite(terra::values(index, mat = FALSE)) |
is.nan(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to logarithm issues in %s", n_invalid, index_type))
}
}
}
# Handle reciprocal operations (CRI indices)
if (index_type %in% c("CRI1", "CRI2", "ARI1", "ARI2")) {
is_invalid <- is.infinite(terra::values(index, mat = FALSE))
if (any(is_invalid, na.rm = TRUE)) {
index[is_invalid] <- NA
if (verbose) {
n_invalid <- sum(is_invalid, na.rm = TRUE)
message(sprintf("Set %d pixels to NA due to reciprocal operations in %s", n_invalid, index_type))
}
}
}
# Handle S2REP special case (Red Edge Position)
if (index_type == "S2REP") {
# Constrain to reasonable wavelength range
index[index < 680 | index > 780] <- NA
}
# Handle BAI special case (very large values possible)
if (index_type == "BAI") {
# Cap extremely large values
index[index > 1000] <- 1000
}
return(index)
}
#' Apply quality filter to remove outliers
#' @keywords internal
apply_quality_filter <- function(ndvi, verbose) {
# Calculate statistics
values <- terra::values(ndvi, mat = FALSE)
values <- values[!is.na(values)]
if (length(values) == 0) return(ndvi)
# Use IQR method for outlier detection
q1 <- quantile(values, 0.25, na.rm = TRUE)
q3 <- quantile(values, 0.75, na.rm = TRUE)
iqr <- q3 - q1
# Conservative outlier bounds (1.5 * IQR)
lower_bound <- q1 - 1.5 * iqr
upper_bound <- q3 + 1.5 * iqr
# For NDVI, also apply reasonable physical bounds
lower_bound <- max(lower_bound, -0.5) # NDVI rarely below -0.5 for vegetation
upper_bound <- min(upper_bound, 1.0) # NDVI cannot exceed 1.0
# Mask outliers
outlier_mask <- (ndvi < lower_bound) | (ndvi > upper_bound)
ndvi[outlier_mask] <- NA
if (verbose) {
n_outliers <- sum(terra::values(outlier_mask, mat = FALSE), na.rm = TRUE)
if (n_outliers > 0) {
message(sprintf("Quality filter removed %d outlier pixels (%.1f%%)",
n_outliers, (n_outliers/terra::ncell(ndvi))*100))
}
}
return(ndvi)
}
#' Apply temporal smoothing
#' @keywords internal
apply_temporal_smoothing <- function(ndvi_stack, verbose) {
n_layers <- terra::nlyr(ndvi_stack)
if (n_layers < 3) {
if (verbose) message("Not enough layers for temporal smoothing (need at least 3)")
return(ndvi_stack)
}
if (verbose) message(sprintf("Applying temporal smoothing to %d layers", n_layers))
# Create smoothed layers using moving average
smoothed_layers <- list()
for (i in 1:n_layers) {
if (i == 1) {
# First layer: average with next
smoothed_layers[[i]] <- (ndvi_stack[[i]] + ndvi_stack[[i + 1]]) / 2
} else if (i == n_layers) {
# Last layer: average with previous
smoothed_layers[[i]] <- (ndvi_stack[[i - 1]] + ndvi_stack[[i]]) / 2
} else {
# Middle layers: 3-point moving average
smoothed_layers[[i]] <- (ndvi_stack[[i - 1]] + ndvi_stack[[i]] + ndvi_stack[[i + 1]]) / 3
}
}
# Combine into stack
result <- terra::rast(smoothed_layers)
names(result) <- names(ndvi_stack)
if (verbose) message("Temporal smoothing completed")
return(result)
}
#' Mask invalid values based on index type
#' @keywords internal
mask_invalid_values <- function(index, index_type, verbose) {
# Define reasonable ranges for different indices based on literature
ranges <- list(
# Basic vegetation indices
"NDVI" = c(-1, 1), "SAVI" = c(-1, 1.5), "MSAVI" = c(0, 2), "OSAVI" = c(-1, 1),
"EVI" = c(-1, 3), "EVI2" = c(-1, 3), "DVI" = c(-2, 2), "RVI" = c(0, 30),
"GNDVI" = c(-1, 1), "WDVI" = c(-2, 2),
# Enhanced indices
"ARVI" = c(-1, 1), "RDVI" = c(-2, 2), "PVI" = c(-2, 2), "IPVI" = c(0, 1),
"TNDVI" = c(0, 1.2), "GEMI" = c(-1, 1), "VARI" = c(-1, 1),
"TSAVI" = c(-1, 1), "ATSAVI" = c(-1, 1), "GESAVI" = c(-1, 1.5),
# Red edge indices
"NDRE" = c(-1, 1), "MTCI" = c(0, 8), "IRECI" = c(0, 8), "S2REP" = c(680, 780),
"PSRI" = c(-1, 1), "CRI1" = c(-5, 5), "CRI2" = c(-5, 5),
"ARI1" = c(-5, 5), "ARI2" = c(-50, 50), "MCARI" = c(0, 2),
# Stress indices
"PRI" = c(-1, 1), "SIPI" = c(0, 2), "CCI" = c(-1, 1), "NDNI" = c(-1, 1),
"CARI" = c(0, 10), "TCARI" = c(0, 5), "MTVI1" = c(-2, 2), "MTVI2" = c(0, 8),
"TVI" = c(0, 150), "NPCI" = c(-1, 1), "RARS" = c(0, 5), "NPQI" = c(-1, 1),
# Water/moisture indices
"NDWI" = c(-1, 1), "MNDWI" = c(-1, 1), "NDMI" = c(-1, 1), "MSI" = c(0, 5),
"NDII" = c(-1, 1), "WI" = c(0, 5), "SRWI" = c(0, 5), "LSWI" = c(-1, 1),
# Specialized indices
"LAI" = c(0, 15), "FAPAR" = c(0, 1), "FCOVER" = c(0, 1), "NBR" = c(-1, 1),
"BAI" = c(0, 1000), "NDSI" = c(-1, 1), "GRVI" = c(-1, 1), "VIG" = c(-1, 1),
"CI" = c(-1, 1), "GBNDVI" = c(-1, 1)
)
if (index_type %in% names(ranges)) {
valid_range <- ranges[[index_type]]
invalid_mask <- (index < valid_range[1]) | (index > valid_range[2])
n_invalid <- sum(terra::values(invalid_mask, mat = FALSE), na.rm = TRUE)
if (n_invalid > 0) {
index[invalid_mask] <- NA
if (verbose) {
message(sprintf("Masked %d invalid pixels outside range [%.3f, %.3f] for %s (%.1f%%)",
n_invalid, valid_range[1], valid_range[2], index_type,
(n_invalid/terra::ncell(index))*100))
}
}
} else {
if (verbose) message(sprintf("No predefined range for %s, skipping range validation", index_type))
}
return(index)
}
#' Validate output quality
#' @keywords internal
validate_output <- function(index, index_type, verbose) {
values <- terra::values(index, mat = FALSE)
n_total <- length(values)
n_valid <- sum(!is.na(values))
# Instead of stopping, just warn and continue
if (n_valid == 0) {
warning(sprintf("All %s values are invalid/NA. This index may not be suitable for your data.", index_type),
call. = FALSE)
if (verbose) {
message(sprintf("Warning: %s produced no valid values - continuing with other indices", index_type))
}
return(invisible(FALSE)) # Return FALSE to indicate failure but don't stop
}
# Warn if very low data coverage but continue
coverage_pct <- (n_valid/n_total) * 100
if (coverage_pct < 10) {
warning(sprintf("Only %.1f%% of %s values are valid. This index may not be suitable for your data.",
coverage_pct, index_type), call. = FALSE)
} else if (coverage_pct < 50 && verbose) {
message(sprintf("Warning: Only %.1f%% of %s values are valid", coverage_pct, index_type))
}
if (verbose && n_valid > 0) {
value_range <- range(values, na.rm = TRUE)
message(sprintf("%s calculation completed:", index_type))
message(sprintf(" - Valid pixels: %d/%d (%.1f%%)", n_valid, n_total, coverage_pct))
message(sprintf(" - Value range: [%.6f, %.6f]", value_range[1], value_range[2]))
# Add statistical summary for verbose output
if (n_valid >= 10) {
values_clean <- values[!is.na(values)]
message(sprintf(" - Mean: %.6f, Std Dev: %.6f", mean(values_clean), sd(values_clean)))
message(sprintf(" - Percentiles (25th, 50th, 75th): [%.4f, %.4f, %.4f]",
quantile(values_clean, 0.25), quantile(values_clean, 0.5), quantile(values_clean, 0.75)))
}
}
return(invisible(TRUE)) # Return TRUE to indicate success
}
#' Select appropriate indices for crop analysis
#' @keywords internal
select_crop_indices <- function(crop_type, analysis_type, growth_stage) {
# Base indices that are always reliable with just Red + NIR
base_indices <- c("NDVI", "SAVI", "DVI")
# Crop-specific indices - ONLY USING RED + NIR COMPATIBLE INDICES
crop_indices <- switch(crop_type,
"corn" = c("NDVI", "EVI2", "SAVI", "MSAVI", "DVI", "RVI"),
"soybeans" = c("NDVI", "EVI2", "SAVI", "OSAVI", "DVI"),
"wheat" = c("NDVI", "SAVI", "MSAVI", "DVI"),
"cotton" = c("NDVI", "EVI2", "SAVI", "MSAVI"),
"rice" = c("NDVI", "SAVI", "DVI"),
"general" = c("NDVI", "EVI2", "SAVI", "MSAVI", "OSAVI", "DVI", "RVI"),
base_indices
)
# Analysis-specific additions - ONLY ROBUST RED+NIR INDICES
if (analysis_type %in% c("stress", "comprehensive")) {
# Only add stress indices that work with Red + NIR
# SIPI needs Green band, so we exclude it
# PRI needs different bands, so we exclude it
# Just use NDVI variations for stress detection
stress_indices <- c("RDVI") # Renormalized DVI can indicate stress
crop_indices <- unique(c(crop_indices, stress_indices))
}
if (analysis_type %in% c("yield", "comprehensive")) {
# Only add yield-relevant indices that are stable with Red + NIR
yield_indices <- c("DVI", "RVI")
crop_indices <- unique(c(crop_indices, yield_indices))
}
if (analysis_type %in% c("growth", "comprehensive")) {
# Growth-related indices that work with Red + NIR
growth_indices <- c("DVI", "RVI")
crop_indices <- unique(c(crop_indices, growth_indices))
}
# Remove any indices that definitely need other bands
problematic_indices <- c("SIPI", "PRI", "GNDVI", "EVI", "ARVI", "VARI", "NDWI", "MNDWI",
"NDMI", "MTCI", "PSRI", "CCI", "NDRE", "CARI", "TCARI")
crop_indices <- setdiff(crop_indices, problematic_indices)
# Ensure we always have at least the basic indices
if (length(crop_indices) == 0) {
crop_indices <- base_indices
}
return(unique(crop_indices))
}
#' Detect vegetation stress
#' @keywords internal
detect_vegetation_stress <- function(vegetation_indices, crop_type, indices, verbose) {
stress_results <- list()
# Only analyze stress indices that are actually present
stress_indices <- intersect(indices, c("SIPI", "NDVI", "EVI", "GNDVI"))
for (idx in stress_indices) {
if (idx %in% names(vegetation_indices)) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) == 0) {
stress_results[[idx]] <- list(error = "No valid values for stress analysis")
next
}
# Define stress thresholds based on literature and crop type
thresholds <- switch(idx,
"NDVI" = list(
healthy = c(0.6, 1.0),
moderate_stress = c(0.4, 0.6),
severe_stress = c(0.0, 0.4)
),
"EVI" = list(
healthy = c(0.4, 1.0),
moderate_stress = c(0.2, 0.4),
severe_stress = c(0.0, 0.2)
),
"SIPI" = list(
healthy = c(1.0, 2.0),
moderate_stress = c(0.8, 1.0),
severe_stress = c(0.5, 0.8)
),
"GNDVI" = list(
healthy = c(0.5, 1.0),
moderate_stress = c(0.3, 0.5),
severe_stress = c(0.0, 0.3)
),
# Default thresholds
list(healthy = c(0.5, 1.0), moderate_stress = c(0.3, 0.5), severe_stress = c(0.0, 0.3))
)
# Calculate stress statistics
healthy_pixels <- sum(values >= thresholds$healthy[1] & values <= thresholds$healthy[2])
moderate_stress_pixels <- sum(values >= thresholds$moderate_stress[1] & values <= thresholds$moderate_stress[2])
severe_stress_pixels <- sum(values >= thresholds$severe_stress[1] & values <= thresholds$severe_stress[2])
total_pixels <- length(values)
stress_results[[idx]] <- list(
healthy_percentage = (healthy_pixels / total_pixels) * 100,
moderate_stress_percentage = (moderate_stress_pixels / total_pixels) * 100,
severe_stress_percentage = (severe_stress_pixels / total_pixels) * 100,
mean_value = mean(values),
median_value = median(values),
std_dev = sd(values),
thresholds_used = thresholds,
total_pixels_analyzed = total_pixels
)
if (verbose) {
message(sprintf("Stress analysis for %s:", idx))
message(sprintf(" Healthy: %.1f%%, Moderate stress: %.1f%%, Severe stress: %.1f%%",
stress_results[[idx]]$healthy_percentage,
stress_results[[idx]]$moderate_stress_percentage,
stress_results[[idx]]$severe_stress_percentage))
}
}
}
return(stress_results)
}
#' Analyze growth stage
#' @keywords internal
analyze_growth_stage <- function(vegetation_indices, crop_type, growth_stage, verbose) {
# Extract key indices for growth analysis
growth_indices <- intersect(names(vegetation_indices), c("NDVI", "EVI", "GNDVI", "DVI"))
if (length(growth_indices) == 0) {
return(list(error = "No suitable indices available for growth stage analysis"))
}
growth_analysis <- list()
for (idx in growth_indices) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) == 0) {
growth_analysis[[idx]] <- list(error = "No valid values")
next
}
growth_analysis[[idx]] <- list(
mean = mean(values),
median = median(values),
std_dev = sd(values),
min = min(values),
max = max(values),
range = max(values) - min(values),
percentiles = quantile(values, c(0.1, 0.25, 0.75, 0.9)),
coefficient_of_variation = sd(values) / mean(values),
n_pixels = length(values)
)
}
# Add growth stage classification if NDVI available
if ("NDVI" %in% names(growth_analysis) && !is.null(growth_analysis$NDVI$mean)) {
ndvi_mean <- growth_analysis$NDVI$mean
# Crop-specific growth stage classification
predicted_stage <- switch(crop_type,
"corn" = {
if (ndvi_mean < 0.3) "emergence" else
if (ndvi_mean < 0.6) "vegetative" else
if (ndvi_mean < 0.8) "reproductive" else "maturity"
},
"soybeans" = {
if (ndvi_mean < 0.4) "emergence" else
if (ndvi_mean < 0.65) "vegetative" else
if (ndvi_mean < 0.8) "reproductive" else "maturity"
},
"wheat" = {
if (ndvi_mean < 0.35) "tillering" else
if (ndvi_mean < 0.7) "stem_elongation" else
if (ndvi_mean < 0.8) "grain_filling" else "maturity"
},
# General classification
{
if (ndvi_mean < 0.3) "early" else
if (ndvi_mean < 0.6) "vegetative" else
if (ndvi_mean < 0.8) "reproductive" else "mature"
}
)
# Calculate confidence based on how close to stage boundaries
stage_boundaries <- c(0.3, 0.6, 0.8)
distances_to_boundaries <- abs(ndvi_mean - stage_boundaries)
min_distance <- min(distances_to_boundaries)
stage_confidence <- max(0, 1 - (min_distance / 0.15)) # Confidence decreases as we approach boundaries
growth_analysis$predicted_growth_stage <- predicted_stage
growth_analysis$stage_confidence <- stage_confidence
growth_analysis$crop_type_used = crop_type
if (verbose) {
message(sprintf("Predicted growth stage: %s (confidence: %.2f)", predicted_stage, stage_confidence))
}
}
return(growth_analysis)
}
#' Analyze yield potential
#' @keywords internal
analyze_yield_potential <- function(vegetation_indices, crop_type, indices, verbose) {
# Select yield-relevant indices
yield_indices <- intersect(indices, c("NDVI", "EVI", "GNDVI", "DVI", "RVI"))
if (length(yield_indices) == 0) {
return(list(error = "No suitable indices available for yield analysis"))
}
yield_analysis <- list()
# Calculate composite yield index using available indices
composite_values <- 0
n_indices <- 0
index_contributions <- list()
for (idx in yield_indices) {
if (idx %in% names(vegetation_indices)) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) > 0) {
# Normalize values to 0-1 range for each index
if (idx %in% c("NDVI", "EVI", "GNDVI")) {
# For these indices, higher values generally indicate better yield potential
normalized <- pmax(0, pmin(1, (values - 0.1) / (0.9 - 0.1)))
} else if (idx == "DVI") {
# DVI can be negative, so handle differently
normalized <- pmax(0, pmin(1, (values + 0.5) / (2.0 + 0.5)))
} else if (idx == "RVI") {
# RVI ranges from near 0 to high values
normalized <- pmax(0, pmin(1, (values - 1) / (10 - 1)))
} else {
# Generic normalization
min_val <- quantile(values, 0.05, na.rm = TRUE)
max_val <- quantile(values, 0.95, na.rm = TRUE)
normalized <- pmax(0, pmin(1, (values - min_val) / (max_val - min_val)))
}
mean_normalized <- mean(normalized, na.rm = TRUE)
composite_values <- composite_values + mean_normalized
n_indices <- n_indices + 1
index_contributions[[idx]] <- list(
mean_normalized = mean_normalized,
raw_mean = mean(values, na.rm = TRUE),
raw_std = sd(values, na.rm = TRUE)
)
if (verbose) {
message(sprintf(" %s contribution: %.3f (raw mean: %.3f)", idx, mean_normalized, mean(values, na.rm = TRUE)))
}
}
}
}
if (n_indices > 0) {
composite_index <- composite_values / n_indices
# Crop-specific yield potential classification
yield_potential_class <- switch(crop_type,
"corn" = {
if (composite_index < 0.3) "Low" else
if (composite_index < 0.6) "Medium" else
if (composite_index < 0.8) "High" else "Very High"
},
"soybeans" = {
if (composite_index < 0.35) "Low" else
if (composite_index < 0.65) "Medium" else
if (composite_index < 0.85) "High" else "Very High"
},
"wheat" = {
if (composite_index < 0.3) "Low" else
if (composite_index < 0.6) "Medium" else
if (composite_index < 0.8) "High" else "Very High"
},
# Default classification
{
if (composite_index < 0.3) "Low" else
if (composite_index < 0.6) "Medium" else
if (composite_index < 0.8) "High" else "Very High"
}
)
yield_analysis <- list(
composite_yield_index = composite_index,
yield_potential_class = yield_potential_class,
indices_used = yield_indices,
n_indices_used = n_indices,
index_contributions = index_contributions,
crop_type = crop_type,
classification_confidence = abs(composite_index - 0.5) * 2 # Higher confidence away from middle
)
if (verbose) {
message(sprintf("Composite yield index: %.3f", composite_index))
message(sprintf("Yield potential class: %s", yield_potential_class))
}
} else {
yield_analysis <- list(error = "No valid data for yield analysis")
}
return(yield_analysis)
}
#' Calculate vegetation statistics
#' @keywords internal
calculate_vegetation_statistics <- function(vegetation_indices, indices, verbose) {
if (!inherits(vegetation_indices, "SpatRaster")) {
return(list(error = "Invalid vegetation indices data"))
}
statistics <- list()
for (idx in indices) {
if (idx %in% names(vegetation_indices)) {
values <- terra::values(vegetation_indices[[idx]], mat = FALSE)
values <- values[!is.na(values)]
if (length(values) > 0) {
statistics[[idx]] <- list(
count = length(values),
mean = mean(values),
median = median(values),
std_dev = sd(values),
min = min(values),
max = max(values),
range = max(values) - min(values),
cv = sd(values) / abs(mean(values)), # Coefficient of variation
percentiles = quantile(values, c(0.05, 0.25, 0.75, 0.95)),
coverage_percent = (length(values) / terra::ncell(vegetation_indices)) * 100
)
# Add histogram data for distribution analysis
if (length(values) >= 100) {
hist_data <- hist(values, breaks = 20, plot = FALSE)
statistics[[idx]]$histogram <- list(
breaks = hist_data$breaks,
counts = hist_data$counts,
density = hist_data$density
)
}
if (verbose) {
stats <- statistics[[idx]]
message(sprintf("%s statistics:", idx))
message(sprintf(" Mean: %.4f +/- %.4f", stats$mean, stats$std_dev))
message(sprintf(" Range: [%.4f, %.4f]", stats$min, stats$max))
message(sprintf(" Coverage: %.1f%% (%d pixels)", stats$coverage_percent, stats$count))
}
} else {
statistics[[idx]] <- list(error = "No valid values for statistics calculation")
if (verbose) message(sprintf("%s: No valid values", idx))
}
} else {
if (verbose) message(sprintf("%s: Index not found in data", idx))
}
}
# Add overall summary
valid_indices <- names(statistics)[!sapply(statistics, function(x) "error" %in% names(x))]
if (length(valid_indices) > 0) {
statistics$summary <- list(
total_indices_calculated = length(valid_indices),
indices_with_valid_data = valid_indices,
total_indices_requested = length(indices),
success_rate = (length(valid_indices) / length(indices)) * 100
)
}
return(statistics)
}
#' Get index formulas
#' @keywords internal
get_index_formulas <- function(indices) {
formulas <- setNames(rep("", length(indices)), indices)
formula_db <- list(
# ==================== BASIC VEGETATION INDICES (10) ====================
"NDVI" = "(NIR - Red) / (NIR + Red)",
"SAVI" = "((NIR - Red) / (NIR + Red + L)) * (1 + L), where L=0.5",
"MSAVI" = "0.5 * (2*NIR + 1 - sqrt((2*NIR + 1)^2 - 8*(NIR - Red)))",
"OSAVI" = "(NIR - Red) / (NIR + Red + 0.16)",
"EVI" = "2.5 * ((NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1))",
"EVI2" = "2.5 * ((NIR - Red) / (NIR + 2.4*Red + 1))",
"DVI" = "NIR - Red",
"RVI" = "NIR / Red",
"GNDVI" = "(NIR - Green) / (NIR + Green)",
"WDVI" = "NIR - 0.5 * Red",
# ==================== ENHANCED/IMPROVED INDICES (12) ====================
"ARVI" = "(NIR - (2*Red - Blue)) / (NIR + (2*Red - Blue))",
"RDVI" = "(NIR - Red) / sqrt(NIR + Red)",
"PVI" = "(NIR - a*Red - b) / sqrt(1 + a^2), where a=1.5, b=10",
"IPVI" = "NIR / (NIR + Red)",
"TNDVI" = "sqrt(((NIR - Red) / (NIR + Red)) + 0.5)",
"GEMI" = "eta * (1 - 0.25*eta) - (Red - 0.125) / (1 - Red), where eta = (2*(NIR^2 - Red^2) + 1.5*NIR + 0.5*Red) / (NIR + Red + 0.5)",
"VARI" = "(Green - Red) / (Green + Red - Blue)",
"TSAVI" = "(a * (NIR - a*Red - b)) / (Red + a*NIR - a*b + X*(1 + a^2)), where a=1.5, b=10, X=0.08",
"ATSAVI" = "(a * (NIR - a*Red - b)) / (a*NIR + Red - a*b + X*(1 + a^2)), where a=1.22, b=0.03, X=0.08",
"GESAVI" = "((NIR - Red) / (NIR + Red + Z)) * (1 + Z), where Z=0.5",
"MTVI" = "1.2 * (1.2*(NIR - Green) - 2.5*(Red - Green))",
"CTVI" = "((NDVI + 0.5) / |NDVI + 0.5|) * sqrt(|NDVI + 0.5|), where NDVI = (NIR - Red) / (NIR + Red)",
# ==================== RED EDGE AND ADVANCED INDICES (10) ====================
"NDRE" = "(NIR - RedEdge) / (NIR + RedEdge)",
"MTCI" = "(RedEdge - Red) / (NIR - Red)",
"IRECI" = "(RedEdge - Red) / (RedEdge / NIR)",
"S2REP" = "705 + 35 * ((Red + RedEdge)/2 - Red) / (RedEdge - Red)",
"PSRI" = "(Red - Green) / RedEdge",
"CRI1" = "(1 / Green) - (1 / Red)",
"CRI2" = "(1 / Green) - (1 / RedEdge)",
"ARI1" = "(1 / Green) - (1 / RedEdge)",
"ARI2" = "NIR * ((1 / Green) - (1 / RedEdge))",
"MCARI" = "((RedEdge - Red) - 0.2*(RedEdge - Green)) * (RedEdge / Red)",
# ==================== STRESS AND CHLOROPHYLL INDICES (12) ====================
"PRI" = "(Green - NIR) / (Green + NIR)",
"SIPI" = "(NIR - Red) / (NIR - Green)",
"CCI" = "(RedEdge - Red) / (RedEdge + Red)",
"NDNI" = "(log(1/NIR) - log(1/SWIR1)) / (log(1/NIR) + log(1/SWIR1))",
"CARI" = "RedEdge * (Red / Green)",
"TCARI" = "3 * ((RedEdge - Red) - 0.2*(RedEdge - Green) * (RedEdge / Red))",
"MTVI1" = "1.2 * (1.2*(NIR - Green) - 2.5*(Red - Green))",
"MTVI2" = "1.5 * (1.2*(NIR - Green) - 2.5*(Red - Green)) / sqrt((2*NIR + 1)^2 - (6*NIR - 5*sqrt(Red)) - 0.5)",
"TVI" = "0.5 * (120*(NIR - Green) - 200*(Red - Green))",
"NPCI" = "(Red - Blue) / (Red + Blue)",
"RARS" = "Red / NIR",
"NPQI" = "(Red - Blue) / (Red + Blue)",
# ==================== WATER AND MOISTURE INDICES (8) ====================
"NDWI" = "(Green - NIR) / (Green + NIR)",
"MNDWI" = "(Green - SWIR1) / (Green + SWIR1)",
"NDMI" = "(NIR - SWIR1) / (NIR + SWIR1)",
"MSI" = "SWIR1 / NIR",
"NDII" = "(NIR - SWIR1) / (NIR + SWIR1)",
"WI" = "NIR / SWIR1",
"SRWI" = "NIR / SWIR1",
"LSWI" = "(NIR - SWIR1) / (NIR + SWIR1)",
# ==================== SPECIALIZED APPLICATIONS (10) ====================
"LAI" = "3.618 * EVI - 0.118, where EVI = 2.5 * ((NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1))",
"FAPAR" = "-0.161 + 1.257 * NDVI, where NDVI = (NIR - Red) / (NIR + Red)",
"FCOVER" = "-2.274 + 4.336*NDVI - 1.33*NDVI^2, where NDVI = (NIR - Red) / (NIR + Red)",
"NBR" = "(NIR - SWIR2) / (NIR + SWIR2)",
"BAI" = "1 / ((0.1 - Red)^2 + (0.06 - NIR)^2)",
"NDSI" = "(Green - SWIR1) / (Green + SWIR1)",
"GRVI" = "(Green - Red) / (Green + Red)",
"VIG" = "(Green - Red) / (Green + Red)",
"CI" = "(Red - Green) / Red",
"GBNDVI" = "(NIR - (Green + Blue)) / (NIR + Green + Blue)"
)
for (idx in indices) {
if (idx %in% names(formula_db)) {
formulas[idx] <- formula_db[[idx]]
} else {
formulas[idx] <- "Formula not available"
}
}
return(formulas)
}
#' Get index typical ranges
#' @keywords internal
get_index_ranges <- function(indices) {
ranges <- setNames(rep("", length(indices)), indices)
range_db <- list(
# ==================== BASIC VEGETATION INDICES (10) ====================
"NDVI" = "[-1, 1]",
"SAVI" = "[-1, 1.5]",
"MSAVI" = "[0, 2]",
"OSAVI" = "[-1, 1]",
"EVI" = "[-1, 3]",
"EVI2" = "[-1, 3]",
"DVI" = "[-2, 2]",
"RVI" = "[0, 30]",
"GNDVI" = "[-1, 1]",
"WDVI" = "[-2, 2]",
# ==================== ENHANCED/IMPROVED INDICES (12) ====================
"ARVI" = "[-1, 1]",
"RDVI" = "[-2, 2]",
"PVI" = "[-2, 2]",
"IPVI" = "[0, 1]",
"TNDVI" = "[0, 1.2]",
"GEMI" = "[-1, 1]",
"VARI" = "[-1, 1]",
"TSAVI" = "[-1, 1.5]",
"ATSAVI" = "[-1, 1.5]",
"GESAVI" = "[-1, 1.5]",
"MTVI" = "[-2, 2]",
"CTVI" = "[-1, 1.5]",
# ==================== RED EDGE AND ADVANCED INDICES (10) ====================
"NDRE" = "[-1, 1]",
"MTCI" = "[0, 8]",
"IRECI" = "[-2, 5]",
"S2REP" = "[680, 780]", # Wavelength in nm
"PSRI" = "[-1, 1]",
"CRI1" = "[-10, 10]",
"CRI2" = "[-10, 10]",
"ARI1" = "[-10, 10]",
"ARI2" = "[-10, 10]",
"MCARI" = "[-2, 2]",
# ==================== STRESS AND CHLOROPHYLL INDICES (12) ====================
"PRI" = "[-1, 1]",
"SIPI" = "[0, 2]",
"CCI" = "[-1, 1]",
"NDNI" = "[-1, 1]",
"CARI" = "[0, 5]",
"TCARI" = "[-2, 5]",
"MTVI1" = "[-2, 2]",
"MTVI2" = "[-2, 5]",
"TVI" = "[-100, 100]",
"NPCI" = "[-1, 1]",
"RARS" = "[0, 5]",
"NPQI" = "[-1, 1]",
# ==================== WATER AND MOISTURE INDICES (8) ====================
"NDWI" = "[-1, 1]",
"MNDWI" = "[-1, 1]",
"NDMI" = "[-1, 1]",
"MSI" = "[0, 5]",
"NDII" = "[-1, 1]",
"WI" = "[0, 10]",
"SRWI" = "[0, 10]",
"LSWI" = "[-1, 1]",
# ==================== SPECIALIZED APPLICATIONS (10) ====================
"LAI" = "[0, 15]",
"FAPAR" = "[0, 1]",
"FCOVER" = "[0, 1]",
"NBR" = "[-1, 1]",
"BAI" = "[0, 1000]",
"NDSI" = "[-1, 1]",
"GRVI" = "[-1, 1]",
"VIG" = "[-1, 1]",
"CI" = "[-1, 1]",
"GBNDVI" = "[-1, 1]"
)
for (idx in indices) {
if (idx %in% names(range_db)) {
ranges[idx] <- range_db[[idx]]
} else {
ranges[idx] <- "Variable"
}
}
return(ranges)
}
#' Get index references
#' @keywords internal
get_index_references <- function(indices) {
references <- setNames(rep("", length(indices)), indices)
ref_db <- list(
# ==================== BASIC VEGETATION INDICES (10) ====================
"NDVI" = "Rouse et al. (1974)",
"SAVI" = "Huete (1988)",
"MSAVI" = "Qi et al. (1994)",
"OSAVI" = "Rondeaux et al. (1996)",
"EVI" = "Huete et al. (1997)",
"EVI2" = "Jiang et al. (2008)",
"DVI" = "Richardson & Wiegand (1977)",
"RVI" = "Birth & McVey (1968)",
"GNDVI" = "Gitelson et al. (1996)",
"WDVI" = "Clevers (1988)",
# ==================== ENHANCED/IMPROVED INDICES (12) ====================
"ARVI" = "Kaufman & Tanre (1992)",
"RDVI" = "Roujean & Breon (1995)",
"PVI" = "Richardson & Wiegand (1977)",
"IPVI" = "Crippen (1990)",
"TNDVI" = "Deering et al. (1975)",
"GEMI" = "Pinty & Verstraete (1992)",
"VARI" = "Gitelson et al. (2002)",
"TSAVI" = "Baret & Guyot (1991)",
"ATSAVI" = "Baret et al. (1992)",
"GESAVI" = "Gilabert et al. (2002)",
"MTVI" = "Haboudane et al. (2004)",
"CTVI" = "Perry & Lautenschlager (1984)",
# ==================== RED EDGE AND ADVANCED INDICES (10) ====================
"NDRE" = "Gitelson & Merzlyak (1994)",
"MTCI" = "Dash & Curran (2004)",
"IRECI" = "Frampton et al. (2013)",
"S2REP" = "Frampton et al. (2013)",
"PSRI" = "Merzlyak et al. (1999)",
"CRI1" = "Gitelson et al. (2002)",
"CRI2" = "Gitelson et al. (2002)",
"ARI1" = "Gitelson et al. (2001)",
"ARI2" = "Gitelson et al. (2001)",
"MCARI" = "Daughtry et al. (2000)",
# ==================== STRESS AND CHLOROPHYLL INDICES (12) ====================
"PRI" = "Gamon et al. (1992)",
"SIPI" = "Penuelas et al. (1995)",
"CCI" = "Barnes et al. (2000)",
"NDNI" = "Serrano et al. (2002)",
"CARI" = "Kim et al. (1994)",
"TCARI" = "Haboudane et al. (2002)",
"MTVI1" = "Haboudane et al. (2004)",
"MTVI2" = "Haboudane et al. (2004)",
"TVI" = "Broge & Leblanc (2001)",
"NPCI" = "Penuelas et al. (1994)",
"RARS" = "Chappelle et al. (1992)",
"NPQI" = "Barnes et al. (1992)",
# ==================== WATER AND MOISTURE INDICES (8) ====================
"NDWI" = "McFeeters (1996)",
"MNDWI" = "Xu (2006)",
"NDMI" = "Gao (1996)",
"MSI" = "Hunt & Rock (1989)",
"NDII" = "Hardisky et al. (1983)",
"WI" = "Gao (1996)",
"SRWI" = "Zarco-Tejada et al. (2003)",
"LSWI" = "Xiao et al. (2002)",
# ==================== SPECIALIZED APPLICATIONS (10) ====================
"LAI" = "Baret & Guyot (1991)",
"FAPAR" = "Myneni & Williams (1994)",
"FCOVER" = "Baret et al. (2007)",
"NBR" = "Lopez Garcia & Caselles (1991)",
"BAI" = "Chuvieco et al. (2002)",
"NDSI" = "Hall et al. (1995)",
"GRVI" = "Tucker (1979)",
"VIG" = "Gitelson et al. (2002)",
"CI" = "Escadafal & Huete (1991)",
"GBNDVI" = "Wang et al. (2007)"
)
for (idx in indices) {
if (idx %in% names(ref_db)) {
references[idx] <- ref_db[[idx]]
} else {
references[idx] <- "Various sources"
}
}
return(references)
}
Try the geospatialsuite package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.