End-to-End Multispectral Image Workflow

Introduction

This vignette demonstrates a complete workflow for processing multispectral satellite imagery using GeoIndexR and terra:

  1. Loading imagery
  2. Identifying bands and mapping sensor aliases
  3. Calculating targeted and comprehensive indices
  4. Validating and summarizing distributions
  5. Visualizing results with thematic color maps
  6. Exporting georeferenced GeoTIFF files

Step 1: Loading Imagery & Inspecting Bands

library(GeoIndexR)
library(terra)

# Load synthetic 6-band multispectral raster
img <- get_example_data()
print(img)
#> class       : SpatRaster
#> size        : 10, 10, 6  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 440000, 440100, 5410000, 5410100  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 31N (EPSG:32631)
#> source(s)   : memory
#> names       :     blue,    green,      red,      nir,    swir1,    swir2
#> min values  : 0.020059, 0.041602, 0.020587, 0.010094, 0.005018,  0.00119
#> max values  : 0.198413, 0.235461, 0.347471, 0.817501, 0.545547, 0.444584
names(img)
#> [1] "blue"  "green" "red"   "nir"   "swir1" "swir2"

Step 2: Band Mapping & Sensor Presets

Satellite images often use sensor-specific band notations (such as B02, B03, B04, B08 for Sentinel-2, or B2, B3, B4, B5 for Landsat 8/9).

GeoIndexR supports sensor presets:

# View pre-configured Sentinel-2 bands
band_mapping("sentinel2")
#>    coastal       blue      green        red   rededge1   rededge2   rededge3 
#>      "B01"      "B02"      "B03"      "B04"      "B05"      "B06"      "B07" 
#>        nir       nir2   rededge4 watervapor     cirrus      swir1      swir2 
#>      "B08"      "B8A"      "B8A"      "B09"      "B10"      "B11"      "B12" 
#>       swir 
#>      "B11"

# View pre-configured Landsat 8/9 bands
band_mapping("landsat8")
#>  coastal     blue    green      red      nir    swir1    swir2      pan 
#>     "B1"     "B2"     "B3"     "B4"     "B5"     "B6"     "B7"     "B8" 
#>   cirrus thermal1 thermal2  thermal     swir 
#>     "B9"    "B10"    "B11"    "B10"     "B6"

If your image layers use non-standard names, you can pass a custom named vector:

my_mapping <- c(
  blue  = "blue",
  green = "green",
  red   = "red",
  nir   = "nir",
  swir  = "swir1"
)

Step 3: Computing Multiple Indices

Compute a multi-index suite covering vegetation, water, urban built-up, and soil:

indices <- geo_indices(
  image   = img,
  indices = c("NDVI", "NDWI", "NDBI", "SAVI", "BSI"),
  bands   = my_mapping
)

print(indices)
#> class       : SpatRaster
#> size        : 10, 10, 5  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 440000, 440100, 5410000, 5410100  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 31N (EPSG:32631)
#> source(s)   : memory
#> names       :      NDVI,      NDWI,      NDBI,      SAVI,       BSI
#> min values  : -0.572015, -0.845765, -0.665743, -0.075831, -0.614885
#> max values  :  0.914035,  0.757753,   0.37077,  0.849635,   0.33727
names(indices)
#> [1] "NDVI" "NDWI" "NDBI" "SAVI" "BSI"

Step 4: Summary Statistics

Use index_summary() to inspect the distributions across all computed layers:

stats <- index_summary(indices)
print(stats)
#> 
#> === GeoIndexR Spectral Summary ===
#> 
#>  index     min    max    mean  median     sd     q05     q25    q75    q95
#>   NDVI -0.5720 0.9140  0.2902  0.1186 0.4983 -0.4446 -0.0436 0.8634 0.8951
#>   NDWI -0.8458 0.7578 -0.2117 -0.2913 0.5321 -0.8124 -0.7196 0.3371 0.6810
#>   NDBI -0.6657 0.3708 -0.2529 -0.3872 0.3182 -0.6394 -0.5262 0.0733 0.2541
#>   SAVI -0.0758 0.8496  0.3019  0.0828 0.3705 -0.0623 -0.0214 0.7211 0.8253
#>    BSI -0.6149 0.3373 -0.2416 -0.3799 0.3072 -0.5770 -0.4656 0.1265 0.2634
#>  na_pct total_cells
#>       1         100
#>       1         100
#>       1         100
#>       1         100
#>       1         100

Step 5: Visualization

plot_index() applies scientifically tailored palettes based on index categories:

# Visualizing NDVI
plot_index(indices, index = "NDVI")

Step 6: Exporting with terra::writeRaster

Since all outputs returned by GeoIndexR are standard terra::SpatRaster objects with full CRS and geotransform preservation, saving results to GeoTIFF or COG is simple:

# Export the multilayer raster to GeoTIFF
terra::writeRaster(indices, "computed_indices.tif", overwrite = TRUE)

This ensures complete compatibility with GDAL, QGIS, ArcGIS, Python rasterio, and web mapping services.