Week 5: Discussion Section

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

# Load raster data representing grain sizes with the three classes clay, silt and sand
grain <- terra::rast(system.file("raster/grain.tif", package = "spData"))

Subset Points in New Zealand/Aotearoa

# Subset New Zealand elevation points to > 3100 meters
nz_height3100 <- nz_height %>% 
  dplyr::filter(elevation > 3100)

# Create template: define the extent, resolution, and CRS based on nz_height3100
nz_template <- rast(terra::ext(nz_height3100), 
                    resolution = 3000, 
                    crs = terra::crs(nz_height3100))

Count Points in Each Grid Cell

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

plot(nz_raster, main = "Number of Elevation Points > 3100 in Each Grid Cell")
plot(st_geometry(nz_height3100), add = TRUE)

Find Maximum Elevation in Each Grid Cell

# function "max" returns maximum elevation value per cell
nz_raster2 <- rasterize(nz_height3100, nz_template, field = "elevation", fun = max)

plot(nz_raster2, main = "Maximum Elevation in Each Grid Cell ")
plot(st_geometry(nz_height3100), add = TRUE)

Aggregate and Resample Raster

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

# Convert the new raster's resolution back to the 3kmx3km resolution of original raster
nz_resample <- resample(nz_raster_low, nz_raster)

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

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

Vectorize Raster

# 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")