Calculating different Vegetation Indices in Google Earth Engine (Sentinel -2 images)

Calculating different Vegetation Indices in Google Earth Engine (Sentinel -2 images)

In previous post, I have published the article for EVI and NDVI calculation from sentinel 2 image on Google earth engine (GEE) platform. In this article, I will be showing how various vegetation indices can be computed on GEE platforms and can be added to image collection. So, let's dive into various vegatation indices directly. 

Greenness (Green Leaf Index)

The Green Leaf Index, also known as the Greenness Index, shows how the reflectance in the green channel compares to the other two visible wavelengths (red and blue).

The formula for calculation of GLI is:

 GLI = (2*GRN-RED-BLU)/(RED+GRN+BLU)

  • GRN pixel values from the GREEN band
  • RED pixel values from the RED band
  • BLU pixel values from the BLUE band

The lower the greenness index, the higher the green reflection relative to the other wavelengths. A result of -1 indicates that the green channel is fully reflected, but the remaining two wavelengths (red and blue) are not.

# GLI (Green Leaf Index)
def getGLI(image):
    # Compute the GLI using an expression.
    GLI = image.expression ('(((GREEN - RED) + (GREEN - BLUE)) / ((2 * GREEN) + RED + BLUE))', {
        'GREEN': image.select ('B3'),  
        'RED': image.select ('B4'),
        'BLUE': image.select ('B2')
        }).rename("GLI")
    image = image.addBands(GLI)

    return(image)

 

Leaf Area Index (LAI)

The leaf area index (LAI), equal to half of the total green leaf area per unit of the horizontal ground surface, is a crucial structural feature of vegetation. Canopy interception, evapotranspiration, and gross photosynthesis are directly related to LAI since leaf surfaces constitute the primary frontier of energy and mass exchange.

This index is used to forecast crop growth and yield and estimate foliage cover. ENVI calculates green LAI using the empirical formula developed by Boegh et al. (2002):

LAI = (3.618*EVI - 0.118) 

  • EVI is Enhanced Vegetation Index value. Details for calculation of EVI is here in my previous post

LAI values in the high range typically vary from 0 to 3.5. The LAI values can approach 3.5 when the scene comprises clouds and other bright objects that cause saturated pixels. Before making an LAI image, you should mask clouds and bright features from your scene.

Soil-Adjusted Vegetation Index (SAVI):

In places with little vegetative cover, SAVI is used to compensate for the Normalized Difference Vegetation Index (NDVI) due to the influence of soil brightness. Using a soil-brightness adjustment factor, SAVI seeks to reduce soil brightness impacts. It produces values between -1.0 to 1.0. SAVI is commonly used in arid locations with little vegetative cover. SAVI obtained from surface reflectance is calculated as a ratio of R and NIR values with a soil brightness correction factor (L) of 0.5 to account for most land cover types.

SAVI = ((NIR - Red) / (NIR + Red + L)) x (1 + L)

  • NIR = pixel values from the NEAR INFRARED band
  • Red = pixel values from the RED band
  • L = amount of green vegetation cover
# SAVI (Soil Adjusted Vegetation Index)
def getSAVI(image):
    # Compute the SAVI using an expression.
    SAVI = image.expression ('(((NIR - RED) / (NIR + RED + L))* (1+L))',{
        'L': 0.5, # Cobertura vegetation 0-1
        'NIR': image.select ('B8'),
        'RED': image.select ('B4')
    }).rename("SAVI")
    
    image = image.addBands(SAVI)

    return(image)

Green Chlorophyll Index (GCI )

Based on near-infrared and green bands, the GCI (Green Chlorophyll Index) quantifies leaf chlorophyll concentration in plants. The chlorophyll value, in general, reflects the vegetation.

 the formula for the GCI:

GCI = NIR / GRN – 1

  • GRN = pixel values from the green band
  • NIR = pixel values from the near infrared band
# GCI (Green Chlorophyll Index)
def getGCI(image):
    # Compute the GCI using an expression.
    GCI = image.expression ('(((NIR) / (GREEN)) - 1)', {
        'NIR': image.select ('B8'),  
        'GREEN': image.select ('B3')
        }).rename("GCI")
    
    image = image.addBands(GCI)

    return(image)

 

SIPI (Structure Insensitive Pigment Index)

The SIPI index maximizes sensitivity to the bulk carotenoids to chlorophyll ratio while reducing the impact of canopy structural variation. It's especially effective in locations where the canopy structure, or leaf area index, is highly variable.

SIPI readings vary from 0 to 2, with healthy green vegetation falling between 0.8 and 1.8.

# SIPI (Structure Insensitive Pigment Index)
def getSIPI(image):
    # Compute the SIPI using an expression.
    SIPI = image.expression ('((NIR - BLUE) / (NIR - RED))',{
        'NIR': image.select ('B8'),
        'BLUE': image.select ('B2'), 
        'RED': image.select ('B4')
    }).rename("SIPI")
    
    image = image.addBands(SIPI)

    return(image)

ARVI (Atmospherically Resistant Vegetation Index)

ARVI is especially useful in areas with high levels of atmospheric aerosols. It corrects for air scattering phenomena that affect red light reflectance using blue light reflectance data.

 

The formula in general:

ARVI = (NIR - RED - y * (RED - BLUE))/ (NIR + RED - y*(RED-BLUE))

  • NIR = pixel values from the NEAR INFRARED band
  • RED = pixel values from the RED band
  • BLUE = pixel values from the BLUE band

 

ARVI values vary from -1 to 1, with green vegetation typically falling between 0.20 and 0.80.

# ARVI (Atmospherically Resistant Vegetation Index)
def getARVI(image):
    # Compute the ARVI using an expression.
    ARVI = image.expression ('((NIR - (2 * RED) + BLUE) / (NIR + (2 * RED) + BLUE))',{
        'NIR': image.select ('B8'),
        'BLUE': image.select ('B2'), 
        'RED': image.select ('B4')
    }).rename("ARVI")
    
    image = image.addBands(ARVI)

    return(image)

NBRI (Normalized Burned Ratio Index)

NBR is used to determine burn intensity and identify burned areas. It is traditionally calculated as a ratio between the NIR and SWIR values.

Given that vegetation reflects significantly in the NIR portion of the electromagnetic spectrum. Still, poorly in the SWIR, the NIR and SWIR parts of the electromagnetic spectrum are a powerful combination of bands to utilize for this index. A fire scar that incorporates scarred woody plant and ground, on the other hand, has been demonstrated to reflect more strongly in the SWIR region of the electromagnetic spectrum and beyond.

# NBRI (Normalized Burned Ratio Index)
def getNBRI(image):
    # Compute the NBRI using an expression.
    NBRI = image.expression ('(NIR - SWIR) / (NIR + SWIR)', {
        'NIR': image.select ('B8'),  
        'SWIR': image.select ('B12')
    }).rename("NBRI")
    
    image = image.addBands(NBRI)

    return(image)

 

These indices are created as a functions in python. Now we will move on to implement this in jupyter-lab environment and apply for image collections.

Sentinel_data = ee.ImageCollection('COPERNICUS/S2') \
    .filterDate("2022-03-01","2022-03-31").filterBounds(aoi) \
    .map(getNBRI).median()

In this case i have calculated NBRI (Normalized Burned Ratio Index) but you can change the function as per your required indices or you can add multiple calculation by adding .map(getARVI) in above expression. We now add visulaization parameters and add layers to map as follow:

color = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301']
pallete = {"min":0, "max":1, 'palette':color}

 

# initialize our map
map1 = geemap.Map()
map1.centerObject(aoi, 8)
map1.addLayer(Sentinel_data.clip(aoi).select('NBRI'), pallete, "NBRI")

map1.addLayerControl()
map1

 

Source code: github_link

Please let us know your thoughts on this with following comment section we will definetly reach out to you on your query and feedback.

 

6 Comments

Alemu | typo correction

May 22, 2022, 10:48 a.m.

Dear Admins, First of all, I would like to thank you for giving such a detailed explanation. Then, my comment is on the SAVI code part. On that part of the code, at the return part, you put EVI instead of SAVI. your code: image = image.addBands(EVI) but, you don't have variable EVI there. I believe you wanted to say image.addBands(SAVI). So, you made a typo on the argument. Thank you again.

Krishna | Typo correction (resolved)

May 30, 2022, 4:27 p.m.

Thanks, Alemu Typo corrected.

Nakoa | Resolving different band resolutions before calculating index

Aug. 10, 2022, 10:24 p.m.

Thank you for this great post! I was curious if it is necessary to resample Sentinel 2 bands to have the same resolution before calculating a vegetation index in GEE? This would apply to indices that use bands that are 10m and 20m resolutions such as the NBRI calculation in the above example with Band 8 being 10m and Band 12 being 20m resolution. Thanks for your help!

Krishna | Resolving different band resolutions before calculating index (reply)

Aug. 11, 2022, 4:27 a.m.

GEE can work with different image resolutions, with the results defaulting to the highest. Answer taken from: https://gis.stackexchange.com/questions/392981/is-resampling-necessary-while-calculating-normalized-difference-built-up-index

Rajeev | General comments

March 28, 2023, 3:23 p.m.

Hi, Thank you very much for this detailed demonstration. I have a question similar to the ones that have been asked before. I wonder what would be the resolution of an index when we combine multiple bands at different resolutions. For eg. 10 and 20 (would it be 10 or 20)? I am curious to know from what pixel size the GEE extracts the info while creating a NDVI (lets say) time series using a point shapefile. Thank you very much.

kampus jatim | General Comment

May 29, 2023, 3:01 p.m.

Thanks for sharing this useful information! Hope that you will continue with the kind of stuff you are doing.

Leave a Reply