Generate NDVI Rasters from USGS EarthExplorer Landsat 8
1707
2
09-23-2017 04:36 AM
by
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?
by
Thanks Dan 
