Week 6: Discussion Section

library(tidyverse)
library(sf)
library(terra)
library(tmap)
library(tmaptools)

# Set directory for folder
pre_fire_dir <- here::here("data", "LC80340322016189-SC20170128091153")

# Create a list of all images that have the extension .tif and contain the word band
pre_fire_bands <- list.files(pre_fire_dir,
                             pattern = glob2rx("*band*.tif$"),
                             full.names = TRUE)
# Create a raster stack
pre_fire_rast <- rast(pre_fire_bands)

# Read mask raster
pre_mask <- rast(here::here("data", "LC80340322016189-SC20170128091153", "LC80340322016189LGN00_cfmask_crop.tif"))

# Set directory for folder
post_fire_dir <- here::here("data", "LC80340322016205-SC20170127160728")

# Create a list of all images that have the extension .tif and contain the word band
post_fire_bands <- list.files(post_fire_dir,
                             pattern = glob2rx("*band*.tif$"),
                             full.names = TRUE)
# Create a raster stack
post_fire_rast <- rast(post_fire_bands)

# Read mask raster
post_mask <- rast(here::here("data", "LC80340322016189-SC20170128091153", "LC80340322016205LGN00_cfmask_crop.tif"))

nbr_fun <- function(nir, swir2){
    (nir - swir2)/(nir + swir2)
}

Rename Bands

bands <- c("Aerosol", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2")
names(pre_fire_rast) <- bands
names(post_fire_rast) <- bands

Mask Clouds and Shadows

# Set all cells with values greater than 0 to NA
pre_mask[pre_mask > 0] <- NA

# Subset raster based on mask
pre_fire_rast <- mask(pre_fire_rast, mask = pre_mask)

# View raster
plot(pre_fire_rast)

# Set all cells with values greater than 0 to NA
post_mask[post_mask > 0] <- NA

# Subset raster based on mask
post_fire_rast <- mask(post_fire_rast, mask = post_mask)

# View raster
plot(post_fire_rast)

Plot True Color Composite

plotRGB(pre_fire_rast, r = 4, g = 3, b = 2, stretch = "lin", colNA = "black")

plotRGB(post_fire_rast, r = 4, g = 3, b = 2, stretch = "lin", colNA = "black")

Plot False Color Composite

plotRGB(pre_fire_rast, r = 7, g = 5, b = 3, stretch = "lin", colNA = "black")

plotRGB(post_fire_rast, r = 7, g = 5, b = 3, stretch = "lin", colNA = "black")

Calculate NBR and dNBR

pre_nbr_rast <- terra::lapp(pre_fire_rast[[c(5, 7)]], fun = nbr_fun)

plot(pre_nbr_rast, main = "Cold Springs Pre-Fire NBR", colNA = "black")

post_nbr_rast <- terra::lapp(post_fire_rast[[c(5, 7)]], fun = nbr_fun)

plot(post_nbr_rast, main = "Cold Springs Post-Fire NBR", colNA = "black")

Calculate and Plot dNBR

diff_nbr <- pre_nbr_rast - post_nbr_rast

tm_shape(diff_nbr) +
  tm_raster(style = "equal", n = 6, 
            palette = get_brewer_pal("YlOrRd", n = 6, plot = FALSE),
            title = "Difference NBR (dNBR)", colorNA = "black") +
  tm_layout(legend.outside = TRUE)

Reclassification

# Set categories for severity levels
categories <- c("Enhanced Regrowth", "Unburned", "Low Severity", "Moderate Severity", "High Severity")

# Create reclassification matrix
rcl <- matrix(c(-Inf, -0.1, 1, # group 1 ranges for Enhanced Regrowth
                -0.1, 0.1, 2, # group 2 ranges for Unburned
                0.1, 0.27, 3, # group 3 ranges for Low Severity
                0.27, 0.66, 4, # group 4 ranges for Moderity Severity
                0.66, Inf, 5), # group 5 ranges for High Severity
                ncol = 3, byrow = TRUE)

# Use reclassification matrix to reclassify dNBR raster
reclassified <- classify(diff_nbr, rcl = rcl)

reclassified[is.nan(reclassified)] <- NA

tm_shape(reclassified) +
  tm_raster(style = "cat",
            labels = c(categories, "Missing"),
            palette = get_brewer_pal("YlOrRd", n = 5, plot = FALSE),
            title = "Severity Level", colorNA = "black")+
  tm_layout(legend.outside = TRUE)