Vegetation Indices (NDVI, MSAVI2 and SAVI) using Landsat 7 on Google Earth Engine - LocateIT/trends.earth GitHub Wiki

vegetation indices LDMS uses the following indices to quantify vegetation cover:

  1. Normalized Difference Vegetation Index (NDVI)
  2. Modified Soil Adjusted Vegetation Index (MSAVI2)
  3. Soil-adjusted Vegetation Index (SAVI)

The following indices are computed using Google Earth Engine. To interact with the code visit here

DATA Description
USGS Landsat 7 Collection 1 Tier 1 and Real-Time data TOA Reflectance ee.ImageCollection("LANDSAT/LE07/C01/T1_RT_TOA")
North Africa countries ee.FeatureCollection("users/derickongeri/Admin")

Import Data

Data is imported and a study period from 1999 to 2020 is defined

// Replace country name with EGYPT, LIBYA, ALGERIA, MAURITANIA, MOROCCO, TUNISIA
var roi = northAfrica.filter(ee.Filter.eq('NAME', 'TUNISIA'))
var countryGeom = northAfrica.first().geometry()

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_RT_TOA')

// Study period
var years = [1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020]

Cloud Masking

Cloud masking applied to remove areas of image experiencing cloud using BQA band bit 4 of Landsat 7

// mask landsat 7 clouds 
function maskL7clouds(image) {
  var qa = image.select('BQA');

  // Bits 4 are clouds.
  var cloudBitMask = 1 << 4;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)

  return image.updateMask(mask);
}

Compute NDVI for Humid Areas

// calculate yearly mean ndvi 
function yearlyNDVIMean(year){
  var datasets = dataset.filterDate(year+'-01-01', year+'-12-31')
    // Pre-filter to get less cloudy granules.
    // .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',50))
    .filterBounds(roi)
    .map(maskL7clouds)
    .mean();
                  
  var ndvi = datasets.expression('(nir- red) / (nir + red)',{
    'nir':datasets.select('B4'),
    'red':datasets.select('B3')
    }
  ).rename('y'+year)
  
  return ndvi
}

Compute SAVI for Desert Areas

SAVI is used to correct Normalized Difference Vegetation Index (NDVI) for the influence of soil brightness in areas where vegetative cover is low. Landsat Surface Reflectance-derived SAVI is calculated as a ratio between the R and NIR values with a soil brightness correction factor (L) defined as 0.5 to accommodate most land cover types.

SAVI = ((NIR - R) / (NIR + R + L)) * (1 + L)

// calculate yearly mean savi
function yearlySAVIMean(year){
  var datasets = dataset.filterDate(year+'-01-01', year+'-12-31')
    // Pre-filter to get less cloudy granules.
    // .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',50))
    .filterBounds(roi)
    .map(maskL7clouds)
    .mean();
                  
  var savi = datasets.expression('((nir - red) / (nir + red + 0.5)) * (1 + 0.5)',{
    'nir':datasets.select('B4'),
    'red':datasets.select('B3')
  }).rename('y'+year)

  
  return savi
}

MSAVI2 Time Series for Arid and Stepic zones

MSAVI2 is soil adjusted vegetation indices that seek to address some of the limitation of NDVI when applied to areas with a high degree of exposed soil surface.It eliminates the need to find the soil line from a feature-space plot or even explicitly specify the soil brightness correction factor:

// calculate year mean msavi2
function yearlyMSAVI2Mean(year){
  var datasets = dataset.filterDate(year+'-01-01', year+'-12-31')
    // Pre-filter to get less cloudy granules.
    // .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',50))
    .filterBounds(roi)
    .map(maskL7clouds)
    .mean();
                  
  
  var msavi2 = datasets.expression('(2 * nir + 1 - ( (2 * nir + 1)**2 - 8 * (nir -red) )**(1/2) ) / 2',{
    'nir':datasets.select('B4'),
    'red':datasets.select('B3')
  }).rename('y'+year)
  
  return msavi2
}

Create Images with bands representing yearly indices

A list of images generated is combined into 1 image:

// create vegetation indices images
var ndviImages = ee.Image.cat(years.map(yearlyNDVIMean))
var msavi2Images = ee.Image.cat(years.map(yearlyMSAVI2Mean))
var saviImages = ee.Image.cat(years.map(yearlySAVIMean))

Visualization Pallete

Create an image visualization based on minimum and maximum values

// indices visualization color ramp 
var ndvi_visualization = {
  min: -0.22789797020331423, 
  max: 0.3912415762453296,
  palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};

var savi_visualization = {
  min: -0.03825978758662746, 
  max: 0.23875073954453754,
  palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};

var msavi2_visualization = {
  min: -0.03159281238913536, 
  max: 0.2486189603805542,
  palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};

Map Layers

Map.addLayer(ndviImages.select('y2019').clip(roi), ndvi_visualization, 'NDVI_y2019')
Map.addLayer(saviImages.select('y2019').clip(roi), savi_visualization, 'SAVI_y2019')
Map.addLayer(msavi2Images.select('y2019').clip(roi), msavi2_visualization, 'MSAVI_y2019')

Export data

Final Image layers can be exported to assets, google drive, google cloud storage or download locally based on google cloud uri provided. Below is an example of downloading locally based on google cloud uri provided

print(ndviImages.clip(roi).getDownloadURL({
  name:'NDVI',
  bands:[{
    id:'y2019',
    scale:300,
  }],
  crs:'EPSG:4326',
  // crs_transform:[1,0, 0, 0, 1, 0],
  region: countryGeom 
}))

For google drive

Export.image.toDrive({image: ndviImages.clip(roi), scale: 30, description: id, fileNamePrefix: 'NDVI',
    region: countryGeom , maxPixels: 1e10});

For gee assets

Export.image.toAsset({
  image:ndviImages.clip(roi),
  maxPixels:1e12
})

For google cloud storage

Export.image.toCloudStorage({
  image:ndviImages.clip(roi),
  maxPixels:1e12,
  scale:30,
  bucket:'ldms',
  region:countryGeom 
}

NB: A higher spatial resolution(30m) warrants a higher file storage size while a lower spatial resolution results in a smaller file storage size.