Chapter 6: Raster-vector interactions

Key Learnings from, and Solutions to the exercises in Chapter 6 of the book Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow.

Geocomputation with R
Textbook Solutions
Author

Aditya Dahiya

Published

December 27, 2024

6.1 Introduction

  • Focuses on interactions between raster and vector data models, first introduced in Chapter 2.
  • Covers three main techniques:
    • Raster cropping and masking using vector objects (Section 6.2).
    • Extracting raster values based on vector data (Section 6.3).
    • Raster-vector conversions, explained in Sections 6.4 and 6.5.
Code 
library(sf)        # Simple Features in R
library(terra)     # Handling rasters in R
library(tidyterra) # For plotting rasters in ggplot2
library(magrittr)  # Using pipes with raster objects
library(tidyverse) # All things tidy; Data Wrangling
library(spData)    # Spatial Datasets
library(patchwork) # Composing plots

6.2 Raster Cropping

  • Raster cropping and masking unify spatial extents of data, reduce memory use, and optimize computational resources, especially for remote sensing rasters and vector boundaries integration.

  • Applications:

    • Crop rasters to fit an area of interest.
    • Mask values outside specified bounds.
    • Essential preprocessing for creating maps and analyses.

  • Projection Alignment:

    • Rasters and vectors must share the same projection.
    • Use st_transform() from sf for re-projection.

  • Typical Workflow: shown in Figure 1

    1. Crop raster to the area of interest: terra::crop().
    2. Mask values outside the area: terra::mask().
    3. Combined operation ensures raster extent fits the area and external values are replaced with NA.

  • terra::crop()

    • Purpose: Reduces the extent of a raster object to match a defined area of interest.

    • Key Arguments:

      • x: The raster object to be cropped (e.g., a SpatRaster).
      • y: The spatial object (e.g., SpatRaster, sf, or extent) defining the cropping area.
      • snap: Adjusts alignment of the cropped raster to match the target (options: "near", "out", "in").

    • Output: A cropped raster limited to the extent defined by y.

  • terra::mask()

    • Purpose: Sets raster cell values outside the defined spatial object to NA or a specified value.

    • Key Arguments:

      • x: The raster object to be masked.
      • mask: The spatial object (e.g., SpatRaster, sf, or extent) defining the masking bounds.
      • inverse: If TRUE, masks values inside the bounds instead of outside.
      • updatevalue: Sets a custom value (e.g., 0) for masked areas instead of NA.

    • Output: A raster with specified cells masked based on the spatial bounds of the second argument.

Aspect crop() |> mask() mask()
Extent Modification Reduces raster extent to vector’s bounding box Retains original raster extent
Output Size Smaller raster (optimized for memory) Larger raster (original extent retained)
Use Case When reducing raster size is desired When full raster extent is needed but irrelevant cells must be excluded
Code 
sysfonts::font_add_google("Saira Condensed", "body_font")
showtext::showtext_auto()

theme_set(
  theme_minimal(
    base_family = "body_font",
    base_size = 20
  )
)

theme_custom <- function(...){
  theme(
    panel.grid = element_line(
      linewidth = 0.3,
      linetype = 2
    ),
    plot.margin = margin(0,1,0,1, "mm"),
    axis.text = element_text(
      size = 8
    ),
    axis.ticks = element_blank(),
    axis.ticks.length = unit(0, "mm"),
    legend.title = element_blank(),
    legend.key.width = unit(2, "mm"),
    legend.text = element_text(
      margin = margin(0,0,0,0.5, "mm")
    ),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box.margin = margin(0,0,0,0,"mm")
  )
}

srtm <- rast(system.file("raster/srtm.tif", package = "spDataLarge"))
zion <- read_sf(system.file("vector/zion.gpkg", 
                            package = "spDataLarge")) |> 
  st_transform(st_crs(srtm))

g1 <- ggplot() +
  geom_spatraster(data = srtm) +
  geom_sf(data = zion, fill = NA) +
  labs(title = "Base Map: Raster and Vector") +
  scale_fill_whitebox_c() +
  coord_sf(
    xlim = c(-113.25, -112.85),
    ylim = c(37.14, 37.51),
    default_crs = 4326
  ) +
  theme_custom()

g2 <- ggplot() +
  geom_spatraster(data = srtm |> crop(zion)) +
  geom_sf(data = zion, fill = NA) +
  labs(title = "terra::crop()") +
  scale_fill_whitebox_c() +
  coord_sf(
    xlim = c(-113.25, -112.85),
    ylim = c(37.14, 37.51),
    default_crs = 4326
  ) +
  theme_custom()

g3 <- ggplot() +
  geom_spatraster(data = srtm |> crop(zion) |> mask(zion)) +
  geom_sf(data = zion, fill = NA) +
  labs(title = "crop() |> mask()") +
  scale_fill_whitebox_c() +
  coord_sf(
    xlim = c(-113.25, -112.85),
    ylim = c(37.14, 37.51),
    default_crs = 4326
  ) +
  theme_custom()

g4 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(data = zion, fill = NA) +
  labs(title = "mask()") +
  scale_fill_whitebox_c() +
  coord_sf(
    xlim = c(-113.25, -112.85),
    ylim = c(37.14, 37.51),
    default_crs = 4326
  ) +
  theme_custom()

g5 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion, inverse = TRUE)) +
  geom_sf(data = zion, fill = NA) +
  labs(title = "mask(inverse = TRUE)") +
  scale_fill_whitebox_c() +
  coord_sf(
    xlim = c(-113.25, -112.85),
    ylim = c(37.14, 37.51),
    default_crs = 4326
  ) +
  theme_custom()

g <- g1 + g2 + g3 + g4 + g5 + 
  patchwork::plot_layout(
    nrow = 1,
    guides = "collect"
  ) +
  patchwork::plot_annotation(
    title = "Cropping and masking rasters with vectors",
    theme = theme(
      plot.title = element_text(
        family = "body_font",
        size = 36, 
        lineheight = 0.3,
        hjust = 0.5,
        face = "bold"
      )
    )
  )

ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_2_1.png"),
  plot = g,
  height = 800,
  width = 2200,
  units = "px",
  bg = "white"
)
Figure 1: The workflows and outputs when cropping and masking rasters with vectors.

6.3 Raster Extraction

  • Raster extraction retrieves raster values at specified locations using a geographic selector (typically vector objects like points, lines, or polygons). Uses the terra::extract() function (terra package).

  • Point Selector: Extracts raster cell values for specified points. Example: Extracting elevation values for 20 sample points in Zion National Park using st_sample() as shown in Figure 2

Code 
sysfonts::font_add_google("Saira Condensed", "body_font")
showtext::showtext_auto()
theme_set(
  theme_minimal(
    base_family = "body_font"
  )
)
# Elevation Raster for Zion National Park and nearby areas
srtm <- rast(system.file("raster/srtm.tif", package = "spDataLarge"))

# Zion National Park Boundaries, converted into CRS of `srtm`
zion <- read_sf(
  system.file(
    "vector/zion.gpkg", 
    package = "spDataLarge"
    )
  ) |> 
  st_transform(st_crs(srtm))
# ----------------------------------------------------------------
# Basic data display
g1 <- ggplot() +
  geom_spatraster(data = srtm) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title = "Raw Data",
    subtitle = "Elevation Raster and Zion National Park boundaries"
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,1,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,1,0, "mm"),
      hjust = 0.5,
      size = 14
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

# ----------------------------------------------------------------
# Mask to display only Zion National Park, and
# Display the elevations of randomly selected 20 points

# Get 20 random points inside Zion National Park
set.seed(21)
random_points <- st_sample(zion, size = 20) |> 
  st_sf()

g2 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = random_points,
    size = 1,
    alpha = 0.5,
    pch = 19
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title = "Sampling 20 Random Points",
    subtitle = "Getting rndom poitns from sf::st_sample()"
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,1,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,1,0, "mm"),
      hjust = 0.5,
      size = 14
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

# --------------------------------------------------------------
# Extract heights for these 20 points
random_points <- random_points |> 
  mutate(
    elevation = terra::extract(srtm, random_points) |> pull(srtm)
  )

g3 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = random_points,
    mapping = aes(colour = elevation),
    stroke = 0.5,
    size = 1,
    alpha = 1,
    pch = 1
  ) +
  geom_sf_text(
    data = random_points,
    mapping = aes(label = elevation),
    size = 4,
    alpha = 1,
    family = "body_font",
    hjust = 0,
    vjust = 1,
    nudge_y = -0.002
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  paletteer::scale_colour_paletteer_c("ggthemes::Red") +
  coord_sf(expand = FALSE) +
  guides(colour = "none") +
  labs(
    title = "Extracting elevation on points from rasters",
    subtitle = "Using terra::extract() to get elevations on 20 points",
    x = NULL, y = NULL
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,1,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,1,0, "mm"),
      hjust = 0.5,
      size = 14
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

g <- g1 + g2 + g3 +
  plot_layout(
    nrow = 1,
    guides = "collect"
  )

ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_3_1.png"),
  plot = g,
  height = 600,
  width = 1500,
  units = "px",
  bg = "white"
)

rm(g, g1, g2, g3)
Figure 2: Using {sf} and {terra} to get 20 random points in Zion National Park, and displaying their elevation (above sea level, in metres).
  • Line Selector: Extracts raster values for each cell touched by a line. Recommendation: Split lines into points (using st_segmentize() and st_cast() from sf package) for accurate extraction along transects, such as creating elevation profiles for hiking routes, as shown in Figure 3
Code 
# Get Latitude and Longitude of these random points
df1 <- random_points |>
  st_as_sfc() |> 
  st_as_sf() |>
  rename(geometry = x) |> 
  cbind(st_coordinates(random_points))

# Tibble of northern, sourthern, eastern and western-most points
df2 <- bind_rows(
  df1[which.min(df1$X), ] |> mutate(dir = "west"),
  df1[which.max(df1$X), ] |> mutate(dir = "east"),
  df1[which.min(df1$Y), ] |> mutate(dir = "south"),
  df1[which.max(df1$Y), ] |> mutate(dir = "north"),
)

# Get two paths: east to west, north to south
paths <- bind_rows(
  df2 |> 
    filter(dir %in% c("west", "east")) |> 
    st_union() |> 
    st_cast("LINESTRING") |> 
    st_as_sf() |> 
    mutate(path = "West to East"),

  df2 |> 
    filter(dir %in% c("south", "north")) |> 
    st_union() |> 
    st_cast("LINESTRING") |> 
    st_as_sf() |> 
    mutate(path = "South to North")
)
rm(df1, df2)

# Mid-point Check: Whether the CRS'es match
st_crs(srtm) == st_crs(paths)

# Plot of the paths
g1 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = paths,
    mapping = aes(colour = path),
    linewidth = 0.5,
    arrow = arrow(length = unit(2, "mm"))
  ) +
  geom_sf(
    data = random_points,
    size = 0.5
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  paletteer::scale_colour_paletteer_d(
    "nbapalettes::cavaliers_retro"
  ) +
  coord_sf(expand = FALSE) +
  guides(colour = "none") +
  labs(
    title = "Selecting 2 Paths",
    subtitle = "sing st_cast() and st_coordinates()",
    x = NULL, y = NULL
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,0,0,0, "mm")
  )

# A {sf} object of various waypoints along the paths
# to be able to extract elevation
path_points <- paths |> 
  # Break into line segments of 100 metres each
  st_segmentize(dfMaxLength = 100) |> 
  
  # Convert each line segment into a central point
  # to be able to extract its elevation
  st_cast("POINT") |> 
  
  # Group by path to be able to calculate distance
  # along the path
  group_by(path) |> 
  
  # Extract only the first column (i.e. the distance from
  # first point alone). Otherwise, st_distance() 
  # returns a distance matrix As expected, distance grows 
  # by approx. less than dfMaxLength amount per points
  mutate(dist = st_distance(x)[, 1]) |> 
  ungroup()

# Extract the actual elevation using terra::extract
path_points <- path_points |> 
  mutate(
    elevation = terra::extract(srtm, path_points) |> 
      pull(srtm),
    dist = as.numeric(dist)
  )

g2 <- path_points |> 
  ggplot(
    mapping = aes(
      x = dist,
      y = elevation,
      colour = path
    )
  ) +
  geom_point(
    size = 0.1
  ) +
  geom_line(
    linewidth = 0.2
  ) +
  facet_wrap(~path, nrow = 1) +
  scale_x_continuous(
    labels = scales::label_number(
      scale = 0.001
    )
  ) +
  scale_y_continuous(
    labels = scales::label_number(
      big.mark = ","
    )
  ) +
  paletteer::scale_colour_paletteer_d(
    "nbapalettes::cavaliers_retro"
  ) +
  labs(
    x = "Distance along the path (km)",
    y = "Elevation (m)",
    title = "Elevation Profile along two paths",
    subtitle = "Using st_segmentize() & terra::extract()"
  ) +
  guides(colour = "none") +
  theme(
    axis.text = element_text(
      size = 12,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,0,0,0, "mm"),
    axis.title.y = element_text(
      margin = margin(0,0,0,5, "mm")
    )
  )

my_design = "ABB"

g <- wrap_plots(g1, g2) +
  plot_layout(design = my_design, guides = "collect") 

ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_3_2.png"),
  plot = g,
  height = 600,
  width = 1500,
  units = "px",
  bg = "white"
)

rm(g, g1, g2, paths, path_points)
Figure 3: Creating elevation profiles along routes using st_segmentize() and terra::extract() from Northern most to Southernmost point; and from Westernmost to Easternmost point (amongst the 20 random points generated).
  • Polygon Selector: Extracts multiple raster values per polygon. Example: Summarizing elevation statistics (min, mean, max) for Zion National Park polygons using group_by() and summarize() from dplyr as shown in Figure 4
Code 
# Creating circles around the 20 random points
# Each circle of 1.5 km radius
circles <- random_points |> 
  st_buffer(dist = 1500) |> 
  mutate(ID = row_number())


g1 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = circles, 
    fill = alpha("white", 0.2)
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  coord_sf(expand = FALSE) +
  guides(colour = "none") +
  labs(
    title = "Circles of 3 km diameter",
    subtitle = "Around 20 random points with st_buffer()",
    x = NULL, y = NULL
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

# A custom function for rounded off mean
mean_r <- function(x){round(mean(x), 0)}

# Getting minimum, mean and maximum height for each area
circles <- circles |> 
  
  # Add minimum maximum and mean heights for each circle
  left_join(
    terra::extract(srtm, circles) |> 
      as_tibble() |> 
      # Grouping by each polygon / circle
      group_by(ID) |> 
      summarise(across(srtm, list(min = min, 
                                  mean = mean_r, 
                                  max = max)))
  )

g2 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = circles, 
    fill = alpha("white", 0.5)
  ) +
  geom_sf_text(
    data = circles,
    mapping = aes(
      label = paste0(
        srtm_max, "\n", srtm_min
      )
    ),
    lineheight = 0.3,
    family = "body_font",
    size = 1.5,
    check_overlap = TRUE
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  coord_sf(expand = FALSE) +
  guides(colour = "none") +
  labs(
    title = "Maximum and Minimum elevations for each circle",
    subtitle = "Using st_buffer(), terra::extract() and dplyr::summarize()",
    x = NULL, y = NULL
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

g3 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = circles, 
    fill = alpha("white", 0.5)
  ) +
  geom_sf_text(
    data = circles,
    mapping = aes(
      label = paste0(scales::number(srtm_mean, big.mark = ","), " m")
    ),
    lineheight = 0.3,
    family = "body_font",
    size = 3,
    check_overlap = TRUE
  ) +
  scale_fill_wiki_c(
    name = "Elevation (m)",
    labels = scales::label_number(big.mark = ",")
  ) +
  coord_sf(expand = FALSE) +
  guides(colour = "none") +
  labs(
    title = "Average elevation for each circle",
    subtitle = "Using summarise()",
    x = NULL, y = NULL
  ) +
  theme(
    legend.key.width = unit(2, "mm"),
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

g <- g1 + g2 + g3 +
  plot_layout(
    guides = "collect"
  )

ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_3_3.png"),
  plot = g,
  height = 600,
  width = 1500,
  units = "px",
  bg = "white"
)
Figure 4: Computing average elevation, minimum and maximum elevations within specified areas, like polygons, or in this case circles of diameter 3 km each around the 20 random points.
  • Categorical Raster Extraction: Counts occurrences of raster categories (e.g., land cover) within polygons. Example: Analyzing land cover types in Zion National Park using terra::extract() as shown in Figure 5
Code 
nlcd <- rast(system.file("raster/nlcd.tif", package = "spDataLarge"))

# Get Zion Park Boundary transformed into CRS of NLCD Raster
zion <- zion |> st_transform(crs = st_crs(nlcd))

# Get polygons (circles) transformed into CRS of NLCD Raster
circles1 <- circles |> 
  select(ID) |> 
  st_transform(crs = st_crs(nlcd))

cat_circles <- terra::extract(nlcd, circles1) |> 
  as_tibble() |> 
  group_by(ID) |> 
  count(levels) |> 
  mutate(perc = n / sum(n))

g1 <- ggplot() +
  geom_spatraster(data = nlcd |> mask(zion)) +
  geom_sf(
    data = zion, 
    linewidth = 0.4, 
    fill = "transparent"
  ) +
  geom_sf(
    data = circles, 
    fill = alpha("white", 0.2)
  ) +
  geom_sf_text(
    data = circles, 
    mapping = aes(label = ID),
    size = 5,
    family = "body_font"
  ) +
  paletteer::scale_fill_paletteer_d("MoMAColors::VanGogh") +
  coord_sf(expand = FALSE) +
  guides(colour = "none", fill = "none") +
  labs(
    title = "Circles of 3 km diameter",
    subtitle = "Around 20 random points with st_buffer()",
    x = NULL, y = NULL
  ) +
  theme(
    axis.text = element_text(
      size = 8,
      margin = margin(0,0,0,0, "mm")
    ),
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm")
  )

g4 <- cat_circles |> 
  ggplot(aes(x = 0, y = n, fill = levels)) +
  geom_col(
    position = position_stack(),
    colour = "white",
    linewidth = 0.05
  ) +
  facet_wrap(~ ID) +
  coord_polar(theta = "y") +
  theme_void(
    base_family = "body_font"
  ) +
  labs(
    fill = NULL,
    title = "Land Cover Pie charts for 20 zones",
    subtitle = "Each zone of diameter 3 km, identified by an ID"
  ) +
  paletteer::scale_fill_paletteer_d("MoMAColors::VanGogh") +
  theme(
    plot.title = element_text(
      margin = margin(0,0,2,0, "mm"),
      hjust = 0.5,
      size = 18
    ),
    plot.subtitle = element_text(
      margin = margin(0,0,0,0, "mm"),
      hjust = 0.5,
      size = 14
    ),
    axis.ticks.length = unit(0, "mm"),
    plot.margin = margin(0,2,0,2, "mm"),
    legend.text = element_text(
      margin = margin(0,0,0,1, "mm"),
      size = 14
    ),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box.margin = margin(0,0,0,0, "mm"),
    strip.text = element_text(
      size = 18,
      margin = margin(0,0,0,0, "mm")
    )
  )

g <- g1 + g4 +
  plot_layout(
    design = "ABB"
  )

ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_3_4.png"),
  plot = g,
  height = 600,
  width = 1500,
  units = "px",
  bg = "white"
)
Figure 5: Computing pie charts of different land use types in circles of diameter 3 km around the 20 random points, and using st_transform() for rasters and vectors while performing raster extraction.

Combining the work done in this section to produce a composite graphic shown in Figure 6 for publishing a code demonstration

Code 
library(magick)

i1 <- image_read(here::here("book_solutions", "images",
                            "chapter6_3_1.png"))
i2 <- image_read(here::here("book_solutions", "images",
                            "chapter6_3_2.png"))
i3 <- image_read(here::here("book_solutions", "images",
                            "chapter6_3_3.png"))
i4 <- image_read(here::here("book_solutions", "images",
                            "chapter6_3_4.png"))


c_imgs <- c(i1, i2, i3, i4)

image_append(
  image_scale(c_imgs, "1500x"),
  stack = TRUE
) |> 
  image_write(
    path = here::here("book_solutions", "images",
                            "chapter6_3_5.png")
  )
Figure 6: A composite image of all the work done: extracting raster values from vector points with {terra} and {sf} for points (locations), linestrings (along journey paths) and within polygons (areas) and compiling results with {magick}

  • Further details: —

    • exactextractr::exact_extract() (exactextractr package):

      • Faster than terra::extract() for large datasets.
      • Computes precise polygon overlap fractions for weighted statistics.

    • terra::extract(exact = TRUE) adds a fraction column for detailed overlap calculations. Useful for weighted means or precise coverage estimates but computationally intensive.

6.4 Rasterization

  • Rasterization converts vector objects into raster format, often for analysis (e.g., terrain) or modeling. It simplifies data by standardizing spatial resolution, serving as geographic data aggregation.

  • Key Function: rasterize() in the {terra} package handles rasterization.

Note

terra::rasterize(): Converts vector data into raster format by transferring values from vector geometries to raster cells.

Key Arguments:

  1. x: Input vector data: A SpatVector or a two-column matrix (coordinates for points). Defines the geometries to be rasterized.
  2. y: Template raster: A SpatRaster that determines the extent, resolution, and CRS of the output raster.
  3. field: Specifies the values to transfer to raster:
    • As a character: Variable name from x (e.g., a column in the vector dataset).
    • As a numeric: Direct values or indices recycled to match nrow(x).
  4. values (for matrix input x): Numeric values used for rasterization, recycled if fewer than nrow(x). Can also be a matrix or data frame.
  5. fun: Summarizing function for overlapping geometries in a cell:
    • For lines and polygons: “min”, “max”, “mean”, “count”, “sum”.
    • For points: Any function returning a single value (e.g., mean, length, min). The fun = "length" means a count of the number of points in that cell.
  6. background: Default cell value for areas not covered by any features in x. Default: NA.
  7. touches: Logical: If TRUE, includes all cells touched by a feature (lines/polygons), not just those where centroids fall.
  8. update: Logical: If TRUE, updates existing values in the template raster.
  9. cover: Logical: For polygons, returns the fraction of a cell covered by a feature. Uses sub-cell presence/absence checks.
  10. by: Groups vector data into layers:
    • For SpatVector: Column name or index.
    • For matrices: Vector defining group membership.
  • Default Behaviour of rasterize():
    • If no fun is specified, the last value is used in case of overlaps (fun = "last").
    • When field is empty, raster cells reflect presence/absence of geometries.

Common Use Cases:

  1. Presence/absence rasters: Determine where features exist.
  2. Summary statistics: Aggregate multiple values in a single cell (e.g., sum, mean).
  3. Fractional coverage: Estimate the proportion of a raster cell covered by polygons.

  • Resolution Considerations:

    • Low resolution (large cells): Misses geographic variability.
    • High resolution (small cells): Increases computation time.
    • Resolution choice depends on intended use or alignment with other rasters.

  • Flexibility: Rasterization varies by:

    1. Template raster properties (extent, resolution, CRS).
    2. Input vector type (points, polygons).
    3. Function arguments (e.g., fun, field).

  • Points Rasterization: can be done in three further methods depicted in Figure 7: —

    1. Presence/Absence Raster:
      • Raster indicates whether cells contain cycle hire points.
    2. Count Raster:
      • Uses fun = "length" to count cycle hire points per grid cell.
    3. Summary Statistic Raster:
      • Calculates sum of a variable (e.g., capacity) using field and fun = sum.
Code 
sysfonts::font_add_google("Saira Condensed", "body_font")
showtext::showtext_auto()

theme_set(
  theme_minimal(
    base_family = "body_font",
    base_size = 12
  )
)

theme_custom <- function(...){
  theme(
    legend.key.height = unit(1, "mm"),
    legend.key.width = unit(6, "mm"),
    plot.title = element_text(
      size = 24,
      margin = margin(0,0,2,0, "mm")
    ),
    plot.subtitle = element_text(
      size = 18,
      margin = margin(0,0,1,0, "mm")
    ),
    axis.ticks.length = unit(0, "mm"),
    axis.text = element_text(
      margin = margin(0,0,0,0, "mm")
    ),
    panel.grid = element_line(linewidth = 0.2),
    legend.position = "bottom",
    legend.title.position = "top",
    legend.text = element_text(margin = margin(1,0,0,0, "mm")),
    legend.title = element_text(margin = margin(2,0,1,0, "mm")),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box.margin = margin(0,0,0,0, "mm")
  )
}

# Get cycle hiring points data for London
cycles <- spData::cycle_hire_osm |> 
  # Transform into CRS of EPSG:27700 - British National Grid
  st_transform("EPSG:27700")

# Create a template raster for using in rasterization
template_raster <- rast(
  # Spat Extent of the Vector data
  terra::ext(cycles),
  # Required resolution of the raster
  resolution = 500,
  # CRS for the raster
  crs = crs(cycles)
)

g1 <- ggplot() +
  geom_sf(
    data = cycles,
    mapping = aes(
      fill = capacity
    ),
    pch = 21,
    size = 0.5,
    linewidth = 0.001
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Red",
    na.value = "transparent"
  ) +
  labs(
    title = "Cycle Hire Points (London)",
    subtitle = "OSM Vector data, plotted with {sf}"
  ) +
  theme_custom()

g2 <- ggplot() +
  geom_spatraster(
    data = cycles |> terra::rasterize(template_raster)
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Red",
    na.value = "transparent",
    name = "Hiring points present?"
  ) +
  labs(
    title = "Presence / Absence Raster",
    subtitle = "Default rasterize(); i.e. (fun = \"last\")"
  ) +
  theme_custom()

g3 <- ggplot() +
  geom_spatraster(
    data = cycles |> terra::rasterize(template_raster,
                                      fun = "length")
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Red",
    na.value = "transparent",
    breaks = seq(0, 10, 2),
    name = "Number of hiring points"
  ) +
  labs(
    title = "Count raster (vector points per cell)",
    subtitle = "rasterize(fun = \"length\")"
  ) +
  theme_custom()

g4 <- ggplot() +
  geom_spatraster(
    data = cycles |> terra::rasterize(
      template_raster,
      field = "capacity",
      fun = "sum")
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Red",
    na.value = "transparent",
    name = "Cycles Capacity"
  ) +
  labs(
    title = "Summary Statistic raster",
    subtitle = "rasterize(field = \"capacity\", fun = \"length\")"
  ) +
  theme_custom()

g <- g1 + g2 + g3 + g4 +
  plot_layout(nrow = 1) +
  plot_annotation(
    tag_levels = "a",
    title = "Vectorizing rasters into points - 3 types",
    theme = theme(
      plot.title = element_text(
      size = 48,
      margin = margin(2,0,2,0, "mm"),
      hjust = 0.5
      )
    )
  ) &
  theme(
    plot.tag = element_text(
      size = 24,
      margin = margin(0,0,-5,0, "mm"),
      face = "bold"
    )
  )


ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_4_1.png"),
  plot = g,
  height = 700,
  width = 2200,
  units = "px",
  bg = "white"
)
Figure 7: Rasterizing points geometry of {sf} with terra::rasterize() with three methods. (a) Base vector object, (b) Presence-Absence Raster, (c) Count Raster - number of objects per cell, and, (d) Summary Statistic Raster - a summary statistic calculated for each cell, using a field (of sf object) and function.

  • Line and Polygon Rasterization is shown in Figure 8. The polygon rasterization can proceed in two ways, depending on the argument touches = TRUE / FALSE as shown in Figure 8 .

    • touches = TRUE: Includes all cells touched by lines/polygons.
    • Default (touches = FALSE): Selects cells where centroids fall inside polygons.
Code 
sysfonts::font_add_google("Saira Condensed", "body_font")
showtext::showtext_auto()

theme_set(
  theme_minimal(
    base_family = "body_font",
    base_size = 12
  )
)

theme_custom <- function(...){
  theme(
    legend.key.height = unit(1, "mm"),
    legend.key.width = unit(6, "mm"),
    plot.title = element_text(
      size = 24,
      margin = margin(0,0,2,0, "mm"),
      lineheight = 0.3
    ),
    plot.subtitle = element_text(
      size = 18,
      margin = margin(0,0,1,0, "mm"),
      lineheight = 0.3
    ),
    axis.ticks.length = unit(0, "mm"),
    axis.text = element_text(
      margin = margin(0,0,0,0, "mm")
    ),
    panel.grid = element_line(linewidth = 0.2),
    legend.position = "none"
  )
}

#----------------------------------------------------------------
# Get vector data on State of California from `us_states` data
california <- us_states |> 
  filter(NAME == "California")

california_border <- california |> 
  st_cast("MULTILINESTRING")

template_raster <- rast(
  ext(california),
  crs = st_crs(california)$wkt, # Use crs(california) or st_crs(california)$wkt
  resolution = 0.25  # Give resolution in degrees
)

crs(california)
st_crs(california)$wkt

g1 <- ggplot() +
  geom_sf(data = california) +
  labs(title = "California Map (vector)\nplotted with {sf}") +
  theme_custom()

g2 <- ggplot() +
  geom_spatraster(
    data = california_border |> rasterize(template_raster)
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Classic Red",
    na.value = "transparent"
  ) +
  labs(title = "Multi-Linestring\nwith terra::rasterize()",
       subtitle = "Cells touched by Linestring are coloured") +
  theme_custom()

g3 <- ggplot() +
  geom_spatraster(
    data = california |> rasterize(template_raster,
                                   touches = FALSE)
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Classic Red",
    na.value = "transparent"
  ) +
  labs(title = "Multipolygon\nterra::rasterize(touches = FALSE)",
       subtitle = "Only colours raster cells whose\ncentroids are inside multi-polygon") +
  theme_custom()

g4 <- ggplot() +
  geom_spatraster(
    data = california |> rasterize(template_raster,
                                   touches = FALSE)
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Classic Red",
    na.value = "transparent"
  ) +
  labs(title = "Multipolygon\nterra::rasterize(touches = TRUE)",
       subtitle = "Colours all raster cells which\nare touched by the multi-polygon") +
  theme_custom()

rast1 <- california |> rasterize(
  template_raster,
  touches = FALSE
  )

rast2 <- california_border |> 
  rasterize(
    (rast1 * 0.5),
    update = TRUE
  )

g5 <- ggplot() +
  geom_spatraster(
    data = rast2
  ) +
  paletteer::scale_fill_paletteer_c(
    "ggthemes::Classic Red",
    na.value = "transparent"
  ) +
  labs(title = "Polygon & Linestring\nterra::rasterize(update = TRUE)",
       subtitle = "Rasterized polygon updated with a rasterized linestring") +
  theme_custom()


g <- g1 + g2 + g3 + g4 + g5 +
  plot_layout(nrow = 1, guides = "collect") +
  plot_annotation(
    tag_levels = "a",
    title = "Techniques for Rasterizing lines and polygons with {terra}",
    theme = theme(
      plot.title = element_text(
      size = 48,
      margin = margin(2,0,2,0, "mm"),
      hjust = 0.5
      )
    )
  ) &
  theme(
    plot.tag = element_text(
      size = 24,
      margin = margin(0,0,-5,4, "mm"),
      face = "bold"
    ),
    plot.margin = margin(0,0,0,0, "mm")
  )


ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_4_2.png"),
  plot = g,
  height = 700,
  width = 2200,
  units = "px",
  bg = "white"
)
Figure 8: Rasterizing Multi-Linestrings and Multi-Polygons from {sf} using terra::rasterize(). The figure show: (a) The base vector map of California, (b) Rasterization of California’s borders (LINESTRING) generated from map using sf::st_cast(), (c) A multi-polygon (map of California) rasterized using terra::rasterize(touches = FALSE) colours only the cells whose centroids fall within the multi-polygon, (d) Using terra::rasterize(touches = TRUE) colours all cells which are touched by the multi-polygon, and, (e) Using terra::rasterize(update = TRUE) to compile and update a previous raster of multi-polygon (values halved to maintain colour levels) with a border multi-linestring.

6.5 Spatial Vectorization

  • Converts raster data (spatially continuous) into vector data (spatially discrete: points, lines, polygons). Opposite process of rasterization. It can be of three types – (1) points, (2) lines, and (3) polygons.
  • (1) Points from Raster Centroids:
    • Use terra::as.points() (Hijmans 2024) for converting raster cells (non-NA) to points. It produces a SpatVector class of object, which can then be converted into an {sf} (Pebesma and Bivand 2023) object using sf::st_as_sf() .
    • Example: Elevation raster visualized as points in Figure 9.
Code 
theme_custom <- function(...){
  theme(
    legend.key.height = unit(1, "mm"),
    legend.key.width = unit(5, "mm"),
    plot.title = element_text(
      size = 24,
      margin = margin(0,0,2,0, "mm")
    ),
    plot.subtitle = element_text(
      size = 18,
      margin = margin(0,0,1,0, "mm")
    ),
    axis.ticks.length = unit(0, "mm"),
    axis.text = element_text(
      margin = margin(0,0,0,0, "mm")
    ),
    panel.grid = element_line(linewidth = 0.2),
    legend.position = "bottom",
    legend.title.position = "top",
    legend.text = element_text(margin = margin(1,0,0,0, "mm")),
    legend.title = element_text(margin = margin(2,0,1,0, "mm")),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box.margin = margin(0,0,0,0, "mm"),
    ...
  )
}

elev <- rast(system.file("raster/elev.tif", package = "spData"))

g10 <- ggplot() +
  geom_spatraster(data = elev) +
  scale_fill_princess_c() +
  labs(
    title = "Base raster",
    subtitle = "Class SpatRaster (elev)",
    fill = NULL
  ) +
  theme_custom()

g11 <- ggplot() +
  geom_sf(
    data = elev |> as.points(values = F) |> st_as_sf(),
    size = 3
  ) +
  scale_colour_princess_c() +
  labs(
    title = "Vectorization without raster values",
    subtitle = "Using terra::as.points(values = FALSE)",
    colour = NULL
  ) +
  theme_custom()

g12 <- ggplot() +
  geom_sf(
    data = elev |> as.points() |> st_as_sf(),
    mapping = aes(colour = elev),
    size = 3
  ) +
  scale_colour_princess_c() +
  labs(
    title = "Vectorization, retaining values as attributes",
    subtitle = "Using terra::as.points() [default behaviour]",
    colour = NULL
  ) +
  theme_custom()


g <- g10 + g11 + g12 +
  plot_layout(nrow = 1) +
  plot_annotation(
    tag_levels = "a",
    title = "Vectorizing rasters into points, with or without values",
    theme = theme(
      plot.title = element_text(
      size = 48,
      margin = margin(1,0,1,0, "mm"),
      hjust = 0.5
      )
    )
  ) &
  theme(
    plot.tag = element_text(
      size = 24,
      margin = margin(0,0,-5,4, "mm"),
      face = "bold"
    ),
    plot.margin = margin(0,0,0,0, "mm")
  )


ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_5_1.png"),
  plot = g,
  height = 700,
  width = 2200,
  units = "px",
  bg = "white"
)
Figure 9: Vectorization of raster into points using terra::as_points() with or without retention of raster values.
  • (2) Contour Lines (Isolines):
    • Represent lines of constant values (e.g., elevation, temperature).
    • Use terra::as.contour() which creates a SpatVector, which then is converted into an {sf} object using sf::st_as_sf() which shows contour lines as LINESTRING or MULTILINESTRING geometry featuers, with an additional column on level (every 100 metres)
    • Or, rasterVis::contourplot() for creating and overlaying contours. For example, visualize DEMs with hillshading and contours in Figure 10
Code 
# Digital elevation showing the southern flank of Mt. Mongón.
dem <- rast(system.file("raster/dem.tif", package = "spDataLarge"))

g13 <- ggplot() +
  geom_spatraster(data = dem) +
  scale_fill_wiki_c() +
  labs(
    title = "Base raster",
    subtitle = "Class SpatRaster (dem)",
    fill = NULL
  ) +
  theme_custom()


dem_contour <- terra::as.contour(dem, nlevels = 10) |> 
  st_as_sf()

g14 <- ggplot() +
  geom_sf(
    data = dem_contour,
    mapping = aes(colour = level),
    linewidth = 1
  ) +
  geom_sf_label(
    data = dem_contour,
    mapping = aes(
      colour = level,
      label = level
    ),
    fill = alpha("white", 0.7),
    fontface = "bold",
    family = "body_font",
    size = 7,
    label.padding = unit(0.1, "lines")
  ) +
  scale_colour_wiki_c(
    breaks = unique(dem_contour$level)
  ) +
  labs(
    title = "Contour Lines with {sf}",
    subtitle = "Using terra::as.contour(nlevels = 10)",
    colour = NULL, x = NULL, y = NULL
  ) +
  theme_custom() +
  theme(
    legend.key.width = unit(7, "mm")
  )

dem_contour <- terra::as.contour(dem, nlevels = 30) |> 
  st_as_sf()

g15 <- ggplot() +
  geom_sf(
    data = dem_contour,
    mapping = aes(colour = level),
    linewidth = 0.5
  ) +
  scale_colour_wiki_c(
    breaks = unique(dem_contour$level)
  ) +
  labs(
    title = "Contour Lines with {sf}",
    subtitle = "Using terra::as.contour(nlevels = 30)",
    fill = NULL, colour = NULL
  ) +
  theme_custom() +
  theme(
    legend.key.width = unit(7, "mm")
  )

g <- g13 + g14 + g15 +
  plot_layout(
    nrow = 1
  ) +
  plot_annotation(
    tag_levels = "a",
    title = "Vectorizing rasters into lines - Contour Lines",
    theme = theme(
      plot.title = element_text(
      size = 48,
      margin = margin(2,0,2,0, "mm"),
      hjust = 0.9
      )
    )
  ) &
  theme(
    plot.tag = element_text(
      size = 24,
      margin = margin(0,0,-5,4, "mm"),
      face = "bold"
    ),
    plot.margin = margin(2,0,2,0, "mm")
  )


ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_5_2_1.png"),
  plot = g,
  height = 800,
  width = 1600,
  units = "px",
  bg = "white"
)

# File name to save the PNG
output_file <- here::here("book_solutions", 
                        "images", 
                        "chapter6_5_2_2.png")

png(filename = output_file, 
    width = 600, 
    height = 800,
    units = "px")

# Generate the contour plot
rasterVis::contourplot(dem, main = "rasterVis::contourplot()")

# Close the graphics device
dev.off()
library(magick)
# Read the two images
img1 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_5_2_1.png"))
img2 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_5_2_2.png"))

c(img1, img2) |> 
  image_append() |> 
  image_write(here::here("book_solutions", 
                        "images", 
                        "chapter6_5_2.png"))

unlink(here::here("book_solutions", 
                  "images", 
                  "chapter6_5_2_1.png"))

unlink(here::here("book_solutions", 
                  "images", 
                  "chapter6_5_2_2.png"))
Figure 10
  • (3) Polygons from Rasters:
    • Use terra::as.polygons() to convert raster cells to polygons with five coordinates per cell.
    • Aggregation: Dissolve boundaries with the same attribute values (aggregate argument). Example: Grain raster converted to polygons in Figure 11.
  • Further enhacement: Polygon Smoothing: Use smoothr::smooth() (Strimas-Mackey 2023) (with method = "chaikin", "ksmooth" or "spline") to remove sharp edges. But remember, Smoothing alters spatial coverage; analyse cautiously.
Code 
grain <- rast(system.file("raster/grain.tif", package = "spData"))

g16 <- ggplot() +
  geom_spatraster(dat = grain) +
  scale_fill_princess_d() +
  labs(
    title = "Base raster; Qualitative",
    subtitle = "With three discrete values",
    fill = NULL
  ) +
  coord_sf(expand = FALSE) +
  theme_custom() +
  theme(legend.position = "none")

g17 <- ggplot() +
  geom_sf(
    data = grain |> as.polygons(aggregate = T) |> st_as_sf(),
    mapping = aes(fill = grain),
    colour = "white",
    linewidth = 1
  ) +
  scale_fill_princess_d() +
  labs(
    title = "Vectorized as aggregated polygons",
    subtitle = "Using terra::as.polygons(aggregate = TRUE)",
    fill = NULL
  ) +
  coord_sf(expand = FALSE) +
  theme_custom() +
  theme(legend.position = "none")

g18 <- ggplot() +
  geom_sf(
    data = grain |> as.polygons(aggregate = F) |> st_as_sf(),
    mapping = aes(fill = grain),
    colour = "white",
    linewidth = 0.5
  ) +
  scale_fill_princess_d() +
  labs(
    title = "Vectorized as polygons",
    subtitle = "Using terra::as.polygons(aggregate = FALSE)",
    fill = NULL
  ) +
  coord_sf(expand = FALSE) +
  theme_custom() +
  theme(legend.position = "none")

g19 <- ggplot() +
  geom_sf(
    data = grain |> 
            as.polygons(aggregate = T) |> 
            st_as_sf() |> 
            smoothr::smooth(),
    mapping = aes(fill = grain),
    colour = "white",
    linewidth = 0.5
  ) +
  scale_fill_princess_d() +
  labs(
    title = "Vectorized as smoothed polygons",
    subtitle = "terra::as.polygons() & smoothr::smooth()",
    fill = NULL
  ) +
  coord_sf(expand = FALSE) +
  theme_custom() +
  theme(legend.position = "none")

g20 <- ggplot() +
  geom_sf(
    data = grain |> 
            as.polygons(aggregate = T) |> 
            st_as_sf() |> 
            smoothr::smooth(method = "spline"),
    mapping = aes(fill = grain),
    colour = "white",
    linewidth = 0.5
  ) +
  scale_fill_princess_d() +
  labs(
    title = "Vectorized as smoothed polygons",
    subtitle = "smoothr::smooth(method = \"spline\")",
    fill = NULL
  ) +
  coord_sf(expand = FALSE, clip = "off") +
  theme_custom() +
  theme(legend.position = "none")

g <- g16 + g17 + g18 + g19 + g20 +
  plot_layout(nrow = 1) +
  plot_annotation(
    title = "Vectorizing rasters into polygons (aggregation)",
    tag_levels = "a",
    theme = theme(
      plot.title = element_text(
      size = 48,
      margin = margin(2,0,2,0, "mm"),
      hjust = 0.5
      )
    )
  ) &
  theme(
    plot.tag = element_text(
      size = 24,
      margin = margin(0,-1,-5,4, "mm"),
      face = "bold"
    ),
    plot.margin = margin(2,0,2,0, "mm")
  )


ggsave(
  filename = here::here("book_solutions", 
                        "images", 
                        "chapter6_5_3.png"),
  plot = g,
  height = 700,
  width = 2200,
  units = "px",
  bg = "white"
)
Figure 11
Code 
library(magick)

# Read the two images
img1 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_4_1.png"))
img2 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_4_2.png"))
img3 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_5_1.png"))
img4 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_5_2.png"))
img5 <- image_read(here::here("book_solutions", 
                        "images", 
                        "chapter6_5_3.png"))

c(img1, img2, img3, img4, img5) |> 
  image_append(stack = TRUE) |> 
  image_write(
    here::here("book_solutions",
               "images",
               "chapter6_5.png")
  )

Compilation Plot with {magick}

A compilation plot of the techniques in the last two sections is shown in Figure 12

Figure 12: A compiled image for all the techniques of raster-vector interaction: Rasterization and Vectorization.

Exercise Solutions

Code 
library(sf)             # Simple Features in R
library(terra)          # Rasters in R (with tidy workflow)
library(tidyterra)      # Tidy workflow and ggplot2 with rasters
library(spData)         # Getting data
library(tidyverse)      # Data Wrangling and {ggplot2}
library(patchwork)      # Compiling Plots

zion_points_path <- system.file(
  "vector/zion_points.gpkg", 
  package = "spDataLarge"
  )

zion_points <- read_sf(zion_points_path)

srtm <- rast(system.file("raster/srtm.tif", package = "spDataLarge"))

ch <- st_combine(zion_points) |>
  st_convex_hull() |> 
  st_as_sf()

sysfonts::font_add_google("Saira Condensed", "body_font")
showtext::showtext_auto()

theme_set(
  theme_minimal(
    base_family = "body_font",
    base_size = 12
  )
)

theme_custom <- function(...){
  theme(
    legend.key.height = unit(1, "mm"),
    legend.key.width = unit(6, "mm"),
    plot.title = element_text(
      size = 24,
      margin = margin(0,0,2,0, "mm"),
      lineheight = 0.3
    ),
    plot.subtitle = element_text(
      size = 18,
      margin = margin(0,0,1,0, "mm"),
      lineheight = 0.3
    ),
    axis.ticks.length = unit(0, "mm"),
    axis.text = element_text(
      margin = margin(0,0,0,0, "mm")
    ),
    panel.grid = element_line(linewidth = 0.2),
    legend.position = "bottom",
    legend.title.position = "top",
    legend.text = element_text(margin = margin(1,0,0,0, "mm")),
    legend.title = element_text(margin = margin(2,0,1,0, "mm")),
    legend.margin = margin(0,0,0,0, "mm"),
    legend.box.margin = margin(0,0,0,0, "mm")
  )
}

E1.

Crop the srtm raster using (1) the zion_points dataset and (2) the ch dataset. Are there any differences in the output maps? Next, mask srtm using these two datasets. Can you see any difference now? How can you explain that?

The Figure 13 (b) shows the srtm raster cropped using zion_points dataset. The Figure 13 (c) shows the raster srtm cropped using ch dataset. There seem to be no differences between the two output cropped maps.

Code 
g1 <- ggplot() +
  geom_spatraster(data = srtm) +
  geom_sf(data = ch, alpha = 0.5, fill = "white") +
  geom_sf(data = zion_points, colour = "blue") +
  scale_fill_wiki_c() +
  coord_sf(
    expand = FALSE,
    xlim = c(-113.25, -112.83),
    ylim = c(37.11, 37.53)
  ) +
  labs(
    title = "Base datasets",
    subtitle = "srtm: Elevation Raster; ch: convex hull;\nzion_points: blue points"
  ) +
  theme_custom()

g2 <- ggplot() +
  geom_spatraster(data = srtm |> crop(zion_points)) +
  # geom_sf(data = ch, alpha = 0.5, fill = "white") +
  geom_sf(data = zion_points, colour = "blue", alpha = 0.5) +
  scale_fill_wiki_c() +
  coord_sf(
    expand = FALSE,
    xlim = c(-113.25, -112.83),
    ylim = c(37.11, 37.53)
  ) +
  labs(
    title = "Cropping by points",
    subtitle = "srtm |> crop(zion_points)"
  ) +
  theme_custom()

g3 <- ggplot() +
  geom_spatraster(data = srtm |> crop(zion_points)) +
  geom_sf(data = ch, alpha = 0.5, fill = "white", colour = "white") +
  # geom_sf(data = zion_points, colour = "blue", alpha = 0.5) +
  scale_fill_wiki_c() +
  coord_sf(
    expand = FALSE,
    xlim = c(-113.25, -112.83),
    ylim = c(37.11, 37.53)
  ) +
  labs(
    title = "Cropping by Convex hull",
    subtitle = "srtm |> crop(ch)"
  ) +
  theme_custom()


g <- g1 + g2 + g3 +
  plot_layout(nrow = 1) +
  plot_annotation(tag_levels = "a")

ggsave(
  filename = here::here("book_solutions", "images",
                        "chapter6-e1-1.png"),
  height = 1000,
  width = 2000,
  units = "px"
)
Figure 13

However, the Figure 14 (a) shows masking of the raster with points - showing that masking by points results in selected points on the raster, thus resulting in no visible result. The Figure 14 (b) shows masking of a raster by a polygon, a more expected result.

Code 
g1 <- ggplot() +
  geom_spatraster(data = srtm |> mask(zion_points)) +
  # geom_sf(data = ch, alpha = 0.5, fill = "white") +
  geom_sf(data = zion_points, colour = "blue", 
          alpha = 0.5, size = 0.5) +
  scale_fill_wiki_c() +
  coord_sf(
    expand = FALSE,
    xlim = c(-113.25, -112.83),
    ylim = c(37.11, 37.53)
  ) +
  labs(
    title = "Masking a raster by points",
    subtitle = "srtm |> mask(zion_points)"
  ) +
  theme_custom()

g2 <- ggplot() +
  geom_spatraster(data = srtm |> mask(ch)) +
  # geom_sf(data = ch, alpha = 0.5, fill = "white") +
  geom_sf(data = zion_points, 
          colour = "blue", alpha = 0.5,
          size = 0.5) +
  scale_fill_wiki_c() +
  coord_sf(
    expand = FALSE,
    xlim = c(-113.25, -112.83),
    ylim = c(37.11, 37.53)
  ) +
  labs(
    title = "Masking a raster by a polygon",
    subtitle = "srtm |> mask(ch)"
  ) +
  theme_custom()


g <- g1 + g2 +
  plot_layout(nrow = 1) +
  plot_annotation(tag_levels = "a")

ggsave(
  filename = here::here("book_solutions", "images",
                        "chapter6-e1-2.png"),
  height = 800,
  width = 1600,
  units = "px"
)
Figure 14

The distinction between masking and cropping can be better understood by examining how they interact with polygons and points:

  1. Masking: When masking with a polygon, the output includes all the raster cells that fall within the interior or lie on the boundaries of the polygon. In contrast, masking with points produces an output that includes only the specific raster cells that align directly with those points, without considering the area around them.

  2. Cropping: Cropping, whether by points or polygons, yields the same result. This is because cropping considers the bounding box, which is the smallest rectangle that can enclose the given geometry. For a set of points and their convex hull (a polygon that encloses the points), the bounding box remains identical. Consequently, the raster output produced by cropping is consistent regardless of whether the input geometry is a polygon or a set of points.

E2.

Firstly, extract values from srtm at the points represented in zion_points. Next, extract average values of srtm using a 90 buffer around each point from zion_points and compare these two sets of values. When would extracting values by buffers be more suitable than by points alone?

Bonus: Implement extraction using the exactextractr package and compare the results.

The Table 1 compares the two values from rasters - extracted from the points and extracted from the buffer zone (90 metres) around the points. It would be better to use the buffer around the points, since it removes the chance error or randomness.

Further, the use of exactextractr package helps us calculate even better weighted mean of the buffer zone, by calculating how much of a raster cell’s proportion fall within the buffer zone.

Code 
df1 <- srtm |> 
  terra::extract(zion_points) |> 
  as_tibble() |> 
  rename(height_on_points = srtm)

df2 <- srtm |> 
  terra::extract(st_buffer(zion_points, dist = 90)) |> 
  group_by(ID) |> 
  summarise(
    mean_height_in_buffer_zone = mean(srtm, na.rm = T)
  )

df3 <- srtm |> 
  exactextractr::exact_extract(st_buffer(zion_points, dist = 90)) |> 
  purrr::imap_dfr(~ .x %>% mutate(ID = .y)) |> 
  as_tibble() |> 
  group_by(ID) |> 
  summarise(
    mean_with_exact_fractions = weighted.mean(value, w = coverage_fraction)
  ) |> 
  mutate(ID = parse_number(as.character(ID)))
Code 
df_all <- structure(list(ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 
29, 30), height_on_points = c(1802L, 2433L, 1886L, 1370L, 1452L, 
1635L, 1380L, 2032L, 1830L, 1860L, 1440L, 2145L, 1942L, 1691L, 
1776L, 2198L, 1820L, 1349L, 1758L, 1424L, 2159L, 1809L, 1826L, 
1550L, 1799L, 2102L, 2118L, 1372L, 1905L, 1574L), mean_height_in_buffer_zone = c(1803.5, 
2418.75, 1875, 1385.66666666667, 1411.33333333333, 1629, 1378.25, 
2040.25, 1812.5, 1863, 1485, 2163.25, 1930, 1689.25, 1770.6, 
2197.75, 1820.25, 1342.25, 1755.5, 1432.25, 2163, 1799.25, 1815.5, 
1557.25, 1818.75, 2084, 2112.75, 1381, 1906.5, 1568.75), mean_with_exact_fractions = c(1802.00133958487, 
2424.45674770282, 1876.08426481862, 1392.71469762041, 1418.32776902309, 
1634.65560885957, 1382.29661439712, 2036.58308765316, 1820.0303014567, 
1862.25146840549, 1474.00216389597, 2156.26597640816, 1928.10327037732, 
1690.96565049072, 1772.48742394502, 2199.49153660086, 1814.37498800353, 
1346.31074580787, 1750.9561283278, 1432.07177081194, 2162.75811241747, 
1799.24374331498, 1818.16628377294, 1555.46331837623, 1815.26746442256, 
2084.60258183394, 2114.31006688883, 1382.04095054598, 1907.24485439106, 
1568.26154119403), difference = c(1.5, -14.25, -11, 15.6666666666667, 
-40.6666666666667, -6, -1.75, 8.25, -17.5, 3, 45, 18.25, -12, 
-1.75, -5.40000000000009, -0.25, 0.25, -6.75, -2.5, 8.25, 4, 
-9.75, -10.5, 7.25, 19.75, -18, -5.25, 9, 1.5, -5.25)), row.names = c(NA, 
-30L), class = c("tbl_df", "tbl", "data.frame"))

df_all |> 
  gt::gt() |> 
  gt::fmt_number(
    columns = -ID,
    decimals = 0
  ) |> 
  gt::cols_label_with(fn = snakecase::to_title_case) |> 
  gt::opt_interactive()
Table 1: A table on the height of points, mean height inside the buffer zone and difference between the two

E3.

Subset points higher than 3100 meters in New Zealand (the nz_height object) and create a template raster with a resolution of 3 km for the extent of the new point dataset. Using these two new objects:

  • Count numbers of the highest points in each grid cell.

  • Find the maximum elevation in each grid cell.

The object above_3100 contains the 19 peaks which are above 3100 metres. Further, the template_raster is a template raster with a resolution of 3 km and the extent of the new point dataset. The results for count and maximum height in each cell os new raster are shown in Figure 15.

Code 
data("nz_height")

above_3100 <- nz_height |> 
  filter(elevation > 3100)

# Create a template raster for using in rasterization
template_raster <- rast(
  # Spat Extent of the Vector data
  terra::ext(above_3100),
  # Required resolution of the raster
  resolution = 3000,
  # CRS for the raster
  crs = crs(above_3100)
)

elevation_raster <- above_3100 |> rasterize(template_raster)
elevation_raster_1 <- above_3100 |> rasterize(template_raster, 
                                            fun = "length")
elevation_raster_2 <- above_3100 |> 
  rasterize(
    template_raster,
    field = "elevation",
    fun = "max")

df1 <- elevation_raster_1 |> 
  values() |> 
  as_tibble() |> 
  mutate(cell_id = row_number()) |> 
  drop_na() |> 
  rename(
    numer_of_peaks = V1_length,
    cell_ID = cell_id
  )

df2 <- elevation_raster_2 |> 
  values() |> 
  as_tibble() |> 
  mutate(cell_id = row_number()) |> 
  drop_na() |> 
  rename(
    maximum_height_of_peaks = max,
    cell_ID = cell_id
  )

g1 <- ggplot() +
  geom_spatraster(data = elevation_raster_1) +
  scale_fill_viridis_c(
    direction = -1,
    na.value = "transparent"
  ) +
  labs(
    title = "Number of peaks per cell",
    subtitle = "Raster created using rasterize(fun = `length`)"
  ) +
  theme_custom() +
  theme(
    legend.position = "bottom"
  ) 

g2 <- ggplot() +
  geom_spatraster(data = elevation_raster_2) +
  scale_fill_princess_c(
    na.value = "transparent"
  ) +
  labs(
    title = "Maximum height per cell",
    subtitle = "Using rasterize(field = `elevation`, fun = `max`)"
  ) +
  theme_custom() +
  theme(
    legend.position = "bottom"
  )

g <- g1 + g2 +
  plot_layout(nrow = 1) +
  plot_annotation(
    tag_levels = "a",
    tag_prefix = "(",
    tag_suffix = ")"
  ) &
  theme(
    plot.tag = element_text(
      size = 36,
      face = "bold"
    )
  )

ggsave(
  filename = here::here("book_solutions", "images",
                        "chapter6-e3.png"),
  height = 800,
  width = 1600,
  units = "px"
)
Figure 15: Figure showing the two rasters (a) Raster with each cell of resolution of 3 km, and colour (fill) depicting the count of peaks, i.e., the number of peaks above 3100 metres within that cell. (b) Raster with each cell of same resolution, where each cell’s colour (fill) depicts the maximum elevation (height) within that each cell.

The count of the number of highest points in each grid cell, and the maximum height in each grid cell is shown in Table 2

Code 
# dput(full_join(df1, df2))

df12 <- structure(list(numer_of_peaks = c(7, 2, 1, 6, 1, 1), cell_ID = c(24L, 
25L, 28L, 31L, 38L, 50L), maximum_height_of_peaks = c(3497, 3194, 
3199, 3724, 3593, 3151)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -6L))

df12 |> 
  relocate(cell_ID) |> 
  gt::gt() |> 
  gt::cols_label_with(fn = snakecase::to_title_case)
Table 2: Number of highest points in each grid cell.
Cell Id Numer of Peaks Maximum Height of Peaks
24 7 3497
25 2 3194
28 1 3199
31 6 3724
38 1 3593
50 1 3151

E4.

Aggregate the raster counting high points in New Zealand (created in the previous exercise), reduce its geographic resolution by half (so cells are 6 x 6 km) and plot the result.

  • Resample the lower resolution raster back to the original resolution of 3 km. How have the results changed?

  • Name two advantages and disadvantages of reducing raster resolution.

The code below aggregates the raster showing the count of high points into a resolution of 6 km, using terra::resample() with a template_raster_2 of resolution 6 km. The results are shown in Figure 16.
Advantages of Reducing Raster Resolution

  1. Improved Computation Efficiency: Fewer cells reduce the data size, speeding up processing and analysis.

  2. Simplification: Provides a coarser overview for large-scale analysis, which may be sufficient for specific applications like regional studies.

Disadvantages of Reducing Raster Resolution

  1. Loss of Detail: Fine spatial variations and features are lost, potentially affecting precision in applications requiring high detail.

  2. Accuracy Trade-Off: Aggregated data may misrepresent localized phenomena, leading to potential inaccuracies in spatial analyses.

Code 
# Create a template raster for using in rasterization
template_raster_2 <- rast(
  # Spat Extent of the Vector data
  terra::ext(above_3100),
  # Required resolution of the raster
  resolution = 6000,
  # CRS for the raster
  crs = crs(above_3100)
)


g1 <- ggplot() +
  geom_spatraster(data = elevation_raster_1) +
  scale_fill_viridis_c(
    direction = -1,
    na.value = "transparent"
  ) +
  labs(
    title = "Number of peaks per cell",
    subtitle = "Raster created using rasterize(fun = `length`)"
  ) +
  theme_custom() +
  theme(
    legend.position = "bottom"
  ) 


g2 <- ggplot() +
  geom_spatraster(
    data = elevation_raster_1 |> 
      # terra::aggregate(fact = 2, fun = mean)
      resample(template_raster_2, method = "bilinear")
  ) +
  scale_fill_viridis_c(
    direction = -1,
    na.value = "transparent"
  ) +
  labs(
    title = "Aggregated raster (res = 6 km)",
    subtitle = "Using resample(method = `bilinear`)"
  ) +
  theme_custom() +
  theme(
    legend.position = "bottom"
  )

g3 <- ggplot() +
  geom_spatraster(
    data = elevation_raster_1 |> 
      resample(template_raster_2, method = "bilinear") |> 
      resample(template_raster, method = "bilinear")
  ) +
  scale_fill_viridis_c(
    direction = -1,
    na.value = "transparent"
  ) +
  labs(
    title = "Disaggregated raster (res = 3 km)",
    subtitle = "Reverting back to high res. leads to loss of data"
  ) +
  theme_custom() +
  theme(
    legend.position = "bottom"
  )

g <- g1 + g2 + g3 +
  plot_layout(nrow = 1) +
  plot_annotation(
    tag_levels = "a",
    tag_prefix = "(",
    tag_suffix = ")"
  ) &
  theme(
    plot.tag = element_text(
      size = 36,
      face = "bold"
    )
  )

ggsave(
  filename = here::here("book_solutions", "images",
                        "chapter6-e4.png"),
  height = 800,
  width = 1800,
  units = "px"
)
Figure 16: Figures showing the (a) Original raster of resolution 3 km, showing number of peaks above 3100 meters in each cell, (b) The raster aggregated to resolution of 6 km, (c) The same aggregted raster disaggregated back to resolution of 3 km, with considerable loss of information and data.

E5.

Polygonize the grain dataset and filter all squares representing clay.

  • Name two advantages and disadvantages of vector data over raster data.

  • When would it be useful to convert rasters to vectors in your work?

The code chunks below show the polygonization, and filtered data of all squares representing clay (using as.polygon(aggregate = FALSE)). The Figure 17 shows the grain dataset and the polygonized grain_sf where each row represents a <POLYGON> (aggregated).

Code 
grain <- rast(system.file("raster/grain.tif", package = "spData"))

grain |>
  terra::as.polygons(aggregate = FALSE) |>
  st_as_sf() |>
  st_cast("POLYGON") |>
  filter(grain == "clay")
Simple feature collection with 10 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -1.5 ymin: -1.5 xmax: 1 ymax: 1.5
Geodetic CRS:  WGS 84
   grain                       geometry
1   clay POLYGON ((-1 1, -1 1.5, -0....
2   clay POLYGON ((-1.5 0.5, -1.5 1,...
3   clay POLYGON ((-0.5 0.5, -0.5 1,...
4   clay POLYGON ((0 0.5, 0 1, 0.5 1...
5   clay POLYGON ((-1.5 0, -1.5 0.5,...
6   clay POLYGON ((0 0, 0 0.5, 0.5 0...
7   clay POLYGON ((0.5 0, 0.5 0.5, 1...
8   clay POLYGON ((-1.5 -0.5, -1.5 0...
9   clay POLYGON ((-1 -0.5, -1 0, -0...
10  clay POLYGON ((0.5 -1.5, 0.5 -1,...
Code 
grain <- rast(system.file("raster/grain.tif", package = "spData"))

grain_sf <- grain |> 
  terra::as.polygons() |> 
  st_as_sf() |> 
  st_cast("POLYGON")

g1 <- ggplot() +
  geom_spatraster(data = grain) +
  scale_fill_princess_d() +
  labs(
    title = "Base raster; Qualitative",
    subtitle = "With three discrete values",
    fill = NULL
  ) +
  coord_sf(expand = FALSE) +
  theme_custom() +
  theme(legend.position = "none")

g2 <- ggplot() +
  geom_sf(
    data = grain_sf,
    mapping = aes(fill = grain),
    colour = "white",
    linewidth = 1
  ) +
  scale_fill_princess_d() +
  labs(
    title = "Vectorized as aggregated polygons",
    subtitle = "terra::as.polygons() |> st_as_sf() |> st_cast(`POLYGON`)",
    fill = NULL
  ) +
  coord_sf(expand = FALSE) +
  theme_custom() +
  theme(legend.position = "none")

g <- g1 + g2 +
  plot_layout(nrow = 1) +
  plot_annotation(
    tag_levels = "a",
    tag_prefix = "(",
    tag_suffix = ")"
  ) &
  theme(
    plot.tag = element_text(
      size = 36,
      face = "bold"
    )
  )

ggsave(
  filename = here::here("book_solutions", "images",
                        "chapter6-e5.png"),
  height = 800,
  width = 1700,
  units = "px"
)
Figure 17: Figure showing the (a) raster data for `grain`, and (b) the polygonized data with as.polygons(aggregate = FALSE)

Advantages of Vector Data Over Raster Data

  1. Precision in Representation: Vector data precisely represents geographic features such as points, lines, and polygons. It is ideal for applications requiring exact boundaries (e.g., property lines, administrative borders).
  2. Efficient Storage for Sparse Data: Vectors store only the coordinates of features, making them more efficient for sparse or irregularly distributed data compared to rasters, which store values for all cells, even those with no data.

Disadvantages of Vector Data Over Raster Data

  1. Complex Analysis: Performing spatial operations (e.g., overlay, buffering) is computationally more complex with vectors, especially with large datasets.
  2. Inefficient for Continuous Data: Vector data struggles to efficiently represent continuous phenomena (e.g., elevation, temperature) that vary over space, for which raster data is better suited.

When to Convert Rasters to Vectors in Your Work

  1. Feature Extraction: If you need to delineate specific features (e.g., contours from a Digital Elevation Model (DEM) or boundaries of high-value areas), converting a raster to vector allows for precise representation and further analysis.

  2. Integration with Vector Data: For projects requiring alignment or combination with existing vector datasets (e.g., administrative boundaries or transportation networks), converting rasters to vectors ensures compatibility.

  3. Visual Clarity in Maps: Vectors are preferred for creating clear and professional maps where discrete features need labeling or editing.

  4. Custom Analysis: If detailed geometric operations like buffering or intersection analysis are needed, converting rasters into vectors (e.g., polygons of land cover classes) enables such workflows.

References

Hijmans, Robert J. 2024. “Terra: Spatial Data Analysis.” https://CRAN.R-project.org/package=terra.
Pebesma, Edzer, and Roger Bivand. 2023. “Spatial Data Science: With Applications in r.” https://doi.org/10.1201/9780429459016.
Strimas-Mackey, Matthew. 2023. “Smoothr: Smooth and Tidy Spatial Features.” https://CRAN.R-project.org/package=smoothr.