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
Accepted Solutions
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.
Peter
- 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. - 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.
- 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
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?
Thanks Dan 
