Getting Started with GeoSpatialSuite

In geospatialsuite: Comprehensive Geospatiotemporal Analysis and Multimodal Integration Toolkit

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  warning = FALSE,
  message = FALSE
)

Introduction

GeoSpatialSuite is a comprehensive R package for geospatial analysis and visualization. It provides universal functions that work with any region, robust error handling, and simplified workflows for complex spatial analysis tasks.

Key Features

• 60+ vegetation indices with automatic band detection
• Universal spatial mapping with one-line functionality
• Reliable terra-based visualization without complex dependencies
• Agricultural applications including CDL crop analysis
• Water quality assessment with multiple indices
• Robust error handling throughout all functions

Installation

# Install from source
# install.packages("geospatialsuite")


# Load the package
library(geospatialsuite)

# For vignette building - adjust as needed
library(terra)
library(sf)
library(ggplot2)
library(dplyr)

Core Dependencies

GeoSpatialSuite works reliably with just the core dependencies:

# Core packages (required)
library(terra)    # Raster operations and plotting
library(sf)       # Vector operations  
library(ggplot2)  # Static mapping
library(dplyr)    # Data manipulation

# Optional packages (for enhanced features)
library(leaflet)      # Interactive mapping
library(viridis)      # Enhanced color schemes
library(RColorBrewer) # Additional color palettes

Quick Start: One-Line Mapping

The fastest way to create a map with GeoSpatialSuite:

# One-line mapping - auto-detects data type and creates appropriate map
quick_map("mydata.shp")
quick_map("satellite_image.tif") 
quick_map(my_spatial_data, interactive = TRUE)

Core Concepts

Data Type Support

GeoSpatialSuite automatically handles various data formats:

Raster Data: - Single files: "data.tif", "image.nc" - Directories: "/path/to/raster/files/" - Multi-band data: Sentinel-2, Landsat, MODIS - SpatRaster objects from terra

Vector Data: - Shapefiles: "boundaries.shp" - CSV with coordinates: "points.csv" - GeoJSON: "features.geojson" - sf objects

Region Boundaries: - US States: "Ohio", "California" - Countries: "Nigeria", "Brazil" - Custom areas: c(xmin, ymin, xmax, ymax) - CONUS: "CONUS"

Universal Region Support

# US States
ohio_boundary <- get_region_boundary("Ohio")
california_boundary <- get_region_boundary("California")

# Countries  
nigeria_boundary <- get_region_boundary("Nigeria")
brazil_boundary <- get_region_boundary("Brazil")

# Continental US
conus_boundary <- get_region_boundary("CONUS")

# Custom bounding box
custom_area <- get_region_boundary(c(-84.5, 39.0, -82.0, 41.0))

# State and county
franklin_county <- get_region_boundary("Ohio:Franklin")

Band Naming and Satellite Support

Automatic Band Detection

geospatialsuite automatically detects bands from major satellite platforms:

| Satellite | Bands Supported | Example Band Names | |-----------|-----------------|-------------------| | Landsat 8/9 | 7 bands | B1-B7 (B4=Red, B5=NIR) | | Sentinel-2 | 10 bands | B01-B12 (B04=Red, B08=NIR, B05=RedEdge) | | MODIS | 7 bands | band1-band7 | | Generic | Any | red, nir, blue, green, swir1, swir2 |

Case-Insensitive Matching

All band names are case-insensitive: - "red" = "RED" = "Red" = "rEd" - "B4"="b4"="B04"="b04"`

# These are all equivalent
ndvi1 <- calculate_vegetation_index(red = red_band, nir = nir_band, index_type = "NDVI")
ndvi2 <- calculate_vegetation_index(Red = red_band, NIR = nir_band, index_type = "NDVI")
ndvi3 <- calculate_vegetation_index(RED = red_band, nir = nir_band, index_type = "NDVI")

Working with Different Satellites

# Landsat 8/9 (auto-detected)
landsat <- rast("LC08_stack.tif")  # Has B1-B7
ndvi <- calculate_vegetation_index(
  spectral_data = landsat,
  index_type = "NDVI",
  auto_detect_bands = TRUE
)

# Sentinel-2 (auto-detected, includes red edge)
sentinel <- rast("S2_stack.tif")  # Has B01-B12
ndre <- calculate_vegetation_index(
  spectral_data = sentinel,
  index_type = "NDRE",  # Red edge index
  auto_detect_bands = TRUE
)

# Custom names (rename first)
custom <- rast("custom_satellite.tif")
names(custom) <- c("red", "nir", "blue", "green")
evi <- calculate_vegetation_index(
  spectral_data = custom,
  index_type = "EVI",
  auto_detect_bands = TRUE
)

For complete documentation, see the full band naming guide in the package documentation.

Basic Workflows

1. Simple Vegetation Analysis

Calculate NDVI for any region:

# Basic NDVI calculation
ndvi <- calculate_vegetation_index(
  red = "red_band.tif",
  nir = "nir_band.tif", 
  index_type = "NDVI"
)

# Enhanced NDVI with quality filtering
ndvi_enhanced <- calculate_ndvi_enhanced(
  red_data = red_raster,
  nir_data = nir_raster,
  quality_filter = TRUE,
  clamp_range = c(-0.2, 1)
)

# Multiple vegetation indices
vegetation_indices <- calculate_multiple_indices(
  red = red_band, 
  nir = nir_band,
  blue = blue_band,
  indices = c("NDVI", "EVI", "SAVI", "GNDVI"),
  output_stack = TRUE
)

2. Water Quality Assessment

Analyze water bodies and quality:

# Water detection indices
ndwi <- calculate_water_index(
  green = green_band,
  nir = nir_band, 
  index_type = "NDWI"
)

# Enhanced water detection  
mndwi <- calculate_water_index(
  green = green_band,
  swir1 = swir1_band,
  index_type = "MNDWI"
)

# Comprehensive water analysis
water_analysis <- analyze_water_bodies(
  green = "green.tif",
  nir = "nir.tif", 
  swir1 = "swir1.tif",
  region_boundary = "Ohio"
)

3. Agricultural Analysis

Monitor crops and agricultural areas:

# Get crop codes
corn_codes <- get_comprehensive_cdl_codes("corn")
grain_codes <- get_comprehensive_cdl_codes("grains")

# Create crop mask
corn_mask <- create_crop_mask(
  cdl_data = "cdl_2023.tif",
  crop_codes = corn_codes,
  region_boundary = "Iowa"
)

# Analyze crop areas
crop_analysis <- analyze_cdl_crops_dynamic(
  cdl_data = "cdl_2023.tif",
  crop_selection = "soybeans", 
  region_boundary = "Ohio",
  analysis_type = "area"
)

# Access results
total_area_hectares <- crop_analysis$total_area_ha
total_area_acres <- total_area_hectares * 2.47105

4. Spatial Data Integration

Extract and join spatial data:

# Extract raster values to points - most common use case
field_data <- universal_spatial_join(
  source_data = "field_sites.csv",      # CSV with coordinates
  target_data = "satellite_image.tif",  # Any raster
  method = "extract",
  buffer_distance = 100,                # 100m buffer around points
  summary_function = "mean"
)

# Zonal statistics for polygons
watershed_stats <- universal_spatial_join(
  source_data = "precipitation.tif",    # Raster data
  target_data = "watersheds.shp",       # Polygon boundaries  
  method = "zonal",
  summary_function = "mean"
)

# Resample raster to different resolution
resampled <- universal_spatial_join(
  source_data = "fine_resolution.tif",
  target_data = "coarse_template.tif",
  method = "resample"
)

Visualization Examples

Static Maps with Terra

# Simple raster plotting
terra::plot(ndvi_raster, main = "NDVI Analysis", 
            col = colorRampPalette(c("brown", "yellow", "green"))(100))

# RGB composite
plot_rgb_raster(
  raster_data = satellite_data,
  r = 4, g = 3, b = 2,           # False color composite
  stretch = "hist",
  title = "False Color Composite"
)

# Custom NDVI visualization
create_spatial_map(
  spatial_data = ndvi_raster,
  color_scheme = "ndvi",
  region_boundary = "Ohio",
  title = "Ohio NDVI Analysis"
)

Interactive Maps

# Interactive point map
interactive_map <- create_interactive_map(
  spatial_data = monitoring_stations,
  fill_variable = "water_quality",
  basemap = "terrain",
  title = "Water Quality Monitoring"
)

# Interactive with custom popup
leaflet_map <- create_spatial_map(
  spatial_data = crop_data,
  fill_variable = "yield",
  interactive = TRUE,
  title = "Crop Yield Distribution"
)

Error Handling and Troubleshooting

Common Issues and Solutions

CRS Mismatches:

# Automatic CRS fixing
result <- calculate_vegetation_index(
  red = red_band,
  nir = nir_band,
  auto_crs_fix = TRUE,  # Automatically handle CRS differences
  verbose = TRUE        # See what's happening
)

Missing Data:

# Handle missing values
interpolated <- spatial_interpolation_comprehensive(
  spatial_data = sparse_data,
  target_variables = "temperature",
  method = "NN",               # Nearest neighbor for gaps
  na_strategy = "nearest"      # Strategy for NAs
)

Large Datasets:

# Process large datasets efficiently
result <- universal_spatial_join(
  source_data = large_points,
  target_data = large_raster,
  method = "extract",
  chunk_size = 500000,    # Smaller chunks for memory
  parallel = FALSE        # Disable parallel for stability
)

Diagnostic Functions

# Test package functionality
test_results <- test_geospatialsuite_package_simple(verbose = TRUE)

# Quick diagnostic
quick_diagnostic()

# Check function availability  
function_status <- test_function_availability(verbose = TRUE)

Data Loading Examples

Loading Various Data Types

# Load raster data
rasters <- load_raster_data("/path/to/raster/files/")
single_raster <- load_raster_data("satellite_image.tif")

# Load with region boundary
ohio_rasters <- load_raster_data(
  "/path/to/files/",
  pattern = "\\.(tif|tiff)$"
)

# Date extraction from filenames
dates <- extract_dates_universal(
  c("ndvi_2023-05-15.tif", "ndvi_2023-06-15.tif")
)

Working with Different Coordinate Systems

# Automatic CRS handling
vegetation_indices <- calculate_multiple_indices(
  spectral_data = "/path/to/bands/",   # Mixed CRS files
  indices = c("NDVI", "EVI", "SAVI"),
  auto_detect_bands = TRUE,
  auto_crs_fix = TRUE,                 # Fix CRS mismatches
  region_boundary = "Michigan"
)

Performance Tips

Memory Management

• Use chunk_size parameter for large datasets
• Process regions separately for continental analysis
• Use scale_factor > 1 to reduce resolution first

Processing Speed

• Disable parallel = FALSE for stability if needed
• Use method = "auto" for optimal processing
• Apply region boundaries early to reduce data size

File Organization

• Organize files by date or region
• Use consistent naming conventions
• Keep related files in same directory

Next Steps

After mastering the basics:

1. Vegetation Analysis - Explore the vegetation indices vignette for comprehensive plant monitoring
2. Agricultural Applications - Learn crop analysis and monitoring workflows
3. Water Quality - Understand water body assessment and quality monitoring
4. Advanced Workflows - Build complex, multi-step analysis pipelines

Getting Help

# Package documentation
help(package = "geospatialsuite")

# Function help
?calculate_vegetation_index
?universal_spatial_join
?quick_map

# List available vegetation indices
list_vegetation_indices(detailed = TRUE)

# List available water indices  
list_water_indices(detailed = TRUE)

# Test package installation
test_geospatialsuite_package_simple()

Example: Complete Basic Workflow

# Complete example: Analyze crop vegetation in Ohio
library(geospatialsuite)

# 1. Calculate vegetation indices
vegetation <- calculate_multiple_indices(
  red = "landsat_red.tif",
  nir = "landsat_nir.tif", 
  blue = "landsat_blue.tif",
  indices = c("NDVI", "EVI", "SAVI"),
  region_boundary = "Ohio",
  output_stack = TRUE
)

# 2. Create crop mask
crop_mask <- create_crop_mask(
  cdl_data = "cdl_ohio_2023.tif",
  crop_codes = get_comprehensive_cdl_codes("corn"),
  region_boundary = "Ohio"
)

# 3. Apply crop mask to vegetation data
crop_vegetation <- terra::mask(vegetation, crop_mask)

# 4. Create visualization
quick_map(crop_vegetation$NDVI, title = "Ohio Corn NDVI")

# 5. Calculate statistics
crop_stats <- terra::global(crop_vegetation, "mean", na.rm = TRUE)
print(crop_stats)

# 6. Save results
terra::writeRaster(crop_vegetation, "ohio_corn_vegetation.tif")

This workflow demonstrates the core GeoSpatialSuite approach: load data, process with region boundaries, visualize, and analyze - all with minimal code and robust error handling.

Acknowledgments

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

geospatialsuite documentation built on Nov. 6, 2025, 1:06 a.m.