Calculating NDVI from Landsat 8 USGS Tiles

11057
9
Jump to solution
09-19-2017 09:41 AM
PeterWilson
Occasional Contributor III

I'm busy writing  a Python module to calculate the NDVI (Normalized Difference Vegetation Index) based on the following post: Using Python to calculate NDVI with multiband imagery. According to the USGS background data on Landsat 8 Product, the tiles are provided as Digitial Numbers 16-bit unsigned integer format. The current examples (tutorials and reference material) that I've found to calculate NDVI use Landsat 5 + 7 which is stored as 8-bit unsigned integer format. Can I still use the same formula to calculate NDVI for Landsat 8 that is now stored as 16-bit unsigned integer format. The formula that I'm using is found in: Making Spatial Decisions Using GIS and Remote Sensing: A Workbook 

NDVI = (IR-R)/(IR+R)

ArcGIS: USGS Landsat 8 GeoTIFF (Bands 4 + 5)

ArcGIS: scale range of Landsat 8 Bands 4 + 5

USGS: Using the USGS Landsat 8 Product

USGS EarthExplorer: Landsat 8 Dataset Selection

USGS EarthExplorer: Landsat 8  Download Options

0 Kudos
1 Solution

Accepted Solutions
DanPatterson_Retired
MVP Emeritus

Peter, if you are checking the script for whether it is working, run it manually and compare results.

Your only other alternative is to generate 'images' that have a known patterns of values from which you could calculate NDVI.  This can be done in numpy

ir = [[1, 64, 128], [128, 192, 255], [255, 128, 1]]  # I have left 0's out for now

r = np.swapaxes(ir, 1, 0)

r 
array([[  1, 128, 255],
       [ 64, 192, 128],
       [128, 255,   1]])

np.divide(ir - r, ir + r)  # or whatever
 
array([[ 0.        , -0.33333333, -0.33159269],
       [ 0.33333333,  0.        ,  0.33159269],
       [ 0.33159269, -0.33159269,  0.        ]])

Obviously, the only catch is to use masked arrays to mask 0 values to prevent division by zero errors. (hope I got stuff right... scale the input ranges to suit)

View solution in original post

9 Replies
DanPatterson_Retired
MVP Emeritus

not sure I follow... 16 bit isn't 0-255 https://en.wikipedia.org/wiki/16-bit 8-bit is https://en.wikipedia.org/wiki/8-bit

PeterWilson
Occasional Contributor III

Hi Dan

I found the following on ResearchGate:

Based on the following I presume that I should be able to calculate the NDVI using same formula, its just that the Landsat 8 is being stored in 16-bit and not 8-bit.

0 Kudos
CodyBenkelman
Esri Regular Contributor

Peter

  1. Agreed with Dan, you should be using the full 16 bit range, not scaling to 8 bit.  The simple difference ratio equation you show

    NDVI = (IR-R)/(IR+R)

    IS the equation for NDVI.  It does not have any requirement to be 8 bit (0-255) so you should use the full range of values available in the data.
  2. Second, unless this is purely an exercise for you to learn more Python, why start from scratch?  ArcGIS can immediately calculate NDVI for you - let us know which version you have (ArcGIS Pro 2.0.1?  ArcMap?) and we can advise for how to do it.  Note the NDVI tools within ArcGIS may automatically apply a color ramp, but you'll presumably want the "scientific" version which outputs the raw floating point values.  
  3. NOTE - this will give you pixel values from -1.0 to +1.0, and you'll presumably need to rescale/convert those for VIEWING the data - but if you're attempting any quantitative analysis, you'll want to do that work with the floating point values for best accuracy
PeterWilson
Occasional Contributor III

Hi Cody

Thanks for the response. I'd like to determine the NDVI for the dry and wet season for my study area and then determine the difference between the two. I'm currently still using ArcGIS 10.3.1, unfortunately no ArcGIS Pro. The link that I provided in my initial post is a blog post by ESRI Australia: Using Python to calculate NDVI with multiband imagery, that has a Python script for determing the NDVI. If you are able to assist in the approach I need to following once I've determined the NDVI for say (i.e. 2017 Feb (Wet) and 2017 Oct (Dry)) ,how I can show the difference in NDVI between the wet and dry season. I'd like to be able to code the following in Python as its part of a larger water project that I'm busy with.

0 Kudos
DanPatterson_Retired
MVP Emeritus

a difference is sa.Minus

0 Kudos
PeterWilson
Occasional Contributor III

Hi Dan

I completed writing my Python module to automatically generate NDVI rasters from Landsat 8 tiles (Downloaded from USGS EarthExplorer) with TOA Reflectance and Sun Angle correction.

NDVI: With TOA Reflectance and Sun Angle correction

NDVI: Value range (between 1 and -1)

'''
Created on 23 Sep 2017

Create NDVI Rasters

with TOA Reflectance

and Sun Angle

correction

@author: PeterW
'''
# import site-packages and modules
import re
import argparse
from pathlib import Path
import arcpy
from arcpy.sa import *


def list_landsat_tiles(landsat_dir):
    """
    Create a list of Landsat 8
    tiles bands 4 & 5.
    """
    # Determine how to prevent nested loops - Big 0 Notation
    landsat_tiles = {}
    ndvi_bands = ['_B4.TIF', '_B5.TIF', '_MTL.txt']
    p = Path(landsat_dir)
    for directory in p.iterdir():
        # add error handling to validate its subdirectories
        for band in ndvi_bands:
            match = '*{0}'.format(band)
            landsat_tile = str(directory.glob(match).next())
            landsat_name = str(directory.stem)
            landsat_key = re.findall('_(\d{6})_\d{8}_(\d{8})', landsat_name)[0]
            landsat_tiles.setdefault(landsat_key, []).append(landsat_tile)
    return landsat_tiles


def remove_zero_values(landsat_tiles):
    """
    Convert zero cell values
    to NoData.
    """
    arcpy.CheckOutExtension('Spatial')
    for k, v in landsat_tiles.iteritems():
        red_band = SetNull(v[0], v[0], 'value=0')
        NIR_band = SetNull(v[1], v[1], 'value=0')
        v[0] = red_band
        v[1] = NIR_band
    arcpy.CheckInExtension('Spatial')


def extract_reflectance_coefficients(landsat_tiles):
    """
    Extract the reflectance
    coefficients from metadata
    txt file.
    """
    for k, v in landsat_tiles.iteritems():
        with open(v[2]) as mlt:
            lines = mlt.read().splitlines()
            reflect_mult = float(lines[187].split('=')[1])
            v[2] = reflect_mult
            reflect_add = float(lines[196].split('=')[1])
            v.append(reflect_add)
            sun_elev = float(lines[76].split('=')[1])
            v.append(sun_elev)


def toa_reflectance_correction(landsat_tiles):
    """
    Correct landsat 8
    bands 4 & 5
    for TOA reflectance
    """
    arcpy.CheckOutExtension('Spatial')
    for k, v in landsat_tiles.iteritems():
        reflect4 = (v[2]*v[0])+v[3]
        v[0] = reflect4
        reflect5 = (v[2]*v[1])+v[3]
        v[1] = reflect5
    arcpy.CheckInExtension('Spatial')


def sun_angle_correction(landsat_tiles):
    """
    Correct Landsat 8
    bands 4 & 5
    for sun angle
    """
    arcpy.CheckOutExtension('Spatial')
    for k, v in landsat_tiles.iteritems():
        sun4 = (v[0]/(Sin(v[4])))
        v[0] = sun4
        sun5 = (v[1]/(Sin(v[4])))
        v[1] = sun5
    arcpy.CheckInExtension('Spatial')


def calculate_ndvi(landsat_tiles, output_dir):
    """
    Generate NDVI from
    preprocessed
    landsat 8 bands 4 & 5
    """
    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension('Spatial')
    for f, v in landsat_tiles.iteritems():
        NDVI_name = '_'.join(f)
        arcpy.AddMessage('Processing {0}.tif NDVI'.format(NDVI_name))
        Num = Float(v[1] - v[0])
        Denom = Float(v[1] + v[0])
        NDVI_raster = Divide(Num, Denom)
        NDVI_output = '{0}\\{1}.tif'.format(output_dir, NDVI_name)
        NDVI_raster.save(NDVI_output)
    arcpy.CheckInExtension('Spatial')


def main(landsat_dir, output_dir):
    """
    Determine NDVI for
    each Landsat tile.
    """
    landsat_tiles = list_landsat_tiles(landsat_dir)
    remove_zero_values(landsat_tiles)
    extract_reflectance_coefficients(landsat_tiles)
    toa_reflectance_correction(landsat_tiles)
    sun_angle_correction(landsat_tiles)
    calculate_ndvi(landsat_tiles, output_dir)


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Calculate the NDVI for Landsat 8 tiles')
    parser.add_argument('--landsat_dir', metavar='path', required=True,
                        help='Input Landsat 8 tile directory')
    parser.add_argument('--output_dir', metavar='path', required=True,
                        help='Output NDVI directory')
    args = parser.parse_args()
    main(landsat_dir=args.landsat_dir,
         output_dir=args.output_dir)

The only question that I have is that my output raster is between 1 and -1 but is not represented as 1 to -1 is that because of the actual values within my tile or is there a problem with my calculation?

DanPatterson_Retired
MVP Emeritus

Peter, if you are checking the script for whether it is working, run it manually and compare results.

Your only other alternative is to generate 'images' that have a known patterns of values from which you could calculate NDVI.  This can be done in numpy

ir = [[1, 64, 128], [128, 192, 255], [255, 128, 1]]  # I have left 0's out for now

r = np.swapaxes(ir, 1, 0)

r 
array([[  1, 128, 255],
       [ 64, 192, 128],
       [128, 255,   1]])

np.divide(ir - r, ir + r)  # or whatever
 
array([[ 0.        , -0.33333333, -0.33159269],
       [ 0.33333333,  0.        ,  0.33159269],
       [ 0.33159269, -0.33159269,  0.        ]])

Obviously, the only catch is to use masked arrays to mask 0 values to prevent division by zero errors. (hope I got stuff right... scale the input ranges to suit)

PeterWilson
Occasional Contributor III

Thanks Dan

0 Kudos
DuminduJayasekera
New Contributor III

Peter,

I have two questions: 

1) I am trying to use the above code with LANDSAT data (LANDSAT 8 OLI_TIRS C1 Level-1) to calculate NDVI for a given state in US. I am getting the following error: 

Runtime error
Traceback (most recent call last):
File "<string>", line 4, in <module>
ImportError: No module named pathlib.

I tried installing pathlib but i cannot install due to security reasons in my official laptop.i tried os.path but did not work too. Do not know the reason.

2). Is there a work around for the above code for using LANDSAT ARD data to calculate NDVI values.? 

Any help is appreciated. 

0 Kudos