---
title: "Getting Started with geospatialsuite"
author: "geospatialsuite Development Team"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with geospatialsuite}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  warning = FALSE,
  message = FALSE,
  eval = TRUE
)
```

# 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
- **Built-in sample data** for learning and testing
- **Flexible data loading** from files, directories, or objects

## Installation

```{r installation, eval=FALSE}
# Install from CRAN
install.packages("geospatialsuite")
```

```{r load-library}
# Load the package
library(geospatialsuite)
library(terra)
```

## Quick Start with Built-in Sample Data

geospatialsuite includes built-in sample data so you can start immediately:

```{r quick-start}
library(geospatialsuite)
library(terra)

# List available sample datasets
list_sample_datasets()

# Load sample raster data
red <- load_sample_data("sample_red.rds")
nir <- load_sample_data("sample_nir.rds")

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

# Visualize
plot(ndvi, main = "NDVI from Sample Data", col = terrain.colors(100))
```

## Core Workflows

### 1. Calculate Vegetation Indices

```{r calculate-indices}
# Load sample spectral bands
red <- load_sample_data("sample_red.rds")
nir <- load_sample_data("sample_nir.rds")
blue <- load_sample_data("sample_blue.rds")

# Calculate NDVI
ndvi <- calculate_vegetation_index(red = red, nir = nir, index_type = "NDVI")
plot(ndvi, main = "NDVI")

# Calculate EVI (requires blue band)
evi <- calculate_vegetation_index(red = red, nir = nir, blue = blue, index_type = "EVI")
plot(evi, main = "EVI")

# Calculate multiple indices at once
indices <- calculate_multiple_indices(
  red = red,
  nir = nir,
  blue = blue,
  indices = c("NDVI", "EVI", "SAVI"),
  output_stack = TRUE
)

# Plot all indices
plot(indices)
```

### 2. Working with Multi-band Data

```{r multiband-data}
# Load multi-band sample raster
multiband <- load_sample_data("sample_multiband.rds")

# Check band names
names(multiband)

# Auto-detect bands and calculate index
ndvi_auto <- calculate_vegetation_index(
  spectral_data = multiband,
  index_type = "NDVI",
  auto_detect_bands = TRUE
)
```

### 3. Spatial Operations

```{r spatial-operations}
# Load sample vector data
points <- load_sample_data("sample_points.rds")
boundary <- load_sample_data("sample_boundary.rds")

# Calculate NDVI
red <- load_sample_data("sample_red.rds")
nir <- load_sample_data("sample_nir.rds")
ndvi <- calculate_vegetation_index(red = red, nir = nir, index_type = "NDVI")

# Extract raster values to points
points_with_values <- universal_spatial_join(
  source_data = points,
  target_data = ndvi,
  method = "extract"
)

# Check result
head(points_with_values)
```

### 4. Data Visualization

```{r visualization}
# Load sample data
points <- load_sample_data("sample_points.rds")

# Quick map (auto-detects everything)
quick_map(points, variable = "ndvi", title = "NDVI at Sample Sites")
```

## Available Sample Datasets

```{r list-datasets}
# List all available datasets
datasets <- list_sample_datasets()
print(datasets[, c("filename", "type", "description")])
```

**Sample datasets include:**

- `sample_red.rds` - Red band SpatRaster (10×10 pixels)
- `sample_nir.rds` - NIR band SpatRaster
- `sample_blue.rds` - Blue band SpatRaster
- `sample_green.rds` - Green band SpatRaster
- `sample_swir1.rds` - SWIR1 band SpatRaster
- `sample_multiband.rds` - Multi-band stack (5 bands)
- `sample_points.rds` - Sample field locations (sf object)
- `sample_boundary.rds` - Study area polygon (sf object)
- `sample_coordinates.csv` - Tabular data with coordinates

---

# Working with Your Own Data

## Loading Raster Files with geospatialsuite

### Single GeoTIFF File

```{r load-single-tif, eval=FALSE}
# Use geospatialsuite's load_raster_data() function
# It provides robust error handling and validation

# Load a single .tif file
my_raster <- load_raster_data("path/to/your/ndvi.tif")

# Result is a list, extract the raster
ndvi_raster <- my_raster[[1]]

# Now use with geospatialsuite functions
# The raster is ready for analysis
summary(ndvi_raster)
```

### Multiple GeoTIFF Files

```{r load-multiple-tif, eval=FALSE}
# Load multiple Landsat bands with geospatialsuite
# Handles validation automatically

landsat_files <- c(
  "LC08_B4_red.tif",
  "LC08_B5_nir.tif",
  "LC08_B3_green.tif"
)

# geospatialsuite loads them with error checking
bands <- load_raster_data(landsat_files, verbose = TRUE)

# Extract individual bands
red_band <- bands[[1]]
nir_band <- bands[[2]]
green_band <- bands[[3]]

# Calculate indices using geospatialsuite
ndvi <- calculate_vegetation_index(
  red = red_band,
  nir = nir_band,
  index_type = "NDVI"
)

gndvi <- calculate_vegetation_index(
  green = green_band,
  nir = nir_band,
  index_type = "GNDVI"
)
```

### Load from Directory

```{r load-directory, eval=FALSE}
# geospatialsuite can load all rasters from a directory
# Perfect for batch processing

# Load all .tif files from Landsat directory
all_bands <- load_raster_data(
  "/path/to/landsat/imagery/",
  pattern = "\\.(tif|tiff)$",
  verbose = TRUE
)

# geospatialsuite finds, loads, and validates all files
# Returns a list of SpatRaster objects ready to use
cat("Loaded", length(all_bands), "raster files\n")
```

## Real-World Landsat Workflow

```{r landsat-workflow, eval=FALSE}
# Complete workflow using geospatialsuite functions

library(geospatialsuite)

# 1. Load Landsat bands using geospatialsuite
landsat_bands <- load_raster_data(
  "landsat/LC08_L2SP_021033_20240715/",
  pattern = "SR_B[2-5].TIF$",
  verbose = TRUE
)

# geospatialsuite loaded them with validation
# Extract bands (scaled values 0-1 after Collection 2 scaling)
blue <- landsat_bands[[1]]   # After scaling
green <- landsat_bands[[2]]
red <- landsat_bands[[3]]
nir <- landsat_bands[[4]]

# 2. Calculate vegetation indices using geospatialsuite
# The package has 60+ pre-programmed indices
veg_indices <- calculate_multiple_indices(
  red = red,
  nir = nir,
  blue = blue,
  green = green,
  indices = c("NDVI", "EVI", "SAVI", "GNDVI"),
  output_stack = TRUE
)

# 3. Visualize using geospatialsuite
quick_map(veg_indices$NDVI, title = "Landsat 8 NDVI - July 15, 2024")
```

## Working with Vector Data

### Loading Shapefiles

```{r load-shapefile, eval=FALSE}
# Load shapefile with sf (standard approach)
library(sf)
field_boundaries <- sf::st_read("data/farm_fields.shp")

# Then use with geospatialsuite's spatial functions
# Calculate NDVI first
ndvi <- calculate_vegetation_index(red = red, nir = nir, index_type = "NDVI")

# Extract NDVI to field boundaries using geospatialsuite
fields_with_ndvi <- universal_spatial_join(
  source_data = field_boundaries,
  target_data = ndvi,
  method = "extract"
)

# geospatialsuite handles CRS mismatches automatically
# Returns field boundaries with extracted NDVI statistics
head(fields_with_ndvi)
```

### Working with GeoPackages

```{r load-geopackage, eval=FALSE}
# Load GeoPackage (modern format, better than shapefile)
farm_data <- sf::st_read("farm_management.gpkg", layer = "fields")
sample_points <- sf::st_read("farm_management.gpkg", layer = "samples")

# Use geospatialsuite's spatial join
samples_with_indices <- universal_spatial_join(
  source_data = sample_points,
  target_data = veg_indices,
  method = "extract",
  buffer_distance = 30  # 30m buffer
)

# geospatialsuite extracted all indices in the stack
# Each index becomes a column in the result
names(samples_with_indices)
```

## Multi-band Raster with Auto-Detection

```{r multiband-auto, eval=FALSE}
# geospatialsuite's auto-detection feature

# Load a stacked multi-band GeoTIFF
multiband_raster <- load_raster_data("sentinel2_stack.tif")[[1]]

# Name the bands (Sentinel-2 example)
names(multiband_raster) <- c("blue", "green", "red", "nir", "swir1")

# Use geospatialsuite's auto-detection
# It finds the right bands automatically!
indices <- calculate_multiple_indices(
  spectral_data = multiband_raster,
  indices = c("NDVI", "EVI", "MNDWI"),
  auto_detect_bands = TRUE,  # This is geospatialsuite's feature!
  output_stack = TRUE
)

# No need to specify which band is which
# geospatialsuite figured it out!
```

## Complete Real-World Example

```{r complete-workflow, eval=FALSE}
# End-to-end agricultural monitoring with geospatialsuite

library(geospatialsuite)
library(sf)

# 1. Load satellite imagery using geospatialsuite
spectral_bands <- load_raster_data(
  "/path/to/satellite/bands/",
  pattern = "B[0-9].tif$",
  verbose = TRUE
)

# 2. Extract bands
red <- spectral_bands[[3]]
nir <- spectral_bands[[4]]
green <- spectral_bands[[2]]

# 3. Calculate indices using geospatialsuite
crop_health <- calculate_multiple_indices(
  red = red,
  nir = nir,
  green = green,
  indices = c("NDVI", "GNDVI", "SAVI"),
  output_stack = TRUE
)

# 4. Load field data
fields <- sf::st_read("farm_data/fields.shp")

# 5. Extract to fields using geospatialsuite
fields_analysis <- universal_spatial_join(
  source_data = fields,
  target_data = crop_health,
  method = "extract"
)

# 6. Visualize using geospatialsuite
quick_map(fields_analysis, 
          variable = "NDVI",
          title = "Field Health Assessment")

# geospatialsuite handled:
# - Loading multiple files
# - Calculating indices
# - Spatial extraction with CRS handling
# - Visualization
```

## Listing Available Indices

```{r list-indices}
# See all vegetation indices geospatialsuite provides
veg_indices <- list_vegetation_indices()
head(veg_indices[, c("Index", "Category", "Description")])

# See water indices
water_indices <- list_water_indices()
head(water_indices)
```

## Getting Help

```{r help, eval=FALSE}
# Package documentation
help(package = "geospatialsuite")

# Function help
?calculate_vegetation_index
?load_raster_data
?universal_spatial_join
?quick_map
```

```{r test-package}
# Test package installation
test_geospatialsuite_package_simple()
```

## Summary

geospatialsuite provides:

**Data Loading:**

- `load_raster_data()` - Load .tif files with validation
- `load_sample_data()` - Access built-in samples

**Analysis:**

- `calculate_vegetation_index()` - 60+ indices
- `calculate_multiple_indices()` - Batch processing
- Auto band detection feature

**Spatial Operations:**

- `universal_spatial_join()` - Extract values, automatic CRS handling

**Visualization:**

- `quick_map()` - One-line mapping
- `create_spatial_map()` - Custom maps

All with robust error handling and simplified workflows!

---

**For more detailed tutorials, see the other vignettes:**

- Vegetation Indices Analysis
- Agricultural Applications
- Spatial Analysis and Integration
- Water Quality Assessment
- Complete Workflows and Case Studies