4
\$\begingroup\$

I've written the following using Python Dictionaries and Pathlib Module. I'd like to improve the first function: list_landsat_bands.

I've created a list of file patterns to match, which I then use Pathlib iterdir to iterate over each directory. The reason is that I need to group the Landsat bands retrieved from the pattern matching using Pathlib glob. (i.e. key = (tile number, date): values [tile_band4, tile_band5, tile_metadata (txt)]. I'd like to find a better way of generating the initial dictionary eliminating the nested for loops. Any suggestion if recursion could be used to improve the time complexity of the function. Any other suggestions for improving the Python module are also welcome.

enter image description here

NDVI: Raster Output

enter image description here

Landsat 8: Input Directory

enter image description here

Landsat 8 Tile: Directory

enter image description here

NDVI Rasters: Output Directory

'''
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)
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.