Raster Analysis in R

March 16, 2026 · 5 min read · Updated March 16, 2026 · advanced

spatial raster terra r gis advanced terrain

This tutorial covers advanced raster analysis techniques using the terra package. Building on earlier tutorials in this series, you will learn how to perform terrain analysis, calculate zonal statistics, and work with multi-band satellite imagery.

Prerequisites

Install and load the required packages:

install.packages(c("terra", "tidyverse", "spData"))
library(terra)
library(tidyverse)

Understanding Raster Data

Raster data represents spatial information as a grid of cells, where each cell holds a value. In R, the terra package provides efficient handling of raster operations:

# Create a simple raster from scratch
elev <- rast(nrows = 100, ncols = 100,
             xmin = 0, xmax = 100, ymin = 0, ymax = 100,
             crs = "EPSG:4326")

# Fill with values (elevation in meters)
set.seed(42)
values(elev) <- runif(ncell(elev), 0, 1000)

print(elev)
# class       : SpatRaster 
# dimensions  : 100, 100, 1
# resolution  : 1, 1 
# extent      : 0 100 0 100 
# coord ref.  : EPSG:4326 
# source      : memory 
# names       : layer 
# min value   : 1.43 
# max value   : 998.97

Terrain Analysis

Terrain analysis extracts topographic features from elevation data. The terrain() function computes multiple terrain attributes:

# Load elevation example
elev <- rast(system.file("ex/elev.tif", package = "terra"))

# Calculate terrain attributes
terrain_attrs <- terrain(elev, v = c("slope", "aspect", "tpi", "tri", "car"))

# Slope: measure of steepness (in degrees)
plot(terrain_attrs$slope, main = "Slope (degrees)")

# Aspect: direction of steepest slope (in degrees, 0 = north)
plot(terrain_attrs$aspect, main = "Aspect (degrees)")

# TPI: Topographic Position Index
# Positive values = ridges, negative = valleys
plot(terrain_attrs$tpi, main = "TPI (Topographic Position)")

# TRI: Terrain Ruggedness Index
plot(terrain_attrs$tri, main = "TRI (Terrain Ruggedness)")

Calculating Slope and Aspect Manually

For more control, compute slope and aspect from elevation directly:

# Calculate slope using focal operations
slope_fun <- function(x) {
  if(any(is.na(x))) return(NA)
  
  # Simple finite difference method
  # This is a simplified version
  dx <- x[5] - x[1]  # horizontal change
  dy <- x[7] - x[3]  # vertical change
  
  # Calculate slope in degrees
  slope <- atan(sqrt(dx^2 + dy^2)) * 180 / pi
  return(slope)
}

slope_manual <- focal(elev, w = 3, fun = slope_fun)

Multi-Band Raster Operations

Satellite imagery often comes as multi-band rasters. Each band represents different wavelengths:

# Create a multi-band raster (simulating satellite imagery)
bands <- list()
band_names <- c("blue", "green", "red", "nir", "swir")

for(i in 1:5) {
  r <- rast(nrows = 50, ncols = 50,
            xmin = 0, xmax = 50, ymin = 0, ymax = 50,
            crs = "EPSG:4326")
  values(r) <- runif(ncell(r), 0, 1000)
  bands[[i]] <- r
}

satellite <- rast(bands)
names(satellite) <- band_names

print(satellite)
# class       : SpatRaster 
# dimensions  : 50, 50, 5

Band Math

Perform calculations across bands:

# Normalized Difference Vegetation Index (NDVI)
# NDVI = (NIR - Red) / (NIR + Red)
ndvi <- (satellite$nir - satellite$red) / (satellite$nir + satellite$red)

# Plot NDVI
plot(ndvi, main = "NDVI", col = rev(terrain.colors(50)))

# Enhanced Vegetation Index (EVI)
# EVI = 2.5 * (NIR - Red) / (NIR + 6*Red - 7.5*Blue + 1)
evi <- 2.5 * (satellite$nir - satellite$red) / 
       (satellite$nir + 6*satellite$red - 7.5*satellite$blue + 1)

# Calculate normalized difference water index (NDWI)
ndwi <- (satellite$green - satellite$nir) / (satellite$green + satellite$nir)

Zonal Statistics

Zonal statistics calculate summary values for raster cells grouped by zones (polygons):

library(spData)
data(us_states)

# Load elevation raster
elevation <- rast(system.file("ex/elev.tif", package = "terra"))

# Convert us_states to SpatVector
us_states_vect <- vect(us_states)

# Calculate zonal statistics for each state
# Mean elevation by state
state_stats <- zonal(elevation, us_states_vect, 
                     fun = c("mean", "min", "max", "sd", "sum"))

print(head(state_stats))
#   zoneID    mean     min    max       sd        sum
# 1       1 267.419 10.1667 572.67 118.0869 267419.40

# Add state names back
state_stats$NAME <- us_states_vect$NAME

# Find highest and lowest average elevation states
state_stats %>%
  arrange(desc(mean)) %>%
  select(NAME, mean) %>%
  head(5)

Weighted Zonal Statistics

For area-weighted calculations:

# Calculate proportion of each state above 500m elevation
elev_high <- elevation > 500

# For weighted statistics, use exactextractr or similar
# Here using terra approach with polygons
above_500 <- crop(elev_high, us_states_vect)
above_pct <- zonal(above_500, us_states_vect, fun = "mean")

# Convert to percentage
above_pct$above_500_pct <- above_pct$above_500 * 100

Raster Classification

Classify raster values into categories:

# Create elevation classes
elevation <- rast(system.file("ex/elev.tif", package = "terra"))

# Method 1: Reclassify using matrix
# Define: from, to, becomes
reclass_mat <- matrix(c(
  0, 200, 1,    # Low: 0-200m -> class 1
  200, 400, 2,  # Medium: 200-400m -> class 2
  400, 600, 3,  # High: 400-600m -> class 3
  600, Inf, 4   # Very high: >600m -> class 4
), ncol = 3, byrow = TRUE)

elev_class <- reclassify(elevation, reclass_mat)

# Method 2: Using cut function
elev_classes <- cut(elevation, breaks = c(0, 200, 400, 600, Inf),
                    labels = c("low", "medium", "high", "very_high"))

plot(elev_classes, main = "Elevation Classes")

Slope and Aspect for Solar Analysis

Combine slope and aspect for solar radiation modeling:

# Calculate solar radiation variables
terrain <- terrain(elev, v = c("slope", "aspect"))

# Classify aspect into cardinal directions
aspect_class <- cut(terrain$aspect,
                    breaks = c(0, 45, 135, 225, 315, 360),
                    labels = c("N", "E", "S", "W", "N"))

# North-facing slopes (0-45 and 315-360) receive less direct sunlight
# in northern hemisphere
north_facing <- terrain$aspect >= 315 | terrain$aspect <= 45

# Create solar exposure index
solar_index <- terrain$slope * cos(terrain$aspect * pi / 180)

plot(solar_index, main = "Solar Exposure Index")

Handling Large Rasters

For working with large raster datasets:

# Check memory usage
info(elev)

# Use chunked processing for very large files
read_chunked <- function(rast_file, chunk_size = 1000) {
  # Process in chunks to manage memory
  r <- rast(rast_file)
  n <- nlyr(r)
  
  for(i in seq(1, n, by = chunk_size)) {
    # Process chunk
    chunk <- r[[i:min(i + chunk_size - 1, n)]]
    # ... do something with chunk
  }
}

# Write to file instead of keeping in memory
writeRaster(elev, "elevation_processed.tif", overwrite = TRUE)

What You Have Learned

This tutorial covered advanced raster analysis techniques:

Technique	Description
Terrain analysis	Calculate slope, aspect, TPI, TRI from elevation
Multi-band operations	Work with satellite imagery bands
NDVI/EVI calculation	Compute vegetation indices
Zonal statistics	Summarize raster values by polygon zones
Raster classification	Group continuous values into categories
Solar analysis	Combine slope/aspect for radiation modeling
Memory management	Handle large raster datasets