Playa

Generate NDVI Rasters from USGS EarthExplorer Landsat 8

Discussion created by Playa on Sep 23, 2017
Latest reply on Sep 23, 2017 by Playa

I have written a Python module that generates NDVI Rasters from the USGS EarthExplorer Landsat 8 with TOA Reflectance and sun angle correction.

 

NDVI Raster: with TOA Reflectance and sun angle correction

 

NDVI Raster: values 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_bands(landsat_dir):
    """
    Create a list of Landsat 8
    tiles bands 4 & 5.
    """

    # Determine how to prevent nested loops - Big 0 Notation
    ndvi_bands = ['_B4.TIF', '_B5.TIF', '_MTL.txt']
    landsat_bands = {}
    p = Path(landsat_dir)
    for directory in p.iterdir():
        for pattern in ndvi_bands:
            try:
                match = '*{0}'.format(pattern)
                landsat_band = directory.glob(match).next()
                landsat_band_name = landsat_band.stem
                landsat_key = re.findall('_(\d{6})_(\d{8})_\d{8}',
                                         str(landsat_band_name))[0]
                landsat_bands.setdefault(landsat_key, []).append(str(landsat_band))
            except (StopIteration, IndexError) as e:
                pattern_name = re.findall('_(\w+)\.', pattern)[0]
                directory_name = str(directory.stem)
                if type(e).__name__ == 'StopIteration':
                    msg = ('Landsat band: {0} not found in directory: {1}.'
                           .format(pattern_name, directory_name))
                    raise StopIteration(msg)
                elif str(type(e).__name__) == 'IndexError':
                    msg = ('Landsat band: {0} has incorrect '
                           'Name (6 digits) or Year (8 digits) format.'
                           .format(landsat_band_name))
                    raise IndexError(msg)
    return landsat_bands


def remove_zero_values(landsat_bands):
    """
    Convert zero cell values
    to NoData.
    """

    arcpy.CheckOutExtension('Spatial')
    for k, v in landsat_bands.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_bands):
    """
    Extract the reflectance
    coefficients from metadata
    txt file.
    """

    for k, v in landsat_bands.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_bands):
    """
    Correct landsat 8
    bands 4 & 5
    for TOA reflectance
    """

    arcpy.CheckOutExtension('Spatial')
    for k, v in landsat_bands.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_bands):
    """
    Correct Landsat 8
    bands 4 & 5
    for sun angle
    """

    arcpy.CheckOutExtension('Spatial')
    for k, v in landsat_bands.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_bands, output_dir):
    """
    Generate NDVI from
    preprocessed
    landsat 8 bands 4 & 5
    """

    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension('Spatial')
    for f, v in landsat_bands.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_bands = list_landsat_bands(landsat_dir)
    print(landsat_bands)
    remove_zero_values(landsat_bands)
    extract_reflectance_coefficients(landsat_bands)
    toa_reflectance_correction(landsat_bands)
    sun_angle_correction(landsat_bands)
    calculate_ndvi(landsat_bands, 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)

 

I would like to verify that my methodology is correct as my output NDVI Rasters value range isn't exactly between 1 and -1 and not sure if its just that its based on the actual tiles values or is every NDVI calculated meant to be exactly between 1 and -1?

Outcomes