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