Krdytkyu

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.

NDVI: Raster Output

Landsat 8: Input Directory

Landsat 8 Tile: Directory

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)