INTELECTUAL PROPERTY NOTICE: This script is based on the student report from Arthur Lukasewitsch during the course “Rapid Response Techniques in Remote Sensing” on the Winter semester 2018 at the Hochschule Muenchen fuer Angewandte Wissenschaften.

Images were searched using the EO Browser: https://apps.sentinel-hub.com/

Please visit the following link (https://bit.ly/2rl6voT), to have an overview of the products and features from the EO Browser in the selected study area.

For this practice, an index for detection of mud or landslides will be presented. The area of interest (AOI) is Montecito in southwest USA, California. On January 2018 a severe rainfall event led to landslides in the region. For more information on the disaster please go to: https://wapo.st/2RtN2Nz.

Figure 1. Area affected by the disaster, Montecito. California. January 2018. Source: https://wapo.st/2RtN2Nz.

Figure 1. Area affected by the disaster, Montecito. California. January 2018. Source: https://wapo.st/2RtN2Nz.

The downloaded images have the following characteristics: * Sensor: Sentinel 2 * Processing level: L1C * Date pre-disaster: 2017-11-23 * Date post-disaster: 2018-01-22

For this practice, the following indexes and preprocessing tasks are undertaken: * Convertion of images from jp2 to tif format * Selection of bands. Bands to be used: 2, 3, 4, 8 * Crop to AOI * Calculation of NDVI, True color and False color composites for visualization of changes * Calculation of Mud Index

Set up your workspace

# Load libraries

libTOload <- c("sf", "raster", "rgdal", "gdalUtils")
lapply(libTOload, library, character.only = TRUE)

# If necessary, install the packages with the 'install.packages()' function. 
# Set up your working directory
w <- setwd("your\\path\\to_your_data\\here")

getwd()

Load raw S2 images and transform

# Read JP2000 files (to import them directly). When unpacking your S2 data from the Archive, take the folder "L*******" to the first level where the rest of your data will be stored. 
# Images post-disaster
s2a_20180122_b2_raw <- readGDAL("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B02.jp2")
s2a_20180122_b3_raw <- readGDAL("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B03.jp2")
s2a_20180122_b4_raw <- readGDAL("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B04.jp2")
s2a_20180122_b8_raw <- readGDAL("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B08.jp2")
# Transform raster from JP2000 to GeoTIFF

# Post-disaster images
gdal_translate("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B02.jp2", "raster\\processed\\s2_geotiff\\s2a_20180122_b2_raw.tif")
gdal_translate("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B03.jp2", "raster\\processed\\s2_geotiff\\s2a_20180122_b3_raw.tif")
gdal_translate("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B04.jp2", "raster\\processed\\s2_geotiff\\s2a_20180122_b4_raw.tif")
gdal_translate("L1C_T11SKU_A013509_20180122T185440\\IMG_DATA\\T11SKU_20180122T185421_B08.jp2", "raster\\processed\\s2_geotiff\\s2a_20180122_b8_raw.tif")
# Load rasters
# Post-disaster images
s2a_20180122_b2 <- raster("raster\\processed\\s2_geotiff\\s2a_20180122_b2_raw.tif")
s2a_20180122_b3 <- raster("raster\\processed\\s2_geotiff\\s2a_20180122_b3_raw.tif")
s2a_20180122_b4 <- raster("raster\\processed\\s2_geotiff\\s2a_20180122_b4_raw.tif")
s2a_20180122_b8 <- raster("raster\\processed\\s2_geotiff\\s2a_20180122_b8_raw.tif")
# Stack bands in one file for easier management
s2a_post_st <- stack(s2a_20180122_b2, s2a_20180122_b3, s2a_20180122_b4, s2a_20180122_b8)

# Take a look of your files
s2a_post_st
## class      : RasterStack 
## dimensions : 10980, 10980, 120560400, 4  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 199980, 309780, 3790200, 3900000  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : s2a_20180122_b2_raw, s2a_20180122_b3_raw, s2a_20180122_b4_raw, s2a_20180122_b8_raw 
## min values :                   0,                   0,                   0,                   0 
## max values :               65535,               65535,               65535,               65535
writeRaster(s2a_post_st, datatype="FLT4S", filename = "raster\\processed\\s2_geotiff\\s2a_post_st.tif", format = "GTiff", overwrite=TRUE)

Crop to AOI

# Load vector (using package 'sf')
aoi <- st_read("vector\\montecito_aoi.shp")

# Check projection
st_crs(aoi)
# Compare with the projection of your raster
crs(s2a_post_st)
# If different, reassign projection to vector (see previous practices for this using 'sp' package, or "https://r-spatial.github.io/sf/reference/st_crs.html" for explanation using the 'sf' package)


# Crop files to AOI extent
s2a_post_cr <- crop(s2a_post_st, extent(aoi)) 
# Export your results
writeRaster(s2a_post_cr, datatype="FLT4S", filename = "raster\\processed\\s2_geotiff\\s2a_post_cr.tif", format = "GTiff", overwrite=TRUE)

Create the different plots that will allow you to see in true color, false color and with vegetation values the status of the AOI before and after the landslide.

# First calculate NDVI

# Compare s2a_post_cr[4] and s2a_post_cr[[4]]. Which data can you access with each slicing method?
# Which band corresponds to NIR and to RED?
ndvi_post <- (s2a_post_cr[[4]] - s2a_post_cr[[3]])/(s2a_post_cr[[4]] + s2a_post_cr[[3]])
# Export your results
writeRaster(ndvi_post, datatype="FLT4S", filename = "raster\\processed\\indexes\\ndvi_post.tif", format = "GTiff", overwrite=TRUE)
# Plot your results
# Can you tell which band is red, green, nir and blue and in which order they were stored? This is important for the plotting.

par(mfrow=c(2,2),  bty = "n") #'bty' parameter removes the box around the plot for the 'plot()' command
plotRGB(s2a_post_cr, r=3, g=2, b=1, axes = TRUE, margins = TRUE, stretch = 'lin', main = "True color post landslide")

plotRGB(s2a_post_cr, r=4, g=3, b=3, axes = TRUE, margins = TRUE, stretch = 'lin', main = "False color post landslide. Enhacement of Vegetation")

plotRGB(s2a_post_cr, r=2, g=4, b=3, axes = TRUE, margins = TRUE, stretch = 'lin', main = "False color post landslide. Enhacement of bare soil")

plot(ndvi_post, main = 'NDVI after the landslide', las = 1) #'las' parameter changes the orientation of the labels

TIP : orientation of the axis labels: ‘las’ numeric in {0,1,2,3}; the style of axis labels. 0: always parallel to the axis [default], 1: always horizontal, 2: always perpendicular to the axis, 3: always vertical.

Calculate the mud index

mud_index = (s2a_post_cr[[4]]*1.5-s2a_post_cr[[4]]*1.1)/(s2a_post_cr[[4]]*1.5+s2a_post_cr[[4]]*1.1)*(s2a_post_cr[[4]]/0.6)
mudid_post_calc <- calc(mud_index, function(x){x[x < 0.016] <-NA;return(x)})
cuts3 = c(0, 300) 
nrcol = rgb(1, 0, 0, 0.8) 
plot(mudid_post_calc, axes = FALSE, breaks=cuts3, col=nrcol, main="Mud index") 

# Export your results
writeRaster(mudid_post_calc, datatype="FLT4S", filename = "raster\\processed\\indexes\\mudid_post_calc.tif", format = "GTiff", overwrite=TRUE)

Take all the rasters and load them in QGIS. For the Mud index, turn the non values to transparent and check how good is the mask. If the mask is not correct, change the parameters of the mud index. How much can you improve it?