Vegetation changes detection after a Wind storm. Case study, Passau - Germany
Part I
Part I
Yrneh Ulloa
2 November 2018
1. Background
On August 19, 2017, a storm destroyed in very little time many forest fields. The disaster took place in the Community of Passau, Germany. It caused more than 40 million Euros in damages. More information can be found in the following links:
In this practice, the areas damaged by the storm will be detected. A mixture of satellite images from different sensors, resolutions will be used. For the detection of changes, the Normalized Difference Vegetation Index and a Sueprvised Classification will be used.
2. Processing of Sentinel Images
Sentinel optical images, i.e. Sentinel-2 (S2), can be downloaded free of charge from https://scihub.copernicus.eu/dhus/#/home (Figure 1). In the search field, different parameters can be specified, like date, cloud coverage, type of Sentinel imagery and location. Once you have created your user account, you are ready to use the data available.
Figure 1. Website of Copernicus for the download of diverse products
3. 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:\\Users\\ulloa-to\\Documents\\Kurse_Seminaer\\Master Fernerkundung\\classes_ulloa\\passau_r")# Check the location of the working directory
getwd()## [1] "C:/Users/ulloa-to/Documents/Kurse_Seminaer/Master Fernerkundung/classes_ulloa/passau_r"When copying the path, take into account the R spelling:
* Either change for /
* Or use copy the path and run readClipboard() in your console
4. 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.
# Install packages function when it is not necessary to specify the folder (because this path was saved in the Global options of R)
install.packages("rgdal")
install.packages("gdalUtils")
# Install multiple packages simultaneously
install.packages(c("sp","raster")) #to install more than one package at the same time # Load the packages
library(sp)
library(rgdal)
library(raster)
library(gdalUtils)5. Load and explore Rasters
Sentinel images can be read with the RGDAL Package. If your version is not recent enough (it should have the JPEG2000 Driver), you can also export the .jp2 files into .tif from another programm (for example, ArcMap, Quantum GIS, GRASS, Erdas, among others).
In this example, Sentinel-2 images were downloaded from the dates: 20.06.2017 (before) and 29.08.2017 (after). The Sentinel-2 scenes used are processed to level 2A. Figure 2 shows the differences between the two level products.
Figure 2. Sentinel 2 products. Source: https://sentinel.esa.int
Sentinel-2, level 2A images are grouped in 3 different sets according to the spatial resolution: 10m, 20m, 60m (see figure 3).
Figure 3. Sentinel 2 bands. Source: https://sentinel.esa.int
For the purpose of this exercise, a RasterStack with the 10m resolution images is enough. This includes the bands 2 (blue), 3 (green), 4 (red) and 8 (nir). Compare the equivalence of the Sentinel-2 bands with the ones from Landsat-7 and Landsat-8 (Figure 4).
Figure 4. Sentinel-2, Landsat-7 and Lansat-8 bands comparison. Source: https://sentinel.esa.int
# Read JP2000 files (to import them directly)
# Images pre-disaster
s2a_20170620_b2_raw <- readGDAL("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B02_10m.jp2")
s2a_20170620_b3_raw <- readGDAL("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B03_10m.jp2")
s2a_20170620_b4_raw <- readGDAL("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B04_10m.jp2")
s2a_20170620_b8_raw <- readGDAL("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B08_10m.jp2")
# Images post-disaster
s2a_20170829_b2_raw <- readGDAL("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B02_10m.jp2")
s2a_20170829_b3_raw <- readGDAL("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B03_10m.jp2")
s2a_20170829_b4_raw <- readGDAL("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B04_10m.jp2")
s2a_20170829_b8_raw <- readGDAL("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B08_10m.jp2")Explore the properties of one of the bands of the image. Compare the before and after images
# Complete description of object
s2a_20170620_b2_raw #band nr.2
s2a_20170829_b2_raw# Number of rows, columns and layers
dim(s2a_20170620_b2_raw)
dim(s2a_20170829_b2_raw)The following functions are useful to see the structure of the rasters. For example, you can compare the projection and the pixel value distribution of both images.
# Structure of object
str(s2a_20170620_b2_raw)
str(s2a_20170829_b2_raw)# Summary of object. Identify the pixel size in the results.
summary(s2a_20170620_b2_raw)## Object of class SpatialGridDataFrame
## Coordinates:
## min max
## x 300000 409800
## y 5290200 5400000
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Grid attributes:
## cellcentre.offset cellsize cells.dim
## x 300005 10 10980
## y 5290205 10 10980
## Data attributes:
## band1
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 0.0
## Mean : 222.8
## 3rd Qu.: 236.0
## Max. :20669.0summary(s2a_20170829_b2_raw)## Object of class SpatialGridDataFrame
## Coordinates:
## min max
## x 300000 409800
## y 5290200 5400000
## Is projected: TRUE
## proj4string :
## [+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Grid attributes:
## cellcentre.offset cellsize cells.dim
## x 300005 10 10980
## y 5290205 10 10980
## Data attributes:
## band1
## Min. : 0.0
## 1st Qu.: 0.0
## Median : 0.0
## Mean : 146.2
## 3rd Qu.: 119.0
## Max. :18412.0# Transform raster from JP2000 to GeoTIFF
# Pre-disaster images
gdal_translate("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B02_10m.jp2", "raster\\s2a_20170620_b2_raw.tif")
gdal_translate("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B03_10m.jp2", "raster\\s2a_20170620_b3_raw.tif")
gdal_translate("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B04_10m.jp2", "raster\\s2a_20170620_b4_raw.tif")
gdal_translate("raster\\L2A_T33UUP_A010415_20170620T100453\\IMG_DATA\\R10m\\L2A_T33UUP_20170620T100031_B08_10m.jp2", "raster\\s2a_20170620_b8_raw.tif")
# Post-disaster images
gdal_translate("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B02_10m.jp2", "raster\\s2a_20170829_b2_raw.tif")
gdal_translate("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B03_10m.jp2", "raster\\s2a_20170829_b3_raw.tif")
gdal_translate("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B04_10m.jp2", "raster\\s2a_20170829_b4_raw.tif")
gdal_translate("raster\\L2A_T33UUP_A011416_20170829T100026\\IMG_DATA\\R10m\\L2A_T33UUP_20170829T100031_B08_10m.jp2", "raster\\s2a_20170829_b8_raw.tif")# Load rasters
# Pre-disaster images
s2a_20170620_b2 <- raster("raster\\s2a_20170620_b2_raw.tif")
s2a_20170620_b3 <- raster("raster\\s2a_20170620_b3_raw.tif")
s2a_20170620_b4 <- raster("raster\\s2a_20170620_b4_raw.tif")
s2a_20170620_b8 <- raster("raster\\s2a_20170620_b8_raw.tif")
# Post-disaster images
s2a_20170829_b2 <- raster("raster\\s2a_20170829_b2_raw.tif")
s2a_20170829_b3 <- raster("raster\\s2a_20170829_b3_raw.tif")
s2a_20170829_b4 <- raster("raster\\s2a_20170829_b4_raw.tif")
s2a_20170829_b8 <- raster("raster\\s2a_20170829_b8_raw.tif")6. Stacking of bands
A RasterStack represents a collection of RasterLayer objects with the same extent and resolution. Organizing RasterLayer objects in a RasterStack can be practical when dealing with multiple layers; for example to summarize their values or in spatial modeling.
# Stack all the S2 bands in one file
s2a_pre_st <- stack(s2a_20170620_b2, s2a_20170620_b3, s2a_20170620_b4, s2a_20170620_b8)
s2a_post_st <- stack(s2a_20170829_b2, s2a_20170829_b3, s2a_20170829_b4, s2a_20170829_b8)Explore the file.
s2a_pre_st
s2a_post_st
# Number of layers or bands of RasterStack file
nlayers(s2a_pre_st)7. Crop and mask
Since the images are very large in size, they will be cropped to an Area of Interest (AOI). An AOI is stored in a shapefile and can be imported to R. A shapefile is a type of vector.
Import the vector file, change its projection to match the one’s of the rasters and plot them together.
# Import vector file
aoi <- readOGR(paste(w, "\\vector", sep=""), "polygon")
# Get the coordinate reference system of the raster (CRS)
crs(s2a_pre_st)
# Store projection from any raster S2-2A in an object
pr <- CRS("+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
# Alternatively: pr <- crs(s2a_pre_st)
# Reproject: change the projection of the shapefile to equal the one from the rasters
aoi_reprj <- spTransform(aoi, pr)# Mask the image to the extent of the shapefile (product to be masked, shapefile used)
s2a_pre_st_masked <- mask(s2a_pre_st, aoi_reprj)
s2a_post_st_masked <- mask(s2a_post_st, aoi_reprj)
# Crop S2A images to the extent of the new shapefile
s2a_pre_cr <- crop(s2a_pre_st_masked, extent(aoi_reprj))
s2a_post_cr <- crop(s2a_post_st_masked, extent(aoi_reprj)) # compare the extent of the images before and after masking and cropping them
extent(s2a_pre_st)
extent(s2a_pre_cr)8. Visualization
Plot RGB
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.
If a multispectral image consists of the three visual primary colour bands (red, green, blue), the three bands may be combined to produce a “true colour” image. The display colour assignment for any band of a multispectral image can be done in an entirely arbitrary manner. In this case, the colour of a target in the displayed image does not have any resemblance to its actual colour. The resulting product is known as a false colour composite image.
Combine the bands in different order to create a true and false composite of the image. The function plotRGB() can only be used for RasterStack or RasterBrick. In the case of Figure 5, the false color composition helps to distinguish vegetation from other land covers.
Remember to print x11() before the plotting function to see the result in an external window.
# The layout( ) function has the form layout(mat) where mat is a matrix object specifying the location of the N figures to plot.
nf <- layout(matrix(c(1,2, 3, 4), # the data elements
2, # number of rows
2, # number of columns
byrow = TRUE)) # fill matrix by rows? (T/F)
# Indicate the 2 images that will be plotted in the 1X2 matrix
plotRGB(s2a_pre_cr, r=3, g=2, b=1, axes=TRUE, stretch="hist", main = "Sentinel True Color Pre-storm")
plotRGB(s2a_pre_cr, r=4, g=3, b=2, axes=TRUE, stretch="hist", main = "Sentinel False Color Pre-storm")
plotRGB(s2a_post_cr, r=3, g=2, b=1, axes=TRUE, stretch="hist", main = "Sentinel True Color Post-storm")
plotRGB(s2a_post_cr, r=4, g=3, b=2, axes=TRUE, stretch="hist", main = "Sentinel False Color Post-storm")
Figure 5. Subset of Sentinel 2A RasterStack with different color composite
Export
Export RasterStack file. Visualize the raw images and your results in QGIS.
# export raster / overwrite it, if already exists
writeRaster(s2a_pre_cr, datatype="FLT4S", filename = "raster\\s2a_pre_cr.tif", format = "GTiff", overwrite=TRUE, options=c("INTERLEAVE=BAND","COMPRESS=LZW"))
writeRaster(s2a_post_cr, datatype="FLT4S", filename = "raster\\s2a_post_cr.tif", format = "GTiff", overwrite=TRUE, options=c("INTERLEAVE=BAND","COMPRESS=LZW"))9. NDVI: Normalized Difference Vegetation Index
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 modeling. 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.
Figure 6. Explanation of vegetation reflectance. Source: https://support.dronedeploy.com/docs/ndvi-band-order
The formula for calculating the NDVI is the following:
\(NDVI = (NIR-RED)/(NIR+RED)\)
In Sentinel images, NIR is the Band8 and RED is Band 4.
Figure 7. NDVI Formula. Comparison of healthy and stressed vegetation. Source: https://sentinel.esa.int
# Identify with names/labels the bands that you want to work with
s2a_post_cr[[4]] # Which band in this one? Form which RasterStack?
# Rename
red_pre <- s2a_pre_cr[[3]]
nir_pre <- s2a_pre_cr[[4]]
red_post <- s2a_post_cr[[3]]
nir_post <- s2a_post_cr[[4]]
# Calculate the NDVI using the cropped bands 8 and 4
ndvi_pre <- (nir_pre - red_pre)/(nir_pre + red_pre)
ndvi_post <- (nir_post - red_post)/(nir_post + red_post)
# You can also define a function that can be used with any set of objects
fun_ndvi <- function(nir, red) {(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.
Plot
Plot both NDVI results from before and after the wind storm
# The layout( ) function has the form layout(mat) where mat is a matrix object specifying the location of the N figures to plot.
nf2 <- layout(matrix(c(1,2), 1, 2, byrow = TRUE))
# Indicate the 2 images that will be plotted in the 1X2 matrix
plot(ndvi_pre, col = rev(terrain.colors(30)), main = 'NDVI before Wind storm')
plot(ndvi_post, col = rev(terrain.colors(30)), main = 'NDVI after Wind storm')
Figure 8. NDVI plots without and with threshold (0.3)
Export NDVI file in format *.tif. Compare in QGIS
# Export rasters / overwrite it, if already exists
writeRaster(ndvi_pre, datatype="FLT4S", filename = "raster\\ndvi_pre.tif", format = "GTiff", overwrite=TRUE)
writeRaster(ndvi_post, datatype="FLT4S", filename = "raster\\ndvi_post.tif", format = "GTiff", overwrite=TRUE)Thresholding
The non-vegetation pixels (<0.3) can be masked out. The function calc is used to make aritmetic calculations on rasters (among others functions). In this case, a function is created to select all values smaller than 0.3 and make them equal to 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.
# The layout( ) function has the form layout(mat) where mat is a matrix object specifying the location of the N figures to plot.
nf2 <- layout(matrix(c(1,2), 1, 2, byrow = TRUE))
# Indicate the 2 images that will be plotted in the 1X2 matrix
plot(veg_pre , main = 'Vegetation threshold before the Wind storm')
plot(veg_post, main = 'Vegetation threshold after the Wind storm')
Figure 8. NDVI plots without and with threshold (0.3)
Export NDVI file in format *.tif. Compare in QGIS
# Export rasters / overwrite it, if already exists
writeRaster(veg_pre, datatype="FLT4S", filename = "raster\\veg_pre.tif", format = "GTiff", overwrite=TRUE)
writeRaster(veg_post, datatype="FLT4S", filename = "raster\\veg_post.tif", format = "GTiff", overwrite=TRUE)Can you easily detect the changes between healthy vegetation and the damaged one? Why? Explore other thresholds.