Overview

Background

In order to measure the impact of a flood, it is necessary to detect the changes in the affected region. After creating a water mask with Sentinel-1 images in SNAP, we can compare them with Optical images. Depending of the thresholds used, there might be areas classified as water which are not water. These are called “False positives”, or areas that were wrongly classified as the target category.

A way to detect the water bodies in Optical images, is to calculate the Normalized Difference Water Index (NDWI) using the green and Near Infra-red bands. The water body has strong absorption and low radiation in the range from visible to infrared wavelengths. The NDWI can enhance the water information effectively in most cases. However, it is sensitive to built-up land and often results in over-estimated water bodies.

The NDWI formula is as follows:

\(NDWI = (Green-NIR)/(Green+NIR)\)

To overcome the shortcomings of NDWI, Xu (2006) developed the Modified Normalized Difference Water Index (MNDWI) that uses the Shortwave Infrared (SWIR) band to replace the NIR band used in NDWI. Many authors have demonstrated that MNDWI is more suitable to enhance water information and can extract water bodies with greater accuracy than NDWI

The MNDWI formula is as follows:

\(MNDWI = (Green-SWIR)/(Green+SWIR)\)

In general, compared to NDWI, water bodies have greater positive values in MNDWI, because water bodies generally absorb more SWIR light than NIR light; soil, vegetation and built-up classes have smaller negative values, because they reflect more SWIR light than green light (Du et al. 2016)

In the case of the Sentinel-2 images, the bands to be used for the NDWI are:

\(NDWI (Sentinel-2) = (B3-B8)/(B3+B8)\)

For Sentinel-2, the green band has the spatial resolution of 10m, while the SWIR band (Band 11) has the spatial resolution of 20 m. Thus, MNDWI needs to be calculated at a spatial resolution of either 10 m or 20 m. The 20-m MNDWI is calculated as:

\(MNDWI (Sentinel-2) = (B3-B11)/(B3+B11)\)

Objectives

On the following exercise, we calculate the NDWI and MDWI of a Sentinel-2 image from the same study area, Gorakhpur, India, in the dry season (21.05.2017).The results will be compared with the water mask from the Sentinel-1 images.

News on the event can be found here: https://www.theguardian.com/world/2017/aug/30/mumbai-paralysed-by-floods-as-india-and-region-hit-by-worst-monsoon-rains-in-years

Data

For this practice, use the following files:

  • Raster: s2a_20170521_b3_cr.tif, s2a_20170521_b8_cr.tif, s2a_20170521_b11_cr.tif, s2a_20171023_b3_cr.tif, s2a_20171023_b8_cr.tif, s2a_20171023_b11_cr.tif

Procedure

1. Install and/or load necessary packages

As for every new session of R, you should set your working directory and load the necessary libraries:

# Install and load the necessary packages 
ipak <- function(pkg){
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    if (length(new.pkg)) 
        install.packages(new.pkg, dependencies = TRUE, repos = "http://cran.us.r-project.org")
    sapply(pkg, require, character.only = TRUE)
}

packages <- c("maptools", "rgdal", "raster", "rgeos", "rasterVis", "RCurl", "devtools", "gdalUtils")
ipak(packages)
##  maptools     rgdal    raster     rgeos rasterVis     RCurl  devtools gdalUtils 
##      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE

2. Load rasters

Set up your 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()

When copying the path, take into account the R spelling:
* Either change  for /
* Or use copy the path and run readClipboard() in your console

For this exercise, Sentinel-2 images before and after the event will be used: 21.05.2017 and 23.10.2017, respectively. The images are already converted into GeoTIFF and cropped to the AOI.

# Read raster files
# Pre - flood
s2a_20170521_b3 <- raster("raster\\s2a_20170521_b3_cr.tif")
s2a_20170521_b8 <- raster("raster\\s2a_20170521_b8_cr.tif")
s2a_20170521_b11 <- raster("raster\\s2a_20170521_b11_cr.tif")

# Post - flood
s2a_20171023_b3 <- raster("raster\\s2a_20171023_b3_cr.tif")
s2a_20171023_b8 <- raster("raster\\s2a_20171023_b8_cr.tif")
s2a_20171023_b11 <- raster("raster\\s2a_20171023_b11_cr.tif")

To look at the properties and characteristics of the rasters, make a summary or just enter the name of the object. Check many images and pay attention to the spatial resolution.

# Summary
summary(s2a_20170521_b3)
##         s2a_20170521_b3_cr
## Min.                  1090
## 1st Qu.               1630
## Median                1776
## 3rd Qu.               1912
## Max.                  2548
## NA's                     0
# Print the name of the object
s2a_20170521_b3
## class      : RasterLayer 
## dimensions : 5837, 5869, 34257353  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 701830, 760520, 2912700, 2971070  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs 
## source     : W:/Ulloa/SHK/Lehre/flood_sentinel/gorakhpur/r_script/raster/s2a_20170521_b3_cr.tif 
## names      : s2a_20170521_b3_cr 
## values     : 1069, 4475  (min, max)

3. Stack and plot

Make a RasterStack with the 10m resolution images

# pre-disaster
s2a_20170521_10mst <- stack(s2a_20170521_b3, s2a_20170521_b8)

# post-disaster
s2a_20171023_10mst <- stack(s2a_20171023_b3, s2a_20171023_b8)

Check the new object and properties inserting the name in the Console.

To acces to each band of the new RasterStack, use [[]]

# Band Green (B3)
green_pre <- s2a_20170521_10mst[[1]]
green_post <- s2a_20171023_10mst[[1]]

# Band NIR (B8)
NIR_pre <- s2a_20170521_10mst[[2]]
NIR_post <- s2a_20171023_10mst[[2]]

Visualize each band

# 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), 2, 2, byrow = TRUE))  

# Indicate the 4 images that will be plotted in the  1X2 matrix
plot(green_pre, main="Green Band Pre-flood")
plot(NIR_pre, main="NIR Band Pre-flood")

plot(green_post, main="Green Band Post-flood")
plot(NIR_post, main="NIR Band Post-flood")
Figure 1. Bands Green and NIR from Sentinel-2, Gorakhpur - India.

Figure 1. Bands Green and NIR from Sentinel-2, Gorakhpur - India.

4. NDWI and MNDWI calculations

Now calculate the index, using a simple band arithmetics.

ndwi_pre = (green_pre-NIR_pre)/(green_pre+NIR_pre)
ndwi_post = (green_post-NIR_post)/(green_post+NIR_post)

McFeeters in Du et al. (2016) indicates that values of NDWI greater than zero are assumed to represent water surfaces, while values less than, or equal, to zero are assumed to be non-water surfaces.

# Threshold: Mask out the values smaller than 0 which are supposed to be non-water surfaces
water_ndwi_pre <- calc(ndwi_pre, function(x){x[x < 0] <- NA;return(x)})
water_ndwi_post <- calc(ndwi_post, function(x){x[x < 0] <- NA;return(x)})

Plot all the results

# With 'par()' you can indicate in how many columns and rows the images have to be arranged
par(mfrow=c(2,3)) 

# Indicate the 6 images that will be plotted in the  2X3 matrix. Pay attention to the order in which the images are printed! They correspond to the values of the 2x3 plotting matrix
plot(ndwi_pre, main="NDWI Pre-flood")
plot(water_ndwi_pre, main="Water pixels NDWI Pre-flood")

cuts=c(-1, 0, 1)
pal= colorRampPalette(c("black", "white"))
plot(ndwi_pre, breaks=cuts, col=pal(2), main= "Water mask NDWI Pre-flood")

plot(ndwi_post, main="NDWI Post-flood")
plot(water_ndwi_post, main="Water pixels NDWI Post-flood")

cuts=c(-1, 0, 1)
pal= colorRampPalette(c("black", "white"))
plot(ndwi_post, breaks=cuts, col=pal(2), main= "Water mask NDWI Post-flood")
Figure 2. NDWI calculations, Gorakhpur - India.

Figure 2. NDWI calculations, Gorakhpur - India.

5. Resample

In order to calculate MNDWI, both bands (B3 & B11) have to have the same resolution.

# Check the resolution of bands
res(s2a_20170521_b3)
## [1] 10 10
res(s2a_20170521_b11)
## [1] 20 20

In order to match both resolutions, B3 has to be resampled from 10x10m to 20x20m.

Resample from 10m to 20m resolution

The change in the spatial resolution of a raster image is known as image resampling. There are two common methods of resampling:

  • “ngb”: Nearest-neighbor; assigns the value of the nearest cell.
  • “bilinear”: Bilinear interpolation; assigns a weighted average of the four nearest cells (the default).

Figure 3. Summary of different types of resampling methods.

# Resample band 3 to the spatial resolution of band 11
s2a_20170521_b3_res <- resample(s2a_20170521_b3, s2a_20170521_b11, method="bilinear") 
s2a_20171023_b3_res <- resample(s2a_20171023_b3, s2a_20171023_b11, method="bilinear")

Check the resolution of the resulting raster

res(s2a_20170521_b3_res)
## [1] 20 20
res(s2a_20171023_b3_res)
## [1] 20 20

Now, stack and calculate the MNDWI with the new B3 band and band 11

# Stack 
s2a_20170521_20mst <- stack(s2a_20170521_b3_res, s2a_20170521_b11)
s2a_20171023_20mst <- stack(s2a_20171023_b3_res, s2a_20171023_b11)

# Band Green (B3) - 20m
green_pre_r <- s2a_20170521_20mst[[1]]
green_post_r <- s2a_20171023_20mst[[1]]

# Band NIR (B8)
SWIR_pre <- s2a_20170521_20mst[[2]]
SWIR_post <- s2a_20171023_20mst[[2]]

# Calculate MNDWI
mndwi_pre = (green_pre_r-SWIR_pre)/(green_pre_r+SWIR_pre)
mndwi_post = (green_post_r-SWIR_post)/(green_post_r+SWIR_post)

# Threshold: Mask out the values smaller than 0 which are supposed to be non-water surfaces
water_mndwi_pre <- calc(mndwi_pre, function(x){x[x < 0] <- NA;return(x)})
water_mndwi_post <- calc(mndwi_post, function(x){x[x < 0] <- NA;return(x)})

Plot the results to visualize changes

# With 'par()' you can indicate in how many columns and rows the images have to be arranged
par(mfrow=c(2,3)) 

# Indicate the 6 images that will be plotted in the  2X3 matrix
plot(mndwi_pre, main="MNDWI Pre-flood")
plot(water_mndwi_pre, main="Water pixels MNDWI, Pre-flood")

cuts=c(-1, 0, 1)
pal= colorRampPalette(c("black", "white"))
plot(mndwi_pre, breaks=cuts, col=pal(2), main= "Water mask MNDWI, Pre-flood")

plot(mndwi_post, main="MNDWI, Post-flood")
plot(water_mndwi_post, main="Water pixels MNDWI, Post-flood")

cuts=c(-1, 0, 1)
pal= colorRampPalette(c("black", "white"))
plot(mndwi_post, breaks=cuts, col=pal(2), main= "Water mask MNDWI, Post-flood")
Figure 4. MNDWI calculations, Gorakhpur - India.

Figure 4. MNDWI calculations, Gorakhpur - India.

6. Export all the relevant images

# Export raster/Overwrite it, if already exists
writeRaster(water_ndwi_pre, datatype="FLT4S", filename = "raster\\water_ndwi_pre.tif", format = "GTiff", overwrite=TRUE)
writeRaster(water_mndwi_pre, datatype="FLT4S", filename = "raster\\water_mndwi_pre.tif", format = "GTiff", overwrite=TRUE)

writeRaster(water_ndwi_post, datatype="FLT4S", filename = "raster\\water_ndwi_post.tif", format = "GTiff", overwrite=TRUE)
writeRaster(water_mndwi_post, datatype="FLT4S", filename = "raster\\water_mndwi_post.tif", format = "GTiff", overwrite=TRUE)

Export the raster files and afterwards visualize them in QGIS Compare them with the unprocessed raw Sentinel-2 images. Load also your SNAP Processed Sentinel-1 images as well as the Kennaugh calculations and compare.

Discuss about the differences in methods and sensing dates of the images before and after the disaster.

Figure 5. Summary of all results, methods and images done so far for water bodies detection in the frame of flood events.

References

Du, Yu., Yihang Zhang, Feng Ling, Qunming Wang, Wenbo Li, Xiaodong Li. 2016. Water Bodies’ Mapping from Sentinel-2 Imagery with Modified Normalized Difference Water Index at 10-m Spatial Resolution Produced by Sharpening the SWIR Band. Remote Sensing, 8(4), 354; https://doi.org/10.3390/rs8040354

Xu, Hanqiu. 2006. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14): 3025-3033