Week 4: Discussion Section - Answer Key

Practice raster manipulations

Author

Allie Caughmam

Plot a Histogram and Boxplot of dem

library(tidyverse)
library(spData)
library(spDataLarge)
library(sf)
library(stars)
library(terra)
library(kableExtra)
dem <- terra::rast(system.file("raster/dem.tif", package = "spDataLarge"))
landsat <- terra::rast(system.file("raster/landsat.tif", package = "spDataLarge"))
srtm <- terra::rast(system.file("raster/srtm.tif", package = "spDataLarge"))
peru <- stars::read_stars(here::here("data", "week4-discussion", "PER_elv.tif"))

Plot a histogram and boxplot of dem

hist(dem,
     main = "Digital Elevation Model Raster Value Distribution",
     xlab = "Value")

boxplot(dem,
        main = "Digital Elevation Model Raster Value Distribution",
        ylab = "Value")

Reclassify Elevation and Find Mean

# define a reclassification matrix
rcl <- matrix(c(-Inf, 300, 0, # values -Inf to 300 = 0
                300, 500, 1,  # values 300 to 500 = 1
                500, Inf, 2), # values 500 to Inf = 2
              ncol = 3, byrow = TRUE)

# apply the matrix to reclassify the raster, making all cells 0 or 1 or 2
dem_rcl <- terra::classify(dem, rcl = rcl)

# assign labels to the numerical categories
levels(dem_rcl) <- tibble::tibble(id = 0:2, 
                                  cats = c("low", "medium", "high"))

# calculate mean elevation for each category using original DEM values
elevation_mean <- terra::zonal(dem, dem_rcl, fun = "mean")
elevation_mean

Find Correlation Between NDWI and NDVI

Define functions for calculating NDWI and NDVI

ndwi_fun <- function(green, nir){
    (green - nir)/(green + nir)
}

ndvi_fun <- function(nir, red){
  (nir - red)/(nir + red)
}

Apply the functions to the appropriate Landsat 8 bands. Landsat 8 bands 2-5 correspond to bands 1-4 for this raster. Bands are as follows:

Band Color
1 blue 30 meter
2 green 30 meter
3 red 30 meter
4 near-infrared 30 meter 
ndwi_rast <- terra::lapp(landsat[[c(2, 4)]],
                         fun = ndwi_fun)
plot(ndwi_rast,
     main = "Zion National Park NDWI")

ndvi_rast <- terra::lapp(landsat[[c(4, 3)]],
                         fun = ndvi_fun)

# stack rasters
combine <- c(ndvi_rast, ndwi_rast)

plot(combine, main = c("NDVI", "NDWI"))

# calculate the correlation between raster layers 
terra::layerCor(combine, fun = cor)

Find Distances Across All peru Cells

# Aggregate by a factor of 20 to reduce resolution and create new raster
peru_agg <- terra::aggregate(rast(peru), fact = 20)
plot(peru_agg)

# Create mask of ocean (NA values)
water_mask <- is.na(peru_agg) # returns TRUE value for NA
# Set all FALSE values to NA
water_mask[water_mask == 0] <- NA
plot(water_mask)

# Find distance from each cell to ocean/coastline (default is unit = "m")
distance_to_coast <- terra::distance(water_mask)
# Convert from meters to kilometers 
distance_to_coast_km <- distance_to_coast/1000

plot(distance_to_coast_km, main = "Distance to the coast (km)")

Change Resolution of srtm

plot(srtm)

rast_template <- terra::rast(terra::ext(srtm), res = 0.01)

srtm_resampl1 <- terra::resample(srtm, y = rast_template, method = "bilinear")
srtm_resampl2 <- terra::resample(srtm, y = rast_template, method = "near")
srtm_resampl3 <- terra::resample(srtm, y = rast_template, method = "cubic")
srtm_resampl4 <- terra::resample(srtm, y = rast_template, method = "cubicspline")
srtm_resampl5 <- terra::resample(srtm, y = rast_template, method = "lanczos")

srtm_resampl_all <- c(srtm_resampl1, srtm_resampl2, srtm_resampl3, srtm_resampl4, srtm_resampl5)
labs <- c("Bilinear", "Near", "Cubic", "Cubic Spline", "Lanczos")
plot(srtm_resampl_all, main = labs)