Raster Analysis in R

March 16, 2026 · 11 min read ·Updated May 29, 2026 ·advanced

spatialrasterterrargisadvancedterrain

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

Loading and inspecting rasters

terra::rast("file.tif") loads a raster. print(r) shows dimensions, resolution, extent, and CRS. plot(r) renders a quick visualization. Multi-band rasters (satellite imagery, time series) return a SpatRaster with multiple layers. nlyr(r) gives the number of layers. r[[2]] extracts the second layer. names(r) shows layer names. summary(r) gives value statistics per layer.

Raster algebra

Raster algebra applies operations element-wise to entire grids. r1 + r2 adds two rasters cell by cell, both must have the same extent and resolution. r > 100 creates a binary raster where cells above 100 are 1 and others are 0. NDVI (Normalized Difference Vegetation Index) is computed as (nir - red) / (nir + red) where nir and red are raster layers. These operations are vectorized in C and run efficiently even on large rasters.

Masking and cropping

crop(r, extent_obj) clips the raster to a bounding box. mask(r, polygon_sf) sets cells outside the polygon to NA. Combine them: crop(mask(r, boundary), boundary) clips and masks in two steps. trim(r) removes NA border rows and columns added by masking. For large rasters, crop before mask to reduce the area that mask must process.

Raster data model

A raster represents continuous spatial data as a grid of cells (pixels), each storing a value. The metadata includes extent (bounding box), resolution (cell size), and CRS. Unlike vector data where geometries can be any shape, raster cells are always rectangular and uniformly sized.

The terra package is the current standard for raster analysis in R. rast() creates or reads raster objects. terra::rast("elevation.tif") reads a GeoTIFF. Multi-band rasters (such as multispectral satellite imagery) load as SpatRaster objects with multiple layers accessible by index or name.

Key raster properties: nrow(r), ncol(r) give cell dimensions; res(r) gives the cell resolution in CRS units; ext(r) returns the spatial extent; crs(r) returns the coordinate reference system. values(r) extracts all cell values as a matrix.

Raster algebra and map calculations

Raster objects support element-wise arithmetic. r1 + r2 adds two rasters cell by cell. r * 2 multiplies all cells by a scalar. sqrt(r) applies a function cell-wise. These operations automatically align rasters with the same extent and resolution. For rasters with different extents, crop() and resample() prepare them first.

The Normalized Difference Vegetation Index (NDVI) illustrates raster algebra: given red and near-infrared bands, ndvi <- (nir - red) / (nir + red) computes NDVI for every cell simultaneously. This is more efficient than looping over cells.

terra::lapp() applies a function to multiple raster layers at once: lapp(c(nir, red), fun = function(x, y) (x - y) / (x + y)). For operations that depend on local neighborhood context (focal operations), focal(r, w = matrix(1, 3, 3), fun = mean) computes a 3x3 moving average.

Raster classification and reclassification

classify(r, rcl) reclassifies raster values. The rcl argument is a matrix with columns [from, to, becomes] for interval-based classification, or [is, becomes] for value-based substitution. Reclassifying elevation into terrain classes (flat, gentle slope, steep) is a common preprocessing step for habitat modeling.

cut(values(r), breaks = c(0, 100, 500, 2000, 9000)) discretizes continuous raster values using R’s standard numeric cut function. Assign the result back with values(r) <- ... to update the raster in place.

Resampling and reprojection

terra::resample(r, template) resamples a raster to match the resolution and extent of a template. The method argument controls interpolation: "bilinear" for continuous data, "near" (nearest neighbor) for categorical data and class maps.

terra::project(r, crs) reprojects a raster to a new CRS. Reprojection involves resampling because cell boundaries rarely align across projections. For analysis requiring both raster and vector layers, project one to match the other rather than projecting both.

Downsampling (reducing resolution) speeds up exploratory analysis. aggregate(r, fact = 10) reduces resolution by a factor of 10, computing mean values over 10x10 cell blocks. Work with aggregated data during development and switch to full resolution for final outputs.

Raster data characteristics

Raster data represents the world as a grid of cells, where each cell has a value. The cell values might represent elevation, temperature, land cover type, or satellite reflectance. The grid has a spatial extent, the geographic area it covers, and a resolution, the size of each cell. Higher resolution means smaller cells and more spatial detail, but also more cells and larger files. The tradeoff between resolution and file size is central to raster data management.

Unlike vector data where geometries represent discrete features, raster data represents continuous fields. Elevation is everywhere; a raster assigns one elevation value per cell rather than one per mountain or valley. Spatial operations on rasters, interpolation, neighborhood statistics, overlay, use this continuous representation. Many natural phenomena, temperature, precipitation, soil properties, are better modeled as rasters than as discrete features.

Coordinate reference systems and alignment

Raster analysis requires that all layers are in the same coordinate reference system and have the same grid alignment. When layers do not align, cells do not correspond to the same geographic locations and overlay operations produce incorrect results. Reprojecting a raster to a new CRS changes the grid geometry, which requires resampling, deciding what value each new cell gets based on the original cells that overlap it.

Resampling method matters. Nearest neighbor resampling assigns each new cell the value of the closest original cell, which preserves discrete categories correctly. Bilinear resampling computes a weighted average of surrounding cells, which is appropriate for continuous data. Using bilinear resampling on categorical data produces nonsensical fractional values. Always choose the resampling method based on whether the data is continuous or categorical.

Map algebra

Map algebra applies mathematical operations to rasters cell by cell. Adding two rasters produces a new raster where each cell is the sum of the corresponding cells in the inputs. Dividing by a constant divides every cell value. Applying a function transforms each cell independently. These operations work on aligned rasters of the same extent and resolution — cells are matched by position, not by spatial proximity.

Neighborhood operations generalize map algebra by computing each cell’s output value from a window of surrounding cells. A focal mean smooths a raster by replacing each cell with the mean of its neighbors. A focal maximum extracts local peaks. The window size and shape determine the spatial scale of the operation. Larger windows smooth more aggressively but reduce spatial detail.

Next steps

Now that you understand raster analysis in r, explore these related topics to deepen your knowledge and apply these techniques in more complex scenarios.