Deriving Spectral Indices with Google Earth Engine in R

R rgee Remote Sensing Google Earth Engine

In this post, I have derived several spectral indices from Sentinel-2 data using Google Earth Engine in R with package rgee.

Abhishek Kumar https://akumar.netlify.app/ (Panjab University, Chandigarh)https://puchd.ac.in/
10-16-2021
Show code 
## set global options for markdown
knitr::opts_chunk$set(comment = "#>", collapse = TRUE, fig.align = 'center',
          fig.width = 7, fig.height = 5, out.width = '90%', dev = "svglite",
          echo=T, message=FALSE, warning=FALSE)

Introduction

Google Earth Engine (GEE) is a cloud-based platform for planetary-scale geo-spatial analysis (Gorelick et al. 2017). The cloud-based computation ability overcomes the limitations of data access, data storage, personal limited computation ability and time. Thus, it provides a platform for large scale analysis of earth observations through access to satellite imagery.

Google offered support only for Python and JavaScript, and GEE was not accessible to R (R Core Team 2021) users until Aybar et al. (2020) developed rgee package (Aybar 2021).

Show code 
library(rgee)
ee_Initialize()

#knitr::write_bib(x = c(.packages(), "blogdown"), "packages.bib")

Data and Pre-processing

Data Product: Sentinel-2 MSI: MultiSpectral Instrument, Level-2A

Sentinel-2 is a wide-swath, high-resolution, multi-spectral imaging mission supporting Copernicus Land Monitoring studies, including the monitoring of vegetation, soil and water cover, as well as observation of inland waterways and coastal areas. The SENTINEL-2 data contain 13 spectral bands representing TOA reflectance scaled by 10000.

Show code 
# define a region of interest
roi <- ee$Geometry$Rectangle(76.85, 30.66, 77.01, 30.75)

# load and pre-process sentinel-2 data
# Filter the image collection
img <- ee$ImageCollection("COPERNICUS/S2_SR")$
  
  # filter the date range we are interested in
  filterDate("2020-01-01", "2020-12-31")$
  
  # filter geographic area we are interested in 
  filterBounds(roi)$
  
  # filter cloud cover less than 20%
  filter(ee$Filter$lt("CLOUD_COVERAGE_ASSESSMENT", 20))$
  
  # filter cloud pixel percentage to less than 20%
  filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", 20))$
  
  # reduce image collection by taking median value
  median()$
  
  # clip to region of interest
  clip(roi)

# print band names
img$bandNames()$getInfo()

#>  [1] "B1"         "B2"         "B3"         "B4"         "B5"        
#>  [6] "B6"         "B7"         "B8"         "B8A"        "B9"        
#> [11] "B11"        "B12"        "AOT"        "WVP"        "SCL"       
#> [16] "TCI_R"      "TCI_G"      "TCI_B"      "MSK_CLDPRB" "MSK_SNWPRB"
#> [21] "QA10"       "QA20"       "QA60"

Median true color composite

Show code 
Map$centerObject(roi, 12)
Map$addLayer(img, 
             visParams = list(min = 0, max = 3000, bands = c("B4", "B3", "B2")),
             name = "true color")

Median false color composite

Show code 
Map$centerObject(roi, 12)
Map$addLayer(img, 
             visParams = list(min = 0, max = 3000, bands = c("B8", "B4", "B3")),
             name = "false color")

Spectral Indices

Vegetation cover

NDVI: Normalized Difference Vegetation Index

Show code 
# calculate indices e.g. NDVI
ndvi <- img$expression(
  expression = "(NIR - RED)/(NIR + RED)",
  opt_map = list("RED" = img$select("B4"),
                 "NIR" = img$select("B8")
                 ))$rename("ndvi")

vis_pal <- c('FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', 
             '74A901', '66A000', '529400', '3E8601', '207401', '056201',
             '004C00', '023B01', '012E01', '011D01', '011301')

Map$centerObject(roi, 12)
Map$addLayer(ndvi, 
             visParams = list(min = 0, max = 1, palette = vis_pal),
             name = "NDVI")

NDVIre: Normalized Difference Vegetation Index red-edge

Show code 
ndvire <- img$expression(
  expression = "(RedEdge1 - RED)/(RedEdge1 + RED)",
  opt_map = list("RED" = img$select("B4"),
                 "RedEdge1" = img$select("B5")
                 ))$rename("ndvire")

vis_pal <- c('FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', 
             '74A901', '66A000', '529400', '3E8601', '207401', '056201',
             '004C00', '023B01', '012E01', '011D01', '011301')

Map$centerObject(roi, 12)
Map$addLayer(ndvire, 
             visParams = list(min = 0, max = 1, palette = vis_pal),
             name = "NDVIre")

EVI: Enhanced vegetation Index

Show code 
# calculate EVI
evi <- img$expression(
  expression = "2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 10000))",
  opt_map = list("BLUE" = img$select("B2"),
                 "RED" = img$select("B4"),
                 "NIR" = img$select("B8")
                 ))$rename("evi")

Map$centerObject(roi, 12)
Map$addLayer(evi, 
             visParams = list(palette = vis_pal),
             name = "EVI")

SAVI: Soil-adjusted Vegetation Index

Show code 
savi <- img$expression(
  expression = "((NIR - RED) / (NIR + RED + 0.5)) * 1.5",
  opt_map = list("RED" = img$select("B4")$multiply(0.0001),
                 "NIR" = img$select("B8")$multiply(0.0001)
                 ))$rename("savi")

Map$centerObject(roi, 12)
Map$addLayer(savi, 
             visParams = list(min = 0, max = 0.5, 
                              palette = c("brown", "yellow", "green")),
             name = "SAVI")

Water bodies

NBR/NDWI: Normalized Burn Ratio / Normalized Difference Water Index

Show code 
ndwi <- img$expression(
  expression = "(NIR - SWIR) / (NIR + SWIR)",
  opt_map = list("NIR" = img$select("B8"),
                 "SWIR" = img$select("B12")
                 ))$rename("ndwi")

Map$centerObject(roi, 12)
Map$addLayer(ndwi, 
             visParams = list(min = 0, max = 0.5, 
                              palette = c("blue", "white", "green")),
             name = "NDWI")

MNDWI: Modified Normalized Difference Water Index

Show code 
mndwi <- img$expression(
  expression = "(GREEN - SWIR1) / (GREEN + SWIR1)",
  opt_map = list("GREEN" = img$select("B3"),
                 "SWIR1" = img$select("B11")
                 ))$rename("mndwi")

Map$centerObject(roi, 12)
Map$addLayer(ndwi, 
             visParams = list(min = 0, max = 0.5, 
                              palette = c("blue", "white", "green")),
             name = "MNDWI")

AWEI_sh: Automated Water Extraction Index - shadow

Show code 
aweish <- img$expression(
  expression = "BLUE + 2.5 * GREEN - 1.5 * (NIR + SWIR1) - 0.25 * SWIR2",
  opt_map = list("BLUE" = img$select("B2"),
                 "GREEN" = img$select("B3"),
                 "NIR" = img$select("B8"),
                 "SWIR1" = img$select("B11"),
                 "SWIR2" = img$select("B12")
                 ))$rename("aweish")

Map$centerObject(roi, 12)
Map$addLayer(aweish, 
             visParams = list(min = 0, max = 1, palette = c("white", "blue")),
             name = "AWEIsh")

Built up areas

NDBI: Normalized Difference Built-up Index

Show code 
ndbi <- img$expression(
  expression = "(SWIR1 - NIR) / (SWIR1 + NIR)",
  opt_map = list("NIR" = img$select("B8"),
                 "SWIR1" = img$select("B11")
                 ))$rename("ndbi")

Map$centerObject(roi, 12)
Map$addLayer(ndbi, 
             visParams = list(min = -1, max = 1, palette = c("white", "brown")),
             name = "NDBI")

BUI: Built-up Index

Show code 
bui <- img$expression(
  expression = "((SWIR1 - NIR) / (SWIR1 + NIR)) - ((NIR - RED)/(NIR + RED))",
  opt_map = list("RED" = img$select("B4"),
                 "NIR" = img$select("B8"),
                 "SWIR1" = img$select("B11")
                 ))$rename("bui")


Map$centerObject(roi, 12)
Map$addLayer(bui, 
             visParams = list(min = -1, max = 0.5, palette = c("green", "white", "brown")),
             name = "BUI")

NDTI: Normalized Difference Tillage Index

Show code 
ndti <- img$expression(
  expression = "(SWIR1 - SWIR2) / (SWIR1 + SWIR2)",
  opt_map = list("SWIR1" = img$select("B11"),
                 "SWIR2" = img$select("B12")
                 ))$rename("ndti")

Map$centerObject(roi, 12)
Map$addLayer(ndti, 
             visParams = list(min = 0, max = 0.5, palette = c("brown", "yellow", "green")),
             name = "NDTI")

BSI: Bare Soil Index

Show code 
bsi <- img$expression(
  expression = "(((SWIR1 + RED) - (NIR + BLUE)) / (SWIR1 + RED)) + (NIR + BLUE)",
  opt_map = list("BLUE" = img$select("B2"),
                 "RED" = img$select("B4"),
                 "NIR" = img$select("B8"),
                 "SWIR1" = img$select("B11"),
                 "SWIR2" = img$select("B12")
                 ))$rename("bsi")

Map$centerObject(roi, 12)
Map$addLayer(bsi, 
             visParams = list(min = 1000, max = 4000, palette = c("white", "brown")),
             name = "BSI")
Aybar, Cesar. 2021. rgee: R Bindings for Calling the Earth Engine API. https://CRAN.R-project.org/package=rgee.
Aybar, Cesar, Qiusheng Wu, Lesly Bautista, Roy Yali, and Antony Barja. 2020. “rgee: An R Package for Interacting with Google Earth Engine.” Journal of Open Source Software 5 (51): 2272. https://doi.org/10.21105/joss.02272.
Gorelick, Noel, Matt Hancher, Mike Dixon, Simon Ilyushchenko, David Thau, and Rebecca Moore. 2017. “Google Earth Engine: Planetary-Scale Geospatial Analysis for Everyone” 202 (December): 18–27. https://doi.org/10.1016/j.rse.2017.06.031.
R Core Team. 2021. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

References

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/abhikumar86/abhikumar86.github.io/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

Kumar (2021, Oct. 16). Abhishek Kumar: Deriving Spectral Indices with Google Earth Engine in R. Retrieved from https://abhikumar86.github.io/posts/2021-10-16-rgee-composites/

BibTeX citation

@misc{kumar2021deriving,
  author = {Kumar, Abhishek},
  title = {Abhishek Kumar: Deriving Spectral Indices with Google Earth Engine in R},
  url = {https://abhikumar86.github.io/posts/2021-10-16-rgee-composites/},
  year = {2021}
}
© Copyright 2021 Abhishek Kumar