Overview

Objectives

In this exercise will be identified and categorized burn areas after a fire event. The selected study area is Kalamos, Greece. On August 2017, several strong fires were reported in the area causing many damages. For this exercise, Landsat 8 images were used and 2 vegetation indexes are used to detected the burned areas: Normalized Difference Vegetation Index (NDVI), and Normalized Burn Ratio (NBR).

News on the fire event can be found in the following link

Data

For this practice, use the following files:

  • Raster: LC08_L1TP_183033_20170730_20170811_01_T1, LC08_L1TP_183033_20170916_20170929_01_T1 (Landsat 8)

  • Vector: aoi_forest_fires_greece.shp

Procedure

1. Set the working environment

To start working in R, you need first to indicate the program where to store the information. This is the “Working Directory”. To set up the working directory, add the path of the folder where your files are located in the following function:

# Set the working directory ("wd") and assign it to an object (here, "w")
w <- setwd("C:\\yourPath\\yourFolder")

When copying the path, take into account the R spelling:

  • Either change backslash for /

  • Or use copy the path and run readClipboard() in your console

2. Install packages

Functions are stored in Packages that have to be first installed and loaded. Install them using the following R commands or in the button “Packages” in the 4th environment (normally located in the right-down side of your R Studio interface). Packages must be installed only once, whereas libraries should be loaded every time you write a new script.

# Define which packages to install
libraries <- c("rgdal", "raster", "rgeos", "RColorBrewer", "dplyr", "sf", "magrittr", "RCurl")
# Install all the packages at once
install.packages(libraries)
# Open all the libraries at once
lapply(libraries, library, character.only=TRUE)

OR use ipak (Recommended)

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)

Use devtools to install ‘RStoolbox’. This library has some very useful functions, and in this case, we will install it from the Github repository, and not from a CRAN mirror.

# First, install and load devtools
library(devtools)
# Second, install and load RStoolbox
install_github("bleutner/RStoolbox")
library(RStoolbox)

3. Perform atmospheric corrections and calibrations

The study area is located in Kalamos, a coastal holiday spot some 45 km (30 miles) northeast of Athens, where fires burned for many days during the month of August 2017. The images to be analyzed are from Landsat 8 captured on the 30.07.2017 (pre-fire) and 29.09.2017 (post-fire).

Landsat 8 products consist of quantized and calibrated scaled Digital Numbers (DN) representing multispectral image data acquired by both the Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS). They are delivered in 16-bit unsigned integer format, and can be rescaled to the Top Of Atmosphere (TOA) reflectance using radiometric rescaling coefficients provided in the product metadata file (MTL file).

Load Landsat 8 Metadata files

# Import meta-data and bands based on MTL file
# pre-fire
metaData_pre <- readMeta("LC08_L1TP_183033_20170730_20170811_01_T1\\LC08_L1TP_183033_20170730_20170811_01_T1_MTL.txt")
# post-fire
metaData_post <- readMeta("LC08_L1TP_183033_20170916_20170929_01_T1\\LC08_L1TP_183033_20170916_20170929_01_T1_MTL.txt")

# Stack the metadata bands
lsMeta_pre <- stackMeta(metaData_pre)
lsMeta_post <- stackMeta(metaData_post)

Perform the radiometric correction

# More information about the function on: https://bleutner.github.io/RStoolbox/rstbx-docu/radCor.html

# Perform TOA correction
l8_boa_ref_toa_pre <- radCor(lsMeta_pre, metaData_pre, method = "apref")
l8_boa_ref_toa_post <- radCor(lsMeta_post, metaData_post, method = "apref")

# Perform at-surface-reflectance correction with DOS (Chavez decay model)
l8_boa_ref_asr_pre <- radCor(lsMeta_pre, metaData_pre, method = "dos")
l8_boa_ref_asr_post <- radCor(lsMeta_post, metaData_post, method = "dos")

Check the object “l8_boa_ref_asr_pre”. Observe the dimensions, projection, number and name of the bands, as well as their values.

l8_boa_ref_asr_pre
## class      : RasterStack 
## dimensions : 7731, 7591, 58686021, 10  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 609285, 837015, 4192785, 4424715  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      :       B1_sre,       B2_sre,       B3_sre,       B4_sre,       B5_sre,       B6_sre,       B7_sre,       B9_sre,       B10_bt,       B11_bt 
## min values :       0.0000,       0.0000,       0.0000,       0.0000,       0.0000,       0.0000,       0.0000,       0.0000,     147.5171,     141.6706 
## max values :   1.00000000,   1.00000000,   1.00000000,   1.00000000,   1.00000000,   1.00000000,   1.00000000,   0.07266292, 319.98734626, 314.41918945

Very often TOA corrections result in negative reflectance values. Landsat atmospheric correction and surface reflectance retrieval algorithms are not ideal for water bodies due to the inherently low level of water leaving radiance, and the consequential very low signal to noise ratio. Similarly, surface reflectance values greater than 1.0 can be encountered over bright targets such as snow and playas. These are known computational artifacts in the Landsat surface reflectance products.

To check if you have negative values after performing an atmosphere correction, export the results of the previous step as a raster file, using writeRaster. You can load the images in QGIS, and analyze the histogram.

4. Brick bands

A RasterBrick is a multilayered object, and the processing can be more efficient than processing a RasterStack. The former refers to a single file, while the latter is a set of many bands or files. A multi-band satellite is a RasterBrick.

Create a `RasterBrick with the previous data, plot, export and compare it.

Now load all the multispectral bands relevant for this practice, namely, bands: Red, NIR, and SWIR 2 for the fire index; and bands: Blue, Green, Red and NIR for True and False color composition.

# Brick of all bands, pre and post
# Here we will use the Bottom of Atmosphere / at-surface reflectance  
boa_pre <- brick(l8_boa_ref_asr_pre)
boa_post <- brick(l8_boa_ref_asr_post)

# Brick of the bands relevant for the fire detection
pre_st <- subset(boa_pre, 4:7)
pre_st <- dropLayer(pre_st, 3) #drop SWIR 1 band
post_st <- subset(boa_post, 4:7)
post_st <- dropLayer(post_st, 3)
pre_rgbn <- subset(boa_pre, 2:5)
post_rgbn <- subset(boa_post, 2:5)
# Check the content of the RasterStack that you created. Do this for each.  
pre_st
## class      : RasterStack 
## dimensions : 7731, 7591, 58686021, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 609285, 837015, 4192785, 4424715  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
## names      : B4_sre, B5_sre, B7_sre 
## min values :      0,      0,      0 
## max values :      1,      1,      1

To display 3-band color image, we use plotRGB. We have to select the index of bands we want to render in the red, green and blue regions. Associating each spectral band (not necessarily a visible band) to a separate primary colour results in a colour composite image.

Combine the bands in different order to create a true and false composite of the image. For Landsat 8, the false color composite is (5,4,3), which highlights the vegetation very good and helps to distinguish vegetation from other land covers. The function plotRGB can only be used for RasterStack or RasterBrick.

# Define the layout of the images to plot. In this case, we choose a 2x2 layout. 
layout(matrix(c(1,2,3,4), nrow=2, ncol=2, byrow = TRUE))

# Now add the plots one by one in each "cell" location of the layout matrix.
plotRGB(pre_rgbn, r=3, g=2, b=1, axes=TRUE, stretch="lin", main = 
          "True Color. Pre - Fire Image")
plotRGB(post_rgbn, r=3, g=2, b=1, axes=TRUE, stretch="lin", main = 
          "True color. Post - Fire Image")
plotRGB(pre_rgbn, r=4, g=3, b=2, axes=TRUE, stretch="lin", main = 
          "False color. Pre - Fire Image")
plotRGB(post_rgbn, r=4, g=3, b=2, axes=TRUE, stretch="lin", main = 
          "False color. Post - Fire Image")
Figure 1. Top of Atmosphere reflectance of the Forest Fires area in Greece.

Figure 1. Top of Atmosphere reflectance of the Forest Fires area in Greece.

Export layers and check the result also in QGIS.

Attention! Take into account that your file is now a RasterBrick. In order to export it as a multilayered file, include the command “Interleave=Band” in the writeRaster function.

writeRaster(pre_rgbn, datatype="FLT4S", filename = "raster_results\\pre_rgbn.tif", format = "GTiff", overwrite=TRUE, options=c("INTERLEAVE=BAND","COMPRESS=LZW"))
writeRaster(post_rgbn, datatype="FLT4S", filename = "raster_results\\post_rgbn.tif", format = "GTiff", overwrite=TRUE, options=c("INTERLEAVE=BAND","COMPRESS=LZW"))

5. Crop and mask images

Mask the rasters to the shapefile of the Area of Interest (AOI). Use for the this purpose, the objects “pre_st” and “post_st”, because are the ones that contain the bands relevant for the vegetation and burn indexes.

# Load vector file (shapefile)
# The path to the folder of the shapefile should be given relative to the working directory/the location of the script (if no wd was set)
aoi <- st_read("vector", "aoi_forest_fires_greece") #aoi_forest_fires_greece is the name of the newly created aoi
## Reading layer `aoi_forest_fires_greece' from data source `W:\Ulloa\SHK_summer2020\Lehre\fire_landsat\r_script\vector' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 748659.9 ymin: 4230907 xmax: 761460.2 ymax: 4246236
## projected CRS:  WGS 84 / UTM zone 34N

Assign or change the coordinate system of the shapefile

# Reproject aoi if needed
aoi <- st_transform(aoi, crs(boa_pre))

Crop and mask the image to the exact boundaries of the aoi

# Mask the image to the shape of the shapefile (product to be masked, shapefile used)
pre_st_masked <- mask(pre_st, aoi)
post_st_masked <- mask(post_st, aoi)

# Crop the image to the extent of the shapefile (product to be cropped, shapefile used)
pre_st_cr <- crop(pre_st_masked, aoi) 
post_st_cr <- crop(post_st_masked, aoi)

Visualize your output. Why does the graph function produce 3 plots?

plot(pre_st_cr)
Figure 2. RasterStack masked of the Pre-Fire event. Bands 4, 5, 7.

Figure 2. RasterStack masked of the Pre-Fire event. Bands 4, 5, 7.

Export the cropped images

Attention! Take into account that your file is now a RasterBrick. In order to export it as a multilayered file, include the command “Interleave=Band” in the writeRaster function.

writeRaster(pre_st_cr, datatype="FLT4S", filename = "raster_results\\pre_st_cr.tif", format = "GTiff", overwrite=TRUE, options=c("INTERLEAVE=BAND","COMPRESS=LZW"))
writeRaster(post_st_cr, datatype="FLT4S", filename = "raster_results\\post_st_cr.tif", format = "GTiff", overwrite=TRUE, options=c("INTERLEAVE=BAND","COMPRESS=LZW"))