Tutorials Guides Reference Articles Cookbooks About

Working with sf and terra

· 6 min read · Updated March 16, 2026 · intermediate
spatial sf terra r gis intermediate raster vector

The sf and terra packages form the modern foundation for spatial data analysis in R. While the previous tutorial introduced basic concepts, this guide dives deeper into advanced operations that you will use daily when working with real spatial data.

Prerequisites

Make sure you have the required packages installed:

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

Advanced Vector Operations with sf

Transforming Coordinate Reference Systems

Transforming between coordinate systems is a common task. The key is understanding when to use geographic versus projected CRS:

# Load example data
library(spData)
data(us_states)

# Check current CRS
st_crs(us_states)
# Geographic: NAD83 / US (latitude and longitude)

# Transform to a projected CRS (Albers Equal Area for USA)
us_albers <- st_transform(us_states, "ESRI:102003")

# Verify transformation
st_crs(us_albers)
# Projected CRS: USA Contiguous Albers Equal Area

# Now area calculations are in square meters
us_albers <- us_albers %>%
  mutate(area_km2 = as.numeric(st_area(geometry)) / 1e6)

print(us_albers %>% select(NAME, area_km2) %>% head())
# # A tibble: 6 × 2
#   NAME         area_km2
# 1 Alabama       127095.
# 2 Alaska        ...

Buffer Operations

Buffers create polygons around geometries at a specified distance:

# Create points for cities
cities <- tibble(
  name = c("Denver", "Phoenix", "Salt Lake City"),
  lat = c(39.7392, 33.4484, 40.7608),
  lon = c(-104.9903, -112.0740, -111.8910)
) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# Create buffers in meters (requires projected CRS)
cities_albers <- st_transform(cities, "ESRI:102003")

# Buffer: 100 km, 200 km, and 300 km
buffers_100 <- st_buffer(cities_albers, dist = 100000)
buffers_200 <- st_buffer(cities_albers, dist = 200000)
buffers_300 <- st_buffer(cities_albers, dist = 300000)

# Visualize
plot(st_geometry(buffers_300), col = "lightblue", border = NA)
plot(st_geometry(buffers_200), col = "lightgreen", border = NA, add = TRUE)
plot(st_geometry(buffers_100), col = "lightyellow", border = NA, add = TRUE)
plot(st_geometry(cities_albers), pch = 19, col = "red", add = TRUE)

Geometric Operations

The sf package provides many geometric operations:

# Union: combine polygons into one
library(spData)
data(us_states)
us_west <- us_states %>% filter(REGION == "West")
west_union <- st_union(us_west)

plot(west_union, main = "Western United States")

# Intersection: keep only overlapping areas
# Create a bounding box for Colorado region
bbox <- st_bbox(c(xmin = -115, ymin = 32, xmax = -100, ymax = 45), crs = 4326)
bbox_sfc <- st_as_sfc(bbox)

# Intersect with western states
west_crop <- st_intersection(us_west, bbox_sfc)
plot(west_crop$geometry, main = "Western States cropped to region")

# Difference: subtract one geometry from another
denver <- cities_albers %>% filter(name == "Denver")
denver_buffer <- st_buffer(denver, dist = 200000)
us_albers <- st_transform(us_states, "ESRI:102003")
colorado <- us_albers %>% filter(NAME == "Colorado")

outside_denver <- st_difference(denver_buffer, st_union(colorado))
plot(outside_denver, main = "200km around Denver, excluding Colorado")

Working with Geometry Collections

Some operations return geometry collections that need special handling:

# Create lines
line1 <- st_linestring(rbind(c(0, 0), c(1, 1), c(2, 0)))
line2 <- st_linestring(rbind(c(0, 1), c(1, 2), c(2, 1)))

# Union lines
lines <- st_sfc(line1, line2, crs = 4326)
combined <- st_line_merge(st_cast(st_union(lines), "LINESTRING"))

print(combined)
# LINESTRING (0 0, 1 1, 2 0, 2 1, 1 2, 0 1)

# Cast between geometry types
points <- st_sfc(st_point(c(0, 0)), st_point(c(1, 1)), crs = 4326)
multipoint <- st_cast(points, "MULTIPOINT")

# Polygons from points (convex hull)
hull <- st_convex_hull(st_union(points))

Advanced Raster Operations with terra

Creating and Manipulating Rasters

The terra package provides efficient raster operations:

# Create a multi-layer raster
r1 <- rast(nrows = 10, ncols = 10, 
           xmin = 0, xmax = 10, ymin = 0, ymax = 10,
           crs = "EPSG:4326")
r2 <- rast(nrows = 10, ncols = 10, 
           xmin = 0, xmax = 10, ymin = 0, ymax = 10,
           crs = "EPSG:4326")

# Fill with values
values(r1) <- runif(ncell(r1), 0, 100)
values(r2) <- runif(ncell(r2), 0, 50)

# Stack rasters
r_stack <- c(r1, r2)
names(r_stack) <- c("temperature", "precipitation")

print(r_stack)
# class       : SpatRaster 
# dimensions : 10, 10, 2

Raster Algebra

Perform calculations on rasters efficiently:

# Arithmetic operations
r_combined <- r1 + r2
r_scaled <- r1 * 2

# Conditional operations
r_hot <- r1 > 50  # TRUE/FALSE raster

# Focal operations (moving window)
r_smooth <- focal(r1, w = 3, fun = mean)

# Apply function to each cell
r_log <- app(r1, function(x) log(x + 1))

Cropping and Masking

Extract portions of rasters based on vector data:

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

# Create a spatial extent
ext <- ext(700000, 710000, 4800000, 4810000)

# Crop to extent
elev_crop <- crop(elevation, ext)

# Or use an sf object for cropping
library(spData)
data(us_states)
colorado <- us_states %>% filter(NAME == "Colorado")
colorado_sf <- st_geometry(colorado) %>% st_transform(crs(elevation))

# Crop to Colorado
elev_colorado <- crop(elevation, colorado_sf)

# Mask: set values to NA outside polygon
elev_masked <- mask(elev_colorado, colorado_sf)

# Inverse mask: keep only values inside polygon
elev_inside <- mask(elev_colorado, colorado_sf, inverse = TRUE)

Aggregating and Disaggregating

Change raster resolution:

# Aggregate: combine cells (increase cell size, decrease resolution)
elev_agg <- aggregate(elevation, fact = 2, fun = mean)

# Disaggregate: split cells (decrease cell size, increase resolution)
elev_disagg <- disaggregate(elevation, fact = 2)

# Resample: change to match another raster
target <- rast(nrows = 50, ncols = 50, 
               ext(elevation), crs = crs(elevation))
elev_resampled <- resample(elevation, target, method = "bilinear")

Zonal Statistics

Calculate statistics for zones defined by polygons:

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

# Get only Alabama
alabama <- us_states %>% filter(NAME == "Alabama") %>%
  st_geometry() %>%
  st_transform(crs(elevation))

# Crop elevation to Alabama
elev_al <- crop(elevation, alabama)

# Calculate zonal statistics
# Note: terra uses SpatVector, so we convert
al_vect <- as(alabama, "SpatVector")
zonal_stats <- zonal(elev_al, al_vect, fun = c("mean", "min", "max", "sd"))

print(zonal_stats)
#   layer    mean     min    max       sd
# 1     1 267.419 10.1667 572.67 118.0869

Combining sf and terra

The real power comes from combining vector and raster operations:

# Extract raster values at points
library(spData)
data(us_states)

# Get centroid points for states
centroids <- st_centroid(us_states) %>%
  st_transform(4326)

# Sample raster values at these points
elevation <- rast(system.file("ex/elev.tif", package = "terra"))
coords <- st_coordinates(centroids)
elevation_values <- extract(elevation, coords[, c("X", "Y")])

# Add to dataframe
us_states <- us_states %>%
  mutate(elevation = elevation_values$elevation)

# Plot
ggplot(us_states) +
  geom_sf(aes(fill = elevation)) +
  scale_fill_viridis_c(name = "Elevation (m)") +
  labs(title = "State Centroids Elevation") +
  theme_minimal()

# Rasterize vector data
# Create a template raster
template <- rast(nrows = 100, ncols = 100, 
                ext(us_states), crs = st_crs(us_states)$proj4string)

# Rasterize state areas
us_raster <- rasterize(st_geometry(us_states), template, 
                       field = as.numeric(st_area(us_states)) / 1e6,
                       fun = "first")

plot(us_raster, main = "State Areas (km²)")

Performance Tips

Working with large spatial datasets requires optimization:

# Use spatial indices for faster operations
# sf uses spatial indexing automatically in most operations

# For very large data, consider:
# 1. Only load needed columns
us_states_small <- us_states %>% select(NAME, geometry)

# 2. Use spatial indexing explicitly
idx <- st_make_index(us_states, "bbox")

# 3. Parallel processing for large operations
library(future)
plan(multisession)

# 4. Write intermediate results
# st_write(us_states, "us_states.gpkg")

What You Have Learned

This tutorial covered advanced spatial operations:

TopicPackageDescription
CRS transformationsfTransform between geographic and projected CRS
BufferssfCreate polygons around geometries
Geometric operationssfUnion, intersection, difference
Raster algebraterraPerform calculations on rasters
Crop and maskterraExtract portions of raster data
Zonal statisticsterraCalculate stats by zone
Extract valuesterra + sfGet raster values at points
Rasterizeterra + sfConvert vectors to rasters

See Also