Overview

Objectives

This script is a continuation of the script “Forest Fires with Landsat 8. Part I”. For this exercise, use Landsat 8 RasterStacks from the previous practice and calculate: Normalized Difference Vegetation Index (NDVI), and Normalized Burn Ratio (NBR).

Data

For this practice, use the following files:

  • Raster: pre_st_cr.tif, post_st_cr.tif

Procedure

1. Install and/or load necessary packages

As for every new session of R, you should set your working directory and load the necessary libraries:

# Install and load the necessary packages 
ipak <- function(pkg){
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    if (length(new.pkg)) 
        install.packages(new.pkg, dependencies = TRUE, repos = "http://cran.us.r-project.org")
    sapply(pkg, require, character.only = TRUE)
}

packages <- c("rgdal", "raster", "rgeos", "RColorBrewer", "dplyr","sf", "magrittr", "RCurl")
ipak(packages)

2. Load rasters calculated in the previous section: “Forest Fires with Landsat 8. Part I”

Set up your working directory.

# Set the working directory ("wd")
w <- setwd("YOURDRIVE\\yourpath\\yourfolder")

Load the required rasters for this section. The correspondent rasters should have been generated with the script of Part I. Since the files were exported as a RasterBrick, use the brick function to load them.

# Load rasters pre and post fire event
pre_st_cr <- brick("raster_results\\pre_st_cr.tif")
post_st_cr <- brick("raster_results\\post_st_cr.tif")

3. Calculate NDVI

Vegetation indices derived from Earth observation satellites are important for a wide range of applications such as vegetation monitoring, drought studies, agricultural activities, climate and hydrologic modelling. Spectral indices make use of the characteristics in the spectrum of the respective material of interest.

For example, the Normalized Difference Vegetation Index (NDVI) makes use of the red and near infrared bands, since the energy reflected in these wavelengths is clearly related to the amount of vegetation cover on the ground surface.

The formula for calculating the NDVI is the following:

\(NDVI = (NIR-RED)/(NIR+RED)\)

NDVI values are standardized and range between -1 to +1. Pixels having NDVI values greater than 0.3 are very probably vegetation.

# Rename bands using the masked images from the previous step.
red_pre = pre_st_cr$pre_st_cr.1 # Pay attention which band was stored in which position!!!
nir_pre = pre_st_cr$pre_st_cr.2
red_post = post_st_cr$post_st_cr.1
nir_post = post_st_cr$post_st_cr.2

#  If you are using the rasters from Part I, the naming is different. For example: 
# red_pre = pre_st_cr$LC08_L1TP_183033_20170730_20170811_01_T1_B4
# red_post = post_st_cr$LC08_L1TP_183033_20170916_20170929_01_T1_B4

# Calculate NDVI for the pre and post fire images
ndvi_pre <- (nir_pre - red_pre)/(nir_pre + red_pre)
ndvi_post <- (nir_post - red_post)/(nir_post + red_post)

Thresholding

The non-vegetation pixels (<0.3) can be masked out. The function calc is used to make arithmetic calculations on rasters (among others applications). In this case, a function is created to select all values smaller than 0.3, and make them equal to no data (NA). The resulting image will have values greater than 0.3. This process is called “Thresholding”.

# Mask out the values smaller than 0.3
veg_pre <- calc(ndvi_pre, function(x){x[x < 0.3] <- NA; return(x)})
veg_post <- calc(ndvi_post, function(x){x[x < 0.3] <- NA; return(x)})

Plot both NDVI outputs. In the second plot all colored pixels are very probably vegetation.

# Indicate the order of the plots. With 'par()' you can indicate in how many columns and rows the images have to be arranged
par(mfrow=c(2,2))

# Set up the limits of the values scale (y axis)
breakpoints <- seq(-1, 1, by=0.2)

# Indicate the 4 images that will be plotted in the  2X2 matrix
plot(ndvi_pre, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI Pre-Fire')
plot(ndvi_post, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI Post-Fire')
plot(veg_pre, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI Threshold, Pre-Fire')
plot(veg_post, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI Threshold Post-Fire')
Figure 1. NDVI of the AOI, with and without thresholding

Figure 1. NDVI of the AOI, with and without thresholding

Export files

writeRaster(ndvi_pre, datatype="FLT4S", filename = "raster_results\\ndvi_pre.tif", format = "GTiff", overwrite=TRUE)
writeRaster(ndvi_post, datatype="FLT4S", filename = "raster_results\\ndvi_post.tif", format = "GTiff", overwrite=TRUE)

4. Calculate Normalized burn ratio (NBR)

The Normalized burn ratio (NBR) is used to identify burned areas. This index uses the near-infrared (NIR) band where plants have a high absorption, and the shortwave infrared (SWIR) band where scarred woody vegetation has a higher reflection.

\(NBR = (NIR-SWIR)/(NIR+SWIR)\)

# The NIR-band was already renamed before. Rename the swir2 band that is needed. 
swir2_pre = pre_st_cr$pre_st_cr.3

swir2_post = post_st_cr$post_st_cr.3

The NBR index was originally developed for use with Landsat TM and ETM+ bands 4 and 7, but it will work with any multispectral sensor with a NIR band between 760 - 900 nm and a SWIR band between 2080 - 2350 nm. Therefore, for Landsat 8, the band used is the SWIR2.

#  Calculate the NBR Index
nbr_pre <- (nir_pre - swir2_pre)/(nir_pre + swir2_pre) 
nbr_post <- (nir_post - swir2_post)/(nir_post + swir2_post) 
# Export the NBR rasters
writeRaster(nbr_pre, datatype="FLT4S", filename = "raster_results\\nbr_pre.tif", format = "GTiff", overwrite=TRUE)
writeRaster(nbr_post, datatype="FLT4S", filename = "raster_results\\nbr_post.tif", format = "GTiff", overwrite=TRUE)

Visualize the NBR rasters before and after the fire.

# Define color palette
nbr_colors <- colorRampPalette(brewer.pal(11, "RdYlGn"))(100) # Check!: http://colorbrewer2.org/

# Define in how many rows and columns are the graphs plotted
par(mfrow=c(1,2))

# Plot
plot(nbr_pre,
     main = "Landsat derived NBR\n Pre-Fire",
     axes = FALSE,
     box = FALSE,
     col = nbr_colors,
     zlim = c(-1, 1))
plot(nbr_post,
     main = "Landsat derived NBR\n Post-Fire",
     axes = FALSE,
     box = FALSE,
     col = nbr_colors,
     zlim = c(-1, 1))
Figure 2. NBR calculations for the Landsat scenes before and after the fire

Figure 2. NBR calculations for the Landsat scenes before and after the fire

5. Calculate Difference Normalized Burn Index (dNBR)

The Normalized Burn Ratio is more meaningful for understanding fire extent and severity when used after calculating the difference between pre and post fire conditions. This difference is best measured using data collected immediately before the fire and then immediately after the fire.

\(dNBR = (prefireNBR - postfireNBR)\)

# Calculate dNBR
dNBR <- (nbr_pre) - (nbr_post)

Reclassification of the ranges of dNBR

Classify the burn severity map according to the following table. It is important to highlight that although those are quantitative values of burn severity, they can have different interpretations and consequences in real life.

Figure 3. Table of Burn severity levels obtained calculating dNBR, proposed by USGS

# Multiply by 1000 to apreciate the differences better
dNBR_scaled <- 1000*dNBR 

# Sets the ranges that will be used to classify dNBR information about the ranges used
NBR_ranges <- c(-Inf, -500, -1, -500, -251, 1, -251, -101, 2, -101, 99, 3, 99, 269, 4, 
                269, 439, 5, 439, 659, 6, 659, 1300, 7, 1300, +Inf, -1) 

# Sets a classification matrix
class.matrix <- matrix(NBR_ranges, ncol = 3, byrow = TRUE)
 
# Classification matrix is used to classify dNBR_scaled
dNBR_reclass <- reclassify(dNBR_scaled, class.matrix, right=NA)

Visualize the dNBR

# Build the legend for the burn severity map

# Build the attribute table for the legend 
dNBR_reclass <- ratify(dNBR_reclass) # define the raster as a categorial variable
rat <- levels(dNBR_reclass)[[1]] # rat=raster attribute table. dataframe

# Creates the text that will be on the legend
rat$legend  <- c("NA", "Enhanced Regrowth, High", "Enhanced Regrowth, Low", "Unburned", 
                 "Low Severity", "Moderate-low Severity", "Moderate-high Severity", "High Severity") 
levels(dNBR_reclass) <- rat 

# Setting the colors for the severity map
my_col=c("white", "darkolivegreen","darkolivegreen4","limegreen", "yellow2", "orange2", "red", "purple")
# Plots the burn severity map with title in two lines; first line: 'Burn Severity'; second line: 'Kalamos, Greece'.
plot(dNBR_reclass,col=my_col,legend=F,box=F,axes=F, main="Burn Severity \n Kalamos, Greece")

# Plots the legend on the right side of the burn severity map
legend(x='right', legend =rat$legend, fill = my_col, y='right')
Figure 4. Categorized raster of burn severity after fire

Figure 4. Categorized raster of burn severity after fire

Export layers

writeRaster(dNBR_reclass, datatype="FLT4S", filename = "raster_results\\dNBR_reclass.tif", format = "GTiff", overwrite=TRUE)

Graph of value distribution per land category

We can make a distribution of the classified NBR values from the previous raster

barplot(dNBR_reclass,
        main = "Distribution of Classified NBR Values",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)
Figure 5. Distribution of pixel values of each category for the selected AOI. Burn severity levels: A =NA, B=Enhanced Regrowth, High, C=Enhanced Regrowth, Low, D=Unburned, E=Low Severity, F=Moderate-low Severity, G=Moderate-high Severity, H=High Severity

Figure 5. Distribution of pixel values of each category for the selected AOI. Burn severity levels: A =NA, B=Enhanced Regrowth, High, C=Enhanced Regrowth, Low, D=Unburned, E=Low Severity, F=Moderate-low Severity, G=Moderate-high Severity, H=High Severity

Graphs with added information

vals <- getValues(dNBR_reclass)

pix_na <- length(subset(vals, vals == -1))
pix_1 <- length(subset(vals, vals == 1))
pix_2 <- length(subset(vals, vals == 2))
pix_3 <- length(subset(vals, vals == 3))
pix_4 <- length(subset(vals, vals == 4))
pix_5 <- length(subset(vals, vals == 5))
pix_6 <- length(subset(vals, vals == 6))
pix_7 <- length(subset(vals, vals == 7))

pix_values <- c(pix_na,pix_1,pix_2,pix_3,pix_4,pix_5,pix_6,pix_7)

bp <- barplot(dNBR_reclass,
        main = "Number of pixels by NBR category",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

text(bp, y=2500, labels=pix_values) #y chosen manually
Figure 6. Distribution of pixel values of each category for the selected AOI with number of pixels per class.

Figure 6. Distribution of pixel values of each category for the selected AOI with number of pixels per class. 

pix_total <- length(vals)

area_na_perc <- round(pix_na/pix_total*100, digits = 5)
area_1_perc <- round(pix_1/pix_total*100, digits = 5)
area_2_perc <- round(pix_2/pix_total*100, digits = 5)
area_3_perc <- round(pix_3/pix_total*100, digits = 5)
area_4_perc <- round(pix_4/pix_total*100, digits = 5)
area_5_perc <- round(pix_5/pix_total*100, digits = 5)
area_6_perc <- round(pix_6/pix_total*100, digits = 5)
area_7_perc <- round(pix_7/pix_total*100, digits = 5)

area_perc_covered <- c(area_na_perc,area_1_perc,area_2_perc,area_3_perc,area_4_perc,area_5_perc,area_6_perc,area_7_perc)

bp <- barplot(dNBR_reclass,
        main = "Percentage of every NBR category",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

text(bp, y=2500, labels=area_perc_covered) #y chosen manually
Figure 7. Distribution of pixel values of each category for the selected AOI with percentage of area covered

Figure 7. Distribution of pixel values of each category for the selected AOI with percentage of area covered 

ex <- extent(dNBR_reclass)
# calculate extent of the raster in square kilometers
x <- (ex[2]-ex[1])/1000
y <- (ex[4]-ex[3])/1000
area_total <- x*y

area_na <- round((area_na_perc/100)*area_total, digits = 5)
area_1 <- round((area_1_perc/100)*area_total, digits = 5)
area_2 <- round((area_2_perc/100)*area_total, digits = 5)
area_3 <- round((area_3_perc/100)*area_total, digits = 5)
area_4 <- round((area_4_perc/100)*area_total, digits = 5)
area_5 <- round((area_5_perc/100)*area_total, digits = 5)
area_6 <- round((area_6_perc/100)*area_total, digits = 5)
area_7 <- round((area_7_perc/100)*area_total, digits = 5)

area_covered <- c(area_na,area_1,area_2,area_3,area_4,area_5,area_6,area_7)

bp <- barplot(dNBR_reclass,
        main = "Area by NBR category",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

text(bp, y=2500, labels=area_covered) #y chosen manually
Figure 8. Distribution of pixel values of each category for the selected AOI with area covered in square kilometers

Figure 8. Distribution of pixel values of each category for the selected AOI with area covered in square kilometers

6. Extra task

Calculate another version of the NBR index, using the 2 SWIR Bands. See link for more information.