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

ipak <- function(pkg){
  newpkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(newpkg)) 
    install.packages(newpkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}


packages <- c("rgdal", "raster","sp", "RColorBrewer", "RStoolbox", "sf")
ipak(packages)
## Lade nötiges Paket: rgdal
## Lade nötiges Paket: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Program Files/R/R-4.1.1/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Program Files/R/R-4.1.1/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Program Files/R/R-4.1.1/library/rgdal/proj
## Lade nötiges Paket: raster
## Lade nötiges Paket: RColorBrewer
## Lade nötiges Paket: RStoolbox
## Lade nötiges Paket: sf
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

Another way of installing the packages one by one, is the following:

# 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.

4. Change detection: Classification of rasters

Read raster data

s2_pre_cr <- stack("data\\s2_change_rasters\\juneCrop_32N.tif")
s2_post_cr <- stack("data\\s2_change_rasters\\julyCrop_32N.tif")

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

s2_pre_cr <- resample(s2_pre_cr,s2_post_cr)
s2_change <- stack(s2_pre_cr, s2_post_cr)
s2_change
## class      : RasterStack 
## dimensions : 357, 512, 182784, 8  (nrow, ncol, ncell, nlayers)
## resolution : 19.9, 19.9  (x, y)
## extent     : 346382.5, 356571.3, 5587894, 5594998  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## names      : juneCrop_32N.1, juneCrop_32N.2, juneCrop_32N.3, juneCrop_32N.4, julyCrop_32N.1, julyCrop_32N.2, julyCrop_32N.3, julyCrop_32N.4 
## min values :       77.35617,      192.66094,       39.42276,      316.78022,       34.02091,       83.35627,       40.24007,      108.95861 
## max values :       5711.047,       6187.704,       6608.060,       6988.164,       9506.598,       8726.718,       9487.597,       9001.882

Read vector data with land cover change classes

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

  • 1 : agriculture - agriculture
  • 2 : forest - forest
  • 3 : river - sediment
  • 4 : urban - urban
  • 5 : urban - sediment
  • 6 : river - river

Import the vector file with these landcover change classes

# Import change samples
change <- readOGR("data\\vector\\wk3_change_data.shp")

# Since the projection did not match, I made a reprojection here.
change_32N <- spTransform(change, crs(s2_change))

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=FALSE, stretch="lin")

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

Figure 2. Change samples on Sentinel RasterStack

Classification with superClass() from the RStoolbox Package

flood_change <- superClass(img = s2_change,  nSamples=100, trainData = change_32N, responseCol = "class_name")
# Display classification results
flood_change
## superClass results
## ************ Validation **************
## $validation
## [1] "No independent validation was performed!"
## 
## *************** Map ******************
## $map
## class      : RasterLayer 
## dimensions : 357, 512, 182784  (nrow, ncol, ncell)
## resolution : 19.9, 19.9  (x, y)
## extent     : 346382.5, 356571.3, 5587894, 5594998  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : class_name 
## values     : 1, 6  (min, max)
## attributes :
##        ID       value
##  from:  1   agri_agri
##   to :  6 urban_urban

Since, the validation data was not included, therefore, no test of the model could be performed. Create a validation data, include it in the model and run again the classifier. Share the results of this step with your colleagues.

Plot your change detection map

plot(flood_change$map, breaks = c(0, 1, 2, 3, 4, 5, 6), col = c("darkolivegreen", "tomato", "yellowgreen", "tan", "blue2", "cyan") , main = 'Change Detection Classification')
Figure 3. Change detection map

Figure 3. Change detection map

5. Landscape Statistics

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

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

change.freq <- flood_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", "cyan"), 
        names.arg = c("agriculture", "forest", "river", "sedimented river", "sedimented urban", "urban"))
Figure 4. Area per land cover class

Figure 4. Area per land cover class 

map_classi <- flood_change$map

writeRaster(map_classi, datatype="FLT4S", filename = "data/wk3_results\\change_detection_classification.tif", format = "GTiff", overwrite=TRUE)