Mapping Luxembourg with terra

Toy Rasters, Real Elevation Data, and Geospatial Visualization

Owen Horton

What Does terra Do?

terra is designed for modern geospatial analysis in R.

It allows users to:

  • Build and manipulate continuous data structures called SpatRasters
  • Work with vector geometries called SpatVectors
  • Calculate spatial summaries
  • Applications include meteorology, geoscience, and civil engineering

Important Functions

Important functions include, but are not limited to…

  • rast() - [saves the grid as a pixellated coordinate system (SpatRaster)]
  • vect() - [saves the grid inside as a system of minimalist lines (SpatVector)]
  • aggregate() - [increases the size of pixels by a multiplicative factor]
  • extract() - [extracts and summarizes data via groups of spatial boundaries]

Building a Toy SpatRaster

# Load terra
library(terra)

# Create an empty raster
toy_raster <- rast(
  ncol = 100,
  nrow = 100,
  xmin = 0,
  xmax = 10,
  ymin = 0,
  ymax = 10 
) 
# Use ncell(), nrows(), hasValues(), dim(), res(), etc to understand 
# your SpatRaster better

plot(toy_raster)

Assigning Raster Values

set.seed(405)

# Fill raster cells with random values
values(toy_raster) <- runif(ncell(toy_raster))

# Plot the resulting SpatRaster
plot(toy_raster, main = "Toy SpatRaster with 10,000 Cells")

Importing Real Geospatial Data

# Premade SpatRaster of elevations at Luxembourg
Lux_R <- rast(system.file("ex/elev.tif", package = "terra"))

# Premade SpatVector of the divisions of Luxembourg
Lux_V <- vect(system.file("ex/lux.shp", package = "terra"))

nrow(Lux_R) # Number of rows: 90
ncol(Lux_R) # Number of columns: 95
dim(Lux_R) # Dimension (ncol, nrow, nlayers): 90, 95, 1
ncell(Lux_R) # Number of cells: 8550
res(Lux_R) # resolution (size of each pixel): 0.008333333

# Understanding the SpatVector
nrow(Lux_V) # Output: 12
ncol(Lux_V) # Output: 6
dim(Lux_V) # Output: 12, 6
ncell(Lux_V) # Output: 72
# Careful, a SpatVector has no resolution, since only one-dimensional lines
# are used to make a SpatVector

Plotting a SpatRaster and SpatVector Together

# Plot the SpatRaster with appropriate titles
plot(Lux_R, 
     main = "Elevation Map of Luxembourg",
     xlab = "Longitude",
     ylab = "Latitude",
     plg = list(title = "Elevation (ft.)", 
                title.srt = -90, 
                title.x = 6.86,
                title.y = 49.8) 
                # Will likely need to adjust the legend label accordingly
    )

# Don't forget the "add = TRUE" part here to stack rasters on top 
# of one another!
plot(Lux_V, add = TRUE)

Aggregating Pixellated Data

# Increase each pixel size by a factor of the two and find the average of the 
# data they encompass to create a new SpatRaster
Lux_R_Aggregated <- aggregate(Lux_R, fact = 2, fun = "mean", na.rm = TRUE)

res(Lux_R_Aggregated) 
# Output: 0.01666667 - The resolution has now increased by a factor of 2! 


# Plot the resulting aggregated SpatRaster/SpatVector
plot(Lux_R_Aggregated, 
     main = "Elevation Map of Luxembourg",
     xlab = "Longitude",
     ylab = "Latitude",
     plg = list(title = "Elevation (ft.)", 
                title.srt = -90, 
                title.x = 6.86,
                title.y = 49.8) 
    )
plot(Lux_V, add = TRUE)

Extracting Data Grouped by Polygon

# extract() requires a SpatRaster, a SpatVector, and a function of some kind
# In this case, we should also add na.rm = TRUE, and specify that we are
# concerned with the elevation variable, not the lat or long values
Lux_V$avg <- extract(Lux_R, Lux_V, mean, na.rm = TRUE)$elevation

# These two libraries pair nicely when creating visualizations for summarized
# geospatial data
library(ggplot2)
library(tidyterra)

# Plot the resulting extracted map using ggplot and tidyterra
ggplot(data = Lux_V) + geom_spatvector(aes(fill = avg)) + 
  scale_fill_terrain_c() + labs(title = "Average Elevations by Division", x = "Longitude", y = "Latitude", fill = "Elev.")

Benefits and Limitations of terra

  • Smoothly renders continuous spatial data in clean plots (rast() and vect())
  • Can perform spatial data analysis (aggregate(), extract(), and more!)
  • Computational time gets expensive as you decrease resolution
  • Can be tedious to manually implement data (Consider rnaturalearth)

Thank You!