This exercise uses Digital Elevation Model from Tandem-X and Landsat images to map a landslide among the many that happened in Nepal after the earthquake in 2015.

News on the event: https://www.nature.com/news/nepal-earthquake-caused-fewer-landslides-than-feared-1.19038

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(rLandsat8)

3. Digital elevation model

The raster files were downloaded from: https://download.geoservice.dlr.de/TDM90/ Characteristics
* Sensor: TanDEM-X DEM
* Spatial resolution: 90m
* Cloud cover < 10%

# Load DLM rasters

DEM <- raster("raster_dlm\\TDM1_DEM__30_N28E084_V01_C\\DEM\\TDM1_DEM__30_N28E084_DEM.tif")

DEM
## class       : RasterLayer 
## dimensions  : 1201, 1201, 1442401  (nrow, ncol, ncell)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : 83.99958, 85.00042, 27.99958, 29.00042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Users\ulloa-to\Documents\Kurse_Seminaer\Rapid_response_tech_RS_natural_disasters\Practice\dlm_landslides\raster_dlm\TDM1_DEM__30_N28E084_V01_C\DEM\TDM1_DEM__30_N28E084_DEM.tif 
## names       : TDM1_DEM__30_N28E084_DEM 
## values      : -32767, 7507.146  (min, max)

Update in the Raster attributes, the parameter Min/Max Values

DEM <- setMinMax(DEM)

DEM # check in attributes the new parameter that is included
## class       : RasterLayer 
## dimensions  : 1201, 1201, 1442401  (nrow, ncol, ncell)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : 83.99958, 85.00042, 27.99958, 29.00042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Users\ulloa-to\Documents\Kurse_Seminaer\Rapid_response_tech_RS_natural_disasters\Practice\dlm_landslides\raster_dlm\TDM1_DEM__30_N28E084_V01_C\DEM\TDM1_DEM__30_N28E084_DEM.tif 
## names       : TDM1_DEM__30_N28E084_DEM 
## values      : -32767, 7923.936  (min, max)

The units of the raster height is in meters. Invalid pixels have a height/value of -32767.0 meters. (Check Product Specification, page 18. More info: https://tandemx-science.dlr.de/pdfs/TD-GS-PS-0021_DEM-Product-Specification_v3.1.pdf)

Elminate the values lower than zero

DEM <- calc(DEM, function(x){x[x < 0] <- NA;return(x)})

Make a histogram to view the distribution of values of your raster

hist(DEM, main="Distribution of elevation values", 
     col= "purple", 
     maxpixels=30000000)

Plot your raster

plot(DEM, main="Digital Elevation Model")

Modify the map

# add a color map with 7 colors
col=terrain.colors(7)

# add breaks to the colormap (6 breaks = 5 segments)
brk <- c(1000, 2000, 3000, 4000, 5000, 6000, 7000)

plot(DEM, col=col, breaks=brk, main="DEM with more breaks")

The region in which the lanslide occurred, is located at a specific height. In QGIS, load the results of your images and identify this height. With the following commands, make a segmentation of your DEM to the height of interest. Basically, reclassify your DEM using the categories you find suitable

# Create vector withh the values you want to categorize your dataset
reclass_df <- c(0, 1000, NA, 
                1000, 2000, NA, 
                2000, 3000, 1,
                3000, 4000, 2,
                4000, 5000, NA,
                5000, 6000, NA,
                6000, Inf, NA)

reclass_df
##  [1]    0 1000   NA 1000 2000   NA 2000 3000    1 3000 4000    2 4000 5000
## [15]   NA 5000 6000   NA 6000  Inf   NA
# Now you force your data into a matrix of rows and columns
reclass_m <- matrix(reclass_df,
                ncol = 3,
                byrow = TRUE)

reclass_m
##      [,1] [,2] [,3]
## [1,]    0 1000   NA
## [2,] 1000 2000   NA
## [3,] 2000 3000    1
## [4,] 3000 4000    2
## [5,] 4000 5000   NA
## [6,] 5000 6000   NA
## [7,] 6000  Inf   NA
# reclassify the raster using the reclass object - reclass_m
DEM_classified <- reclassify(DEM,
                     reclass_m)

# view reclassified data
barplot(DEM_classified,
        main = "Number of pixels in each class", col=col) 
legend(x="bottomright", y=1.5, legend=c("altitude 2000-3000 m.a.s.l.", "altitude 3000-4000 m.a.s.l."), fill=col)

TASK: you don’t like how it looks like? Change the values of the matrix (the third column)

# assign all pixels that equal 0 to NA or no data value
DEM_classified[DEM_classified == 0] <- NA
# plot
plot(DEM_classified,
     col = c("red", "blue", "green", "chocolate", "yellowgreen", "black", "white"), main = "Elevation at which the landslide took place")

To have an idea of how steep the area of the landslide is, you can calculate some topographic information. Many physical processed in the environment are dependent on the properties of terrain. For example, the slope determines the capacity of drainage of an area. This is relevant for the formation of catchments and the ocurrence of landslides.

slope <- terrain(DEM, "slope", unit = "degrees") # the resultant information is in degrees

summary(slope) 
##                slope
## Min.    1.260363e-02
## 1st Qu. 1.743228e+01
## Median  2.640834e+01
## 3rd Qu. 3.529007e+01
## Max.    8.760489e+01
## NA's    4.415100e+04
aspect <- terrain(DEM, "aspect")

aspect 
## class       : RasterLayer 
## dimensions  : 1201, 1201, 1442401  (nrow, ncol, ncell)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : 83.99958, 85.00042, 27.99958, 29.00042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : aspect 
## values      : 3.020345e-05, 6.283183  (min, max)
# Try plotting them and compare with your other plots. How important are these variables for the landslide event?

4. Ancillary data

Load vector files with the location of the reported landslides. You can find the Shapefiles at: https://data.humdata.org/dataset/nepal-landslides-location

# Load vector locations
locations <- readOGR(paste(w, "\\vector\\npl-landslidesp-unosat-150505", sep=""), "npl_landslidesp_unosat_150505")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\ulloa-to\Documents\Kurse_Seminaer\Rapid_response_tech_RS_natural_disasters\Practice\dlm_landslides\vector\npl-landslidesp-unosat-150505", layer: "npl_landslidesp_unosat_150505"
## with 4 features
## It has 10 fields
# Check CRS
crs(locations)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Besides the location of the disaster, it is also possible to integrate socio-geographic data like the administrative borders of the districts. In this case, the information is available at: https://data.humdata.org/dataset/administrative-bounadries-of-nepal

# Load vector locations
adm_borders <- readOGR(paste(w, "\\vector\\npl_admbnda_nd_20181210_shp\\npl_admbnda_nd_20181210_SHP", sep=""), "npl_admbnda_adm1_nd_20181210")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\ulloa-to\Documents\Kurse_Seminaer\Rapid_response_tech_RS_natural_disasters\Practice\dlm_landslides\vector\npl_admbnda_nd_20181210_shp\npl_admbnda_nd_20181210_SHP", layer: "npl_admbnda_adm1_nd_20181210"
## with 7 features
## It has 12 fields
# Check CRS
crs(adm_borders)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Check the vectors. They cover the whole country. You can check the administrative level in the Attribute Table

# Plot borders
blue_saturated <- rgb(0, 0, 255, 0, maxColorValue = 255) # Set the blue color pf the polygon to full saturation. the alpha channel (4th value) is set to "Not transparent"
plot(adm_borders, col= blue_saturated, # here you add the blue color previously defined
     border = "blue", # define color border
     lwd = 3, # define border thickness
     main = "Administrative borders of Nepal with the locations reported of the landslides")

# Plot locations of the landslides
plot(locations, col="red", lwd = 3, add=TRUE)

Plot the DEM with the vector data about the disaster and administrative borders to have an idea of the landscape

# Plot the raster with both vector files on top
# First, add the raster
plot(DEM, main = "Location of landslides events - Nepal", 
     sub = "Administrative borders and Digital Elevation Model") # add here a lower subtitle

# Second, add the vectors
plot(locations, col="red", lwd = 3, add=TRUE)

bluetrans <- rgb(0, 0, 255, 130, maxColorValue = 255) # Set the blue color of the polygon to almost transparent
plot(adm_borders, col= bluetrans, #here you add the new transparent blue
     border = "blue", # define color border
     lwd = 3, # define border thickness
     add=TRUE)

TASK: Also plot the elevation segmentation, slope and aspect with the ancillary data to understand more the relationship of the landslides with the geophysical variables

5. Load Landsat images and calculate NDVI

Now load the Landsat 8 images correspondent to before and after the event. The images where downloaded from the USGS server (https://earthexplorer.usgs.gov/) with the following characteristics:
* Sensor: Landsat 8 OLI/TIRS C1 Level 1
* Date range: 1.10.2014 - 31.10.2015
* Path: 141
* Row: 40
* Cloud cover < 10%

The images selected are:
* Before the event: LC08_L1TP_141040_20141004_20170418_01_T1
* After the event: LC08_L1TP_141040_20151007_20170403_01_T1

The following part you can do it on your own. Basically, plot images of the following steps:
* Load Landsat images
* Make atmospheric correction
* Create a RasterStack with the bands you need for a true color visualization and a vegetation index calculation
* Load the AOI and crop & mask the RasterStack to its boundaries
* Calculate NDVI

par(mfrow=c(1,2))

# Indicate the 2 images that will be plotted in the  1X2 matrix
plotRGB(lspre_st_cropped, r=4, g=2, b=1, axes=TRUE, stretch="hist", main = "Nepal Landslide (before)")
plotRGB(lspost_st_cropped, r=4, g=2, b=1, axes=TRUE, stretch="hist", main = "Nepal Landslide (after)")

Take into account that for calculating the NDVI with Landsat 8, another bandset should be used. In this case, band 5 = NIR and band 4 = RED.

par(mfrow=c(1,2))

# Set up the limits of the values scale (y axis)
breakpoints <- seq(-1, 1, by=0.2)

plot(ndvi_pre, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI Pre-Landslide')

plot(ndvi_post, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI Post-Landslide')

# Change in NDVI
ndvi_change <- ndvi_post - ndvi_pre

plot(ndvi_change, col = colorRampPalette(c("red", "yellow", "lightgreen"))(length(breakpoints)-1), breaks=breakpoints, main = 'NDVI change - Landslide')

Finally, in order to be able to compare the data, change the DEM and Vector files to the same projection of your Landsat 8 image

For reprojecting a Raster, use:

# Reproject
DEM <- projectRaster(DEM,  crs = "+proj=utm +zone=45 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0")

# Check
proj4string(DEM)
## [1] "+proj=utm +zone=45 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

TASK: Crop the resultant files to the AOI and plot them all together.
TASK: Plot the DLM + ancillary data + optical data and compare