EDS 223: Discussion Week 6 - Vector Raster Interactions

introduction

The following exercises are modified from Chapter 6 of Geocomputation with R by Rovin Lovelace.

prerequisites

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

zion_points = read_sf(system.file("vector/zion_points.gpkg", package = "spDataLarge"))
srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))

ch = st_combine(zion_points) |>
  st_convex_hull() |> 
  st_as_sf()

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

plot(srtm)
plot(st_geometry(zion_points), add = TRUE)
plot(ch, add = TRUE)

Exercises

1. 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?

srtm_crop1 = crop(srtm, zion_points)
srtm_crop2 = crop(srtm, ch)
plot(srtm_crop1)
plot(srtm_crop2)

srtm_mask1 = mask(srtm, zion_points)
srtm_mask2 = mask(srtm, ch)
plot(srtm_mask1)

plot(srtm_mask2)

2. 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:

nz_height3100 = nz_height %>% 
  filter(elevation > 3100)

nz_template = rast(ext(nz_height3100), resolution = 3000, crs = crs(nz_height3100))

2a. Count numbers of the highest points in each grid cell.

nz_raster = rasterize(nz_height3100, nz_template, 
                       field = "elevation", fun = "length")
plot(nz_raster)
plot(st_geometry(nz_height3100), add = TRUE)

2b. Find the maximum elevation in each grid cell.

nz_raster2 = rasterize(nz_height3100, nz_template, 
                       field = "elevation", fun = max)
plot(nz_raster2)
plot(st_geometry(nz_height3100), add = TRUE)

3. Aggregate the raster counting high points in New Zealand (created in the previous exercise), reduce its geographic resolution by half (so cells are 6 by 6 km) and plot the result. Resample the lower resolution raster back to the original resolution of 3 km.

nz_raster_low = raster::aggregate(nz_raster, fact = 2, fun = sum, na.rm = TRUE)

nz_resample = resample(nz_raster_low, nz_raster)
plot(nz_raster_low)

plot(nz_resample) # the results are spread over a greater area and there are border issues

plot(nz_raster)

4. Polygonize the grain dataset and filter all squares representing clay.

grain_poly = as.polygons(grain) %>% 
  st_as_sf()

plot(grain)

plot(grain_poly)

clay = grain_poly %>% 
  filter(grain == "clay")
plot(clay)