Raster Analysis as Easy as Making Toast 🍞

Making Raster Analysis Without Losing Your Mind

code
GIS
analysis
Author

Someone Who’s Been There

Published

June 24, 2025

Terra Package: Making Raster Analysis as Easy as Making Toast 🍞

Or: How I Learned to Stop Worrying and Love Spatial Data

Introduction: Why Terra is Your New Best Friend

Remember the good old days when working with raster data in R felt like trying to assemble IKEA furniture blindfolded? Those days are GONE!

The terra package has swooped in like a superhero to save us from the dark ages of complicated raster operations. Think of Terra as the Swiss Army knife of spatial data - except this knife actually works and doesn’t require a PhD in mechanical engineering to operate.

library(terra)
library(ggplot2)
library(dplyr)

# Set up some fancy options because we're fancy people
options(digits = 3)

Chapter 1: Creating Our First Raster (It’s Easier Than Making Instant Noodles)

Let’s start by creating a raster from scratch. Why? Because we can, and because it’s surprisingly therapeutic.

# Create a simple raster - it's like painting, but with numbers!
my_first_raster <- rast(nrows = 10, ncols = 10, 
                        xmin = 0, xmax = 10, 
                        ymin = 0, ymax = 10)

# Fill it with some random values (because life is random)
values(my_first_raster) <- runif(100, min = 0, max = 100)

# Let's see what we've created
plot(my_first_raster, main = "My Beautiful Random Raster 🎨")

Fun Fact: Creating this raster was easier than explaining to your grandparents why you need another streaming subscription.

Chapter 2: Loading Real Data (Because Random Numbers Are Only Fun for So Long)

Now let’s work with some real data. We’ll create a simple elevation model because mountains are cool and everyone loves mountains.

# Create a more interesting raster - let's make some fake mountains!
rows <- 50
cols <- 50

# Create coordinate matrices
x_coords <- matrix(rep(1:cols, rows), nrow = rows, byrow = TRUE)
y_coords <- matrix(rep(1:rows, cols), nrow = rows)

# Create some "elevation" data using a fancy mathematical formula
# (Don't worry, the math fairy did the heavy lifting)
elevation_data <- 100 * exp(-((x_coords - 25)^2 + (y_coords - 25)^2) / 200) +
                  50 * exp(-((x_coords - 15)^2 + (y_coords - 35)^2) / 150) +
                  runif(rows * cols, -5, 5)

# Create the raster
elevation_raster <- rast(nrows = rows, ncols = cols,
                        xmin = 0, xmax = 100, 
                        ymin = 0, ymax = 100)
values(elevation_raster) <- as.vector(elevation_data)

# Plot it with style
plot(elevation_raster, 
     main = "Mount Awesome and Little Sister Peak 🏔️",
     col = terrain.colors(50))

Chapter 3: Basic Raster Operations (The Fun Stuff!)

Now for the magic! Terra makes raster operations so easy, you’ll think you’re cheating.

# Basic statistics - because numbers tell stories
cat("📊 Raster Statistics:\n")
cat("Min elevation:", global(elevation_raster, "min", na.rm = TRUE)[[1]], "m\n")
cat("Max elevation:", global(elevation_raster, "max", na.rm = TRUE)[[1]], "m\n")
cat("Mean elevation:", global(elevation_raster, "mean", na.rm = TRUE)[[1]], "m\n")

# Create a slope raster (because we're getting fancy now)
slope_raster <- terrain(elevation_raster, "slope")

# Plot both side by side
par(mfrow = c(1, 2))
plot(elevation_raster, main = "Elevation 🗻", col = terrain.colors(50))
plot(slope_raster, main = "Slope 📐", col = heat.colors(50))
par(mfrow = c(1, 1))

Chapter 4: Raster Algebra (It’s Like Regular Algebra, But Cooler)

Here’s where Terra really shines. Raster algebra is as easy as basic math!

# Let's classify our elevation into zones
# Because everything is better when organized into categories

# Create elevation zones
flat_areas <- elevation_raster < 20
hilly_areas <- (elevation_raster >= 20) & (elevation_raster < 60)
mountainous_areas <- elevation_raster >= 60

# Combine them into a classification raster
terrain_classes <- flat_areas * 1 + hilly_areas * 2 + mountainous_areas * 3

# Give it proper names
levels(terrain_classes) <- data.frame(ID = 1:3, 
                                     terrain = c("Flat", "Hilly", "Mountainous"))

# Plot the classification
plot(terrain_classes, 
     main = "Terrain Classification 🗺️",
     col = c("lightgreen", "yellow", "brown"),
     legend = TRUE)

Chapter 5: The Grand Finale - Extracting Values

Let’s create some random points and extract values from our raster. It’s like playing darts, but with coordinates!

# Create some random sampling points
set.seed(42) # For reproducible randomness (an oxymoron, but whatever)
sample_points <- spatSample(elevation_raster, 10, "random", as.points = TRUE)

# Extract elevation values at these points
extracted_values <- extract(elevation_raster, sample_points)

# Create a fancy summary
results <- data.frame(
  Point_ID = 1:10,
  Elevation = round(extracted_values[,2], 1),
  Classification = ifelse(extracted_values[,2] < 20, "Flat ⛳", 
                         ifelse(extracted_values[,2] < 60, "Hilly ⛰️", "Mountain 🏔️"))
)

# Display results
knitr::kable(results, caption = "🎯 Sample Point Results")

# Plot everything together
plot(elevation_raster, main = "Sample Points on Elevation Map 📍", 
     col = terrain.colors(50))
points(sample_points, pch = 19, col = "red", cex = 1.5)
text(sample_points, labels = 1:10, pos = 3, col = "white", font = 2)

Chapter 6: Pro Tips (Because We’re Professionals Now)

Here are some Terra tricks that will make you look like a GIS wizard:

# 1. Quick summary function
summary_stats <- function(raster) {
  cat("🔍 Quick Raster Summary:\n")
  cat("Dimensions:", dim(raster)[1], "x", dim(raster)[2], "cells\n")
  cat("Resolution:", res(raster), "\n")
  cat("Extent:", as.vector(ext(raster)), "\n")
  cat("Value range:", range(values(raster), na.rm = TRUE), "\n")
}

summary_stats(elevation_raster)

# 2. Easy resampling (because sometimes size matters)
small_raster <- aggregate(elevation_raster, fact = 5, fun = "mean")
cat("\n📏 Original size:", ncell(elevation_raster), "cells")
cat("\n📏 Resampled size:", ncell(small_raster), "cells")

# 3. Masking (like Instagram filters for rasters)
high_elevation_mask <- elevation_raster > 50
masked_raster <- mask(elevation_raster, high_elevation_mask)

plot(masked_raster, 
     main = "Only the High Stuff 🎿", 
     col = terrain.colors(50))

Conclusion: You’re Now a Terra Master! 🎓

Congratulations! You’ve just learned that working with raster data doesn’t have to feel like performing brain surgery with oven mitts.

What we’ve accomplished: - ✅ Created rasters from scratch - ✅ Performed basic operations - ✅ Did raster algebra like a boss - ✅ Extracted values like a data detective - ✅ Applied professional techniques

Key takeaways: 1. Terra makes everything easier (seriously, everything) 2. Raster operations are just fancy math 3. Visualization is your friend 4. You’re now 87% cooler than people who don’t use Terra

Next Steps 🚀

Now that you’re a Terra expert, here are some ideas for your next adventure: - Try working with real satellite imagery - Explore temporal raster analysis - Combine Terra with other spatial packages - Build interactive maps with leaflet - Take over the world (optional)

Remember: With great Terra power comes great responsibility. Use your newfound raster superpowers wisely!