Week 5: Discussion Section

NoteSource Materials

The following materials are modified from Geocomputation with R by Robin Lovelace.

Learning Objectives

• Use terra functions aggregate() and resample() to create a new raster
• Use terra functions as.polygons() to convert a raster to polygons (vector)

1. Get Started

• Create a version-controlled R Project
• Create a Quarto document

Let’s load all necessary packages:

library(tidyverse)
library(sf)
library(terra)
library(spData)

You will be working with the following datasets: - High points in New Zealand, obtained from spData - Artificial grain dataset (with three classes: clay, silt, and sand), obtained from spData

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

3. High elevation points in New Zealand/Aotearoa

1. Subset points higher than 3100 meters in nz_height

TipSolution

nz_height3100 <- nz_height |>  
  dplyr::filter(elevation > 3100) # rows where column elevation > 3100 meters

2. Create a template raster with rast(), where the resolution is 3 km x 3 km, with the same extent and CRS as the previous subset

TipSolution

nz_template <- rast(terra::ext(nz_height3100), # Set extent based on nz_height3100
                    resolution = 3000, 
                    crs = terra::crs(nz_height3100)) # Set CRS based on nz_height3100

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

TipSolution

nz_raster <- rasterize(nz_height3100, nz_template, # Convert vector points to raster data
                       fun = "length") # "length" returns a count per cell

plot(nz_raster, main = "Number of elevation points > 3100 in each grid cell")
plot(st_geometry(nz_height3100), add = TRUE)

4. Find the maximum elevation in each grid cell

TipSolution

nz_raster2 <- rasterize(nz_height3100, nz_template, # Convert vector points to raster data
                        fun = max) # "max" returns maximum elevation value per cell

plot(nz_raster2, main = "Maximum elevation in each grid cell ")
plot(st_geometry(nz_height3100), add = TRUE)

5. With the raster template, complete the following:

• Aggregate raster that counts the highest points in New Zealand/Aotearoa, to reduce its geographic resolution by half
• Resample back to the original resolution

TipSolution

nz_raster_low <- aggregate(nz_raster, 
                           fact = 2, # Reduce resolution by combining 2 cells in each direction into larger cells
                           fun = sum, na.rm = TRUE) # Sum values of all cells for resulting elevation value

# Convert resolution back to original 3kmx3km resolution
nz_resample <- resample(nz_raster_low, nz_raster)

plots <- c(nz_raster, nz_resample)
labs <- c("Original 3 x 3 km", "Resample 3 x 3 km")
plot(plots, main = labs)

plot(nz_raster_low, main = "Resample 6 x 6 km")

4. Vectorize grain raster

6. Polygonize grain and filter to only keep squares that represent clay

TipSolution

# Convert raster data to polygon vector data
grain_poly <- as.polygons(grain) %>% 
  st_as_sf()

plot(grain, main = "Grain (Raster)")

plot(grain_poly, main = "Grain (Vector)")

# Subset polygons to only clay
clay <- grain_poly %>% 
  dplyr::filter(grain == "clay")

plot(clay, main = "Clay")