4 The terra package for raster and vector data

The terra package (Hijmans 2022) has functions to create, read, manipulate, and write raster and vector data. Raster data is commonly used to represent spatially continuous phenomena by dividing the study region into a grid of equally sized rectangles (referred to as cells or pixels) with the values of the variable of interest. In terra, multilayer raster data is represented with the SpatRaster class. The SpatVector class is used to represent vector data such as points, lines and polygons, and their attributes. In this chapter, we provide examples on how to read, create, and operate with raster and vector data using terra.

4.1 Raster data

The rast() function can be used to create and read raster data. The writeRaster() function allows us to write raster data. For example, here we use rast() to read raster data representing elevation in Luxembourg from a file from the terra package, and assign it to an object r of class SpatRaster (Figure 4.1).

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

FIGURE 4.1: Elevation raster in Luxembourg obtained from terra.

We can also use rast() to create a SpatRaster object r by specifying the number of columns, the number of rows, as well as the minimum and maximum x and y values.

r <- rast(ncol = 10, nrow = 10,
          xmin = -150, xmax = -80, ymin = 20, ymax = 60)
r

class       : SpatRaster 
dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 7, 4  (x, y)
extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 

We can get the size of the raster with several functions.

nrow(r) # number of rows
ncol(r) # number of columns
dim(r) # dimension
ncell(r) # number of cells

We use values() to set and access values of the raster.

values(r) <- 1:ncell(r)

We can also create a multilayer object using c().

r2 <- r * r
s <- c(r, r2)

Layers can be subsetted with [[ ]].

plot(s[[2]]) # layer 2

Many generic functions can be used to operate with rasters as shown in the following examples.

plot(min(s))
plot(r + r + 10)
plot(round(r))
plot(r == 1)

4.2 Vector data

The class SpatVector of terra allows us to represent vector data such as points, lines, and polygons, as well as the attributes that describe these geometries. We can use vect() to read a shapefile, and writeVector() to write a SpatVector to a file. Here, we obtain the map with the divisions of Luxembourg that is in the shapefile file lux.shp of terra.

pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)

We can also use the vect() function to create a SpatVector. For example, here we create a SpatVector that contains longitude and latitude coordinates of point locations, and attributes representing character names and numeric values for each of the points.

# Longitude and latitude values
long <- c(-0.118092, 2.349014, -3.703339, 12.496366)
lat <- c(51.509865, 48.864716, 40.416729, 41.902782)
longlat <- cbind(long, lat)

# CRS
crspoints <- "+proj=longlat +datum=WGS84"

# Attributes for points
d <- data.frame(
place = c("London", "Paris", "Madrid", "Rome"),
value = c(200, 300, 400, 500))

# SpatVector object
pts <- vect(longlat, atts = d, crs = crspoints)

pts
plot(pts)

4.3 Cropping, masking, and aggregating raster data

The terra package provides several functions to work with raster data. Here, we show examples on how to crop, mask, and aggregate raster data by using a raster file representing temperature data. First, we use the worldclim_country() function of the geodata package (Hijmans et al. 2023) to download global temperature data from the WorldClim database. Specifically, we download monthly average temperature in degree Celsius by specifying the country (country = "Spain"), the variable mean temperature (var = "tavg"), the resolution (res = 10), and the path where to download the data to as a temporary file (path = tempdir()). Figure 4.2 shows maps of the monthly average temperature in Spain.

library(terra)
r <- geodata::worldclim_country(country = "Spain", var = "tavg",
                                res = 10, path = tempdir())
plot(r)

FIGURE 4.2: Monthly average temperature in Spain obtained from geodata.

We can average the temperature raster data over months with the mean() function (Figure 4.3).

r <- mean(r)
plot(r)

FIGURE 4.3: Raster representing the average annual temperature in Spain.

We also download the map of Spain with the rnaturalearth package (South 2017), and delete the Canary Islands region (Figure 4.4).

# Map
library(ggplot2)
map <- rnaturalearth::ne_states("Spain", returnclass = "sf")
map <- map[-which(map$region == "Canary Is."), ] # delete region
ggplot(map) + geom_sf()

We obtain the spatial extent of the map with terra:ext(). Then, we can use crop() to remove the part of the raster that is outside the spatial extent (Figure 4.4).

# Cropping
sextent <- terra::ext(map)
r <- terra::crop(r, sextent)
plot(r)

We can use the mask() function to convert all values outside the map to NA (Figure 4.4).

# Masking
r <- terra::mask(r, vect(map))
plot(r)

The aggregate() function of terra can be used to aggregate groups of cells of a raster in order to create a new raster with a lower resolution (i.e., larger cells). The argument fact of aggregate() denotes the aggregation factor expressed as number of cells in each direction (horizontally and vertically), or two integers denoting the horizontal and vertical aggregation factor. Argument fun specifies the function used to aggregate values (e.g., mean). Figure 4.4 shows a low resolution raster of average annual temperatures in Spain obtained using the aggregate() function.

# Aggregating
r <- terra::aggregate(r, fact = 20, fun = "mean", na.rm = TRUE)
plot(r)

FIGURE 4.4: Map of Spain excluding the Canary Islands region (top-left), and cropped raster (top-right), masked raster (bottom-left), and low resolution raster (bottom-right) representing the average annual temperature in the map.

4.4 Extracting raster values at points

Given a raster of class SpatRaster, we can extract the raster values at a set of points with the extract() function of terra. Here, we provide an example of the use of extract() using a raster representing the elevation of Luxembourg, and a vector file with the divisions of Luxembourg from files in the terra package.

library(terra)
# Raster (SpatRaster)
r <- rast(system.file("ex/elev.tif", package = "terra"))
# Polygons (SpatVector)
v <- vect(system.file("ex/lux.shp", package = "terra"))

We use the terra functions centroids() to obtain the centroids of the division polygons, and crds() to obtain their coordinates.

points <- crds(centroids(v))

Figure 4.5 shows the elevation raster, the polygons, and the points of the centroid locations.

plot(r)
plot(v, add = TRUE)
points(points)

FIGURE 4.5: Elevation raster, and division and centroids of polygons in Luxembourg.

Now, we obtain the values of the raster at points using extract() passing as first argument the SpatVector object with the raster, and as second argument a data frame with the points.

# data frame with the coordinates
points <- as.data.frame(points)
valuesatpoints <- extract(r, points)
cbind(points, valuesatpoints)

       x     y ID elevation
1  6.009 50.07  1       444
2  6.127 49.87  2       295
3  5.887 49.80  3       382
4  6.165 49.93  4       404
5  5.915 49.94  5       414
6  6.378 49.79  6       320
7  6.312 49.55  7       193
8  6.346 49.69  8       228
9  5.964 49.64  9       313
10 6.024 49.52 10       282
11 6.168 49.62 11       328
12 6.114 49.76 12       221

4.5 Extracting and averaging raster values

We can also use extract() to obtain the values of raster objects of class SpatRaster within polygons of class SpatVector. By default, cells extracted within each polygon are cells that have centers covered by the polygon. We can set the argument weights = TRUE to get, apart from the cell values, the percentage of each cell covered by the polygon, and this can be used to compute area-weighted averages. The argument fun of extract() can be used to specify a function (e.g., mean) that summarizes the extracted values by polygon.

# Extracted raster cells within each polygon
head(extract(r, v, na.rm = TRUE))

  ID elevation
1  1       547
2  1       485
3  1       497
4  1       515
5  1       515
6  1       515

# Extracted raster cells and percentage of area
# covered within each polygon
head(extract(r, v, na.rm = TRUE, weights = TRUE))

  ID elevation  weight
1  1        NA 0.04545
2  1        NA 0.10909
3  1       529 0.24545
4  1       542 0.46364
5  1       547 0.68182
6  1       535 0.11818

Thus, the average raster values by polygon are obtained with extract() as follows.

# Average raster values by polygon
v$avg <- extract(r, v, mean, na.rm = TRUE)$elevation

The area-weighted average raster values by polygon are obtained with extract() setting weights = TRUE.

# Area-weighted average raster values by polygon (weights = TRUE)
v$weightedavg <- extract(r, v, mean, na.rm = TRUE,
                 weights = TRUE)$elevation

Figure 4.6 shows maps of the elevation averages and area-weighted averages in each of the Luxembourg divisions created with the ggplot2 (Wickham, Chang, et al. 2022) and tidyterra (Hernangomez 2023b) packages.

library(ggplot2)
library(tidyterra)

# Plot average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = avg)) +
  scale_fill_terrain_c()

# Plot area-weighted average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = weightedavg)) +
  scale_fill_terrain_c()

FIGURE 4.6: Average and area-weighted average of elevation values in each of the divisions of Luxembourg.