1. Set working directory

# Set the working directory ("wd") and assign it to an object (here, "w")
w <- setwd(".\\your\\path")

# Check the location of the working directory
getwd()

2. Load and install packages

# 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 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(RStoolbox)
library(RColorBrewer)

3. Create change data

For detecting change, you can make classes for the observed change in the landscape. For example:

In QGIS, create at least 10-15 polygons per land class. These polygons correspond to the change data that will be used for the supervised classfication.

Use the Sentinel-2 RasterStack image

Figure 1. Creation of change data.

Figure 1. Creation of change data.

4. Change detection: Classification of rasters

Read raster data

s2_pre_cr <- stack("raster\\s2a_pre_cr.tif")
s2_post_cr <- stack("raster\\s2a_post_cr.tif")

Merge the before and after Sentinel-2 images for the classification

s2_change <- stack(s2_pre_cr, s2_post_cr)
s2_change
## class       : RasterStack 
## dimensions  : 735, 1029, 756315, 8  (nrow, ncol, ncell, nlayers)
## resolution  : 10, 10  (x, y)
## extent      : 391960, 402250, 5387690, 5395040  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : s2a_pre_cr.1, s2a_pre_cr.2, s2a_pre_cr.3, s2a_pre_cr.4, s2a_post_cr.1, s2a_post_cr.2, s2a_post_cr.3, s2a_post_cr.4 
## min values  :            1,          121,           83,          129,             1,             1,             1,           140 
## max values  :        20646,        19098,        17962,        16386,         18292,         17293,         16599,         15701

Read vector data with land cover change classes

In this example, we created a multipolygon shapefile with 5 classes that corresponds to the land covers identifiable on the raster image:

  • 1 : forest_forest
  • 2 : forest_forestdestroyed
  • 3 : agriculture_agriculture
  • 4 : urban_urban
  • 5 : water_water

Import the vector file with these landcover change classes

# Import change samples
change <- readOGR(paste(w, "\\vector", sep=""), "change_data")

Plot the change data on top of the raster subset.

# Plot the raster with both vector files on top
# First, add the raster
plotRGB(s2_pre_cr, r=3, g=2, b=1, axes=TRUE, stretch="lin")

# Second, add the vectors
plot(change, col="yellow", add=TRUE)
Figure 2. Change samples on Sentinel RasterStack

Figure 2. Change samples on Sentinel RasterStack

Classification with superClass() from the RStoolbox Package

windstorm_change <- superClass(img = s2_change,  nSamples=100, trainData = change, responseCol = "class_name")

# Display classes:
windstorm_change$classMapping
##   classID                   class
## 1       1 agriculture_agriculture
## 2       2  forest_destroyedforest
## 3       3           forest_forest
## 4       4             urban_urban
## 5       5             water_water
# Display classification results
windstorm_change
## superClass results
## ************ Validation **************
## $validation
## [1] "No independent validation was performed!"
## 
## *************** Map ******************
## $map
## class       : RasterLayer 
## dimensions  : 735, 1029, 756315  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : 391960, 402250, 5387690, 5395040  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : class_name 
## values      : 1, 5  (min, max)
## attributes  :
##  ID                   value
##   1 agriculture_agriculture
##   2  forest_destroyedforest
##   3           forest_forest
##   4             urban_urban
##   5             water_water

Plot your change detection map

plot(windstorm_change$map, breaks = c(0, 1, 2, 3, 4, 5), col = c("darkolivegreen", "tomato", "yellowgreen", "tan", "blue2") , main = 'Change Detection Classification')
Figure 3. Change detection map. After the windstorm, many areas of the forest where damaged and are show here in red color

Figure 3. Change detection map. After the windstorm, many areas of the forest where damaged and are show here in red color

5. Landscape Statistics

Generate statistics of the area that was classified to quantify the damages

windstorm_change.freq <- freq(windstorm_change$map, useNA= "no")

change.freq <- windstorm_change.freq[, "count"]*10^2*1e-06
barplot(change.freq, main = "Area (km2) of landcover change in our study area",
        col= c("darkolivegreen", "tomato", "yellowgreen", "tan", "blue2"), 
        names.arg = c("Forest", "Forest destroyed", "Agriculture", "Urban", "Water"))
Figure 4. Area per land cover class

Figure 4. Area per land cover class