#! /usr/bin/env python

'''
    FILE: estimate_landsat_emissivity.py

    PURPOSE: Estimates a Landsat Emissivity product from ASTER Emissivity and
             NDVI.  The results are meant to be used for generation of a
             Surface Temperature product.

    PROJECT: Land Satellites Data Systems (LSDS) Science Research and
             Development (LSRD) at the USGS EROS

    LICENSE: NASA Open Source Agreement 1.3
'''

import os
import sys
import math
import json
import logging
from collections import namedtuple

import numpy as np
import astropy.convolution.convolve as ap_convolve
from astropy.convolution import Box2DKernel

from osgeo import gdal, osr

from espa import Metadata, Sys, Geo, Dataset
from st_exceptions import NoTilesError, InaccessibleTileError, MissingBandError

# Import local modules
import st_utilities as util
import emissivity_utilities as emis_util
import camel_utilities as cm_util

logger = logging.getLogger(__name__)
ASTER_EMISSIVITY_MEAN_SCALE_FACTOR = 0.001
NDVI_SCALE_FACTOR = 0.01
PQA_SNOW = 5                      # Bit 5 in the Level 1 BQA Pixel Band
PQA_WATER = 7                     # Bit 7 in the Level 1 BQA Pixel Band
PQA_SINGLE_BIT = 0x01             # 00000001

CoefficientInfo = namedtuple('CoefficientInfo',
                             ('estimated_1', 'estimated_2', 'estimated_3',
                              'snow_emissivity', 'water_emissivity',
                              'nominal_emissivity', 'vegetation_coeff',
                              'camel_1', 'camel_2', 'camel_3', 'max_ndvi',
                              'min_ndvi', 'bare_soil_coeff'))


def sensor_coefficients(satellite, band, multiple_bands, st_data_dir):
    """Determines the sensor specific coefficients

    Args:
        satellite <str>: Satellite we are currently processing
        band <str>: Band we are currently processing
        multiple_bands <bool>: more than one band implies Split Window
        st_data_dir <str>: Location of the ST data files

    Returns:
        <CoefficientInfo>: Populated coefficients
    """

    # Read coefficient file.
    coefficient_file = os.path.join(st_data_dir, 'coefficients.json')
    with open(coefficient_file) as coefficient_data:
        all_coeff = json.load(coefficient_data)
        coefficient_data.close()

    if satellite == 'LANDSAT_4': # From JPL
        coeff = all_coeff['Single_Channel']['L4']
    elif satellite == 'LANDSAT_5': # From JPL
        coeff = all_coeff['Single_Channel']['L5']
    elif satellite == 'LANDSAT_7': # From JPL
        coeff = all_coeff['Single_Channel']['L7']
    elif satellite == 'LANDSAT_8':
        if multiple_bands: # Split Window, from RIT
            if band == 'band10':
                coeff = all_coeff['Split_Window']['L8_B10']
            elif band == 'band11':
                coeff = all_coeff['Split_Window']['L8_B11']
            else:
                raise Exception('Unsupported band')
        else: # Single Channel, from JPL
            if band == 'band10':
                coeff = all_coeff['Single_Channel']['L8_B10']
            elif band == 'band11': # not used in Single Channel
                coeff = all_coeff['Single_Channel']['L8_B11']
            else:
                raise Exception('Unsupported band')
    elif satellite == 'LANDSAT_9':
        if multiple_bands: # Split Window, from RIT
            if band == 'band10':
                coeff = all_coeff['Split_Window']['L9_B10']
            elif band == 'band11':
                coeff = all_coeff['Split_Window']['L9_B11']
            else:
                raise Exception('Unsupported band')
        else: # Single Channel, from JPL
            if band == 'band10':
                coeff = all_coeff['Single_Channel']['L9_B10']
            elif band == 'band11': # not used in Single Channel
                coeff = all_coeff['Single_Channel']['L9_B11']
            else:
                raise Exception('Unsupported band')
    else:
        raise Exception('Unsupported satellite')

    return CoefficientInfo(
        max_ndvi=all_coeff['max_ndvi'],
        min_ndvi=all_coeff['min_ndvi'],
        estimated_1=coeff['estimated_1'],
        estimated_2=coeff['estimated_2'],
        estimated_3=coeff['estimated_3'],
        snow_emissivity=coeff['snow_emissivity'],
        water_emissivity=coeff['water_emissivity'],
        nominal_emissivity=coeff['nominal_emissivity'],
        bare_soil_coeff=coeff['bare_soil_coeff'],
        vegetation_coeff=coeff['vegetation_coeff'],
        camel_1=coeff['camel_coeff_1'],
        camel_2=coeff['camel_coeff_2'],
        camel_3=coeff['camel_coeff_3'])


def generate_landsat_ndvi(src_info, no_data_value):
    """Generate Landsat NDVI

    Args:
        src_info <SourceInfo>: Information about the source data
        no_data_value <int>: No data (fill) value to use

    Returns:
        <numpy.2darray>: Generated NDVI band data
    """

    logger.info('Building NDVI band for Landsat data')

    # NIR ----------------------------------------------------------------
    nir_data = Dataset.extract_raster_data(src_info.band_type.nir.name, 1)
    nir_no_data_locations = np.where(nir_data == src_info.band_type.nir.fill_value)
    nir_data = nir_data * src_info.band_type.nir.scale_factor \
        + src_info.band_type.nir.add_offset

    # RED ----------------------------------------------------------------
    red_data = Dataset.extract_raster_data(src_info.band_type.red.name, 1)
    red_no_data_locations = np.where(red_data == src_info.band_type.red.fill_value)
    red_data = red_data * src_info.band_type.red.scale_factor \
        + src_info.band_type.red.add_offset

    # NDVI ---------------------------------------------------------------
    ndvi_data = ((nir_data - red_data) / (nir_data + red_data))

    # Cleanup no data locations
    ndvi_data[nir_no_data_locations] = no_data_value
    ndvi_data[red_no_data_locations] = no_data_value

    # Memory cleanup
    del red_data
    del nir_data
    del nir_no_data_locations
    del red_no_data_locations

    # Turn all negative values to zero
    # Use a really small value so that we don't have negative zero (-0.0)
    ndvi_data[ndvi_data < 0.0000001] = 0

    return ndvi_data.astype(np.float32)


def get_snow_locations(src_info):
    """Generate Landsat snow locations

    Args:
        src_info <SourceInfo>: Information about the source data

    Returns:
        list(<int>): Locations where we decided snow exists
    """

    logger.info('Finding snow locations using pixel QA band')

    # Use the QA_pixel band snow bit to detect snow
    qa_data = Dataset.extract_raster_data(src_info.pixel_qa, 1)
    qa_snow_shifted_mask = np.right_shift(qa_data, PQA_SNOW)
    qa_snow_mask = np.bitwise_and(qa_snow_shifted_mask, PQA_SINGLE_BIT)

    # Save the locations for the snow pixels
    snow_locations = np.where(qa_snow_mask == 1)

    return snow_locations


def extract_aster_data(url, filename, first_band, multiple_bands,
                       keep_temporary):
    """Extracts the internal band(s) data for later processing

    Args:
        url <str>: URL to retrieve the file from
        filename <str>: Base HDF filename to extract from
        first_band <bool>: Is this the first band to be processed?
        multiple_bands <bool>: Is more than one band being processed?
        keep_temporary <bool>: Keep any temporary products generated

    Returns:
        <numpy.2darray>: Mean Band 13 data
        <numpy.2darray>: Mean Band 14 data
        <numpy.2darray>: NDVI band data
        <int>: Samples in the data
        <int>: Lines in the data
        <2x3:float>: GDAL Affine transformation matrix
                     [0] - Map X of upper left corner
                     [1] - Pixel size in X direction
                     [2] - Y rotation
                     [3] - Map Y of upper left corner
                     [4] - X rotation
                     [5] - Pixel size in Y direction
        <bool>: True if the ASTER tile is available, False otherwise
    """

    # Get accessible tile file
    tile_info = emis_util.locate_aster_ged_tile(url=url, filename=filename)
    h5_file_path = tile_info.h5_file_path

    # There are cases where the emissivity data will not be available
    # (for example, in water regions).
    aster_b13_data = []
    aster_b14_data = []
    aster_ndvi_data = []
    samps = 0
    lines = 0
    geo_transform = []
    if not os.path.exists(h5_file_path):
        # The ASTER tile is not available, so don't try to process it
        return (aster_b13_data, aster_b14_data, aster_ndvi_data, samps, lines,
                geo_transform, False)

    # Define the sub-dataset names
    emis_ds_name = ''.join(['HDF5:"', h5_file_path,
                            '"://Emissivity/Mean'])
    ndvi_ds_name = ''.join(['HDF5:"', h5_file_path,
                            '"://NDVI/Mean'])
    lat_ds_name = ''.join(['HDF5:"', h5_file_path,
                           '"://Geolocation/Latitude'])
    lon_ds_name = ''.join(['HDF5:"', h5_file_path,
                           '"://Geolocation/Longitude'])

    logger.debug(emis_ds_name)
    logger.debug(ndvi_ds_name)
    logger.debug(lat_ds_name)
    logger.debug(lon_ds_name)

    aster_b13_data = Dataset.extract_raster_data(emis_ds_name, 4)
    aster_b14_data = Dataset.extract_raster_data(emis_ds_name, 5)
    aster_ndvi_data = Dataset.extract_raster_data(ndvi_ds_name, 1)
    aster_lat_data = Dataset.extract_raster_data(lat_ds_name, 1)
    aster_lon_data = Dataset.extract_raster_data(lon_ds_name, 1)

    # Determine the minimum and maximum latitude and longitude
    x_min = aster_lon_data.min()
    x_max = aster_lon_data.max()
    y_min = aster_lat_data.min()
    y_max = aster_lat_data.max()

    del aster_lon_data
    del aster_lat_data

    # Determine the resolution and dimensions of the ASTER data
    (x_res, y_res, samps, lines) = (
        emis_util.data_resolution_and_size(lat_ds_name,
                                           x_min, x_max, y_min, y_max))

    # Build the geo transform
    geo_transform = [x_min, x_res, 0, y_max, 0, -y_res]

    return (aster_b13_data, aster_b14_data, aster_ndvi_data, samps, lines,
            geo_transform, True)


def generate_estimated_emis_tile(coefficients, tile_name, aster_b13_data,
                                 aster_b14_data, samps, lines, transform, wkt,
                                 no_data_value, no_data_locations):
    """Generate emissivity values for the tile

    Args:
        coefficients <CoefficientInfo>: coefficients for the math
        tile_name <str>: Filename to create for the tile
        aster_b13_data <numpy.2darray>: Unscaled band 13 ASTER tile data
        aster_b14_data <numpy.2darray>: Unscaled band 14 ASTER tile data
        samps <int>: Samples in the data
        lines <int>: Lines in the data
        transform <2x3:float>: GDAL Affine transformation matrix
                               [0] - Map X of upper left corner
                               [1] - Pixel size in X direction
                               [2] - Y rotation
                               [3] - Map Y of upper left corner
                               [4] - X rotation
                               [5] - Pixel size in Y direction
        wkt <str>: Well-Known-Text describing the projection
        no_data_value <float>: Value to use for fill
        no_data_locations [<int>]: Locations to be set to nodata
    """

    # Scale the data
    aster_b13_data = aster_b13_data * ASTER_EMISSIVITY_MEAN_SCALE_FACTOR
    aster_b14_data = aster_b14_data * ASTER_EMISSIVITY_MEAN_SCALE_FACTOR

    # ------------------------------------------------------------
    # Create the estimated Landsat EMIS data
    emis_data = (coefficients.estimated_1 * aster_b13_data +
                 coefficients.estimated_2 * aster_b14_data +
                 coefficients.estimated_3)
    emis_data[no_data_locations] = no_data_value 
    del no_data_locations

    # Create the estimated Landsat EMIS raster output tile
    logger.debug('Creating an estimated Landsat EMIS tile {}'.
                 format(tile_name))
    Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                             tile_name,
                             emis_data,
                             samps, lines,
                             transform,
                             wkt,
                             no_data_value,
                             gdal.GDT_Float32)
    del emis_data


def generate_aster_ndvi_tile(tile_name, ndvi_data, samps, lines, transform,
                             wkt, no_data_value, no_data_locations):
    """Generate ASTER NDVI values for the tile

    Args:
        tile_name <str>: Filename to create for the tile
        ndvi_data <numpy.2darray>: Unscaled NDVI Band ASTER tile data
        samps <int>: Samples in the data
        lines <int>: Lines in the data
        transform <2x3:float>: GDAL Affine transformation matrix
                               [0] - Map X of upper left corner
                               [1] - Pixel size in X direction
                               [2] - Y rotation
                               [3] - Map Y of upper left corner
                               [4] - X rotation
                               [5] - Pixel size in Y direction
        wkt <str>: Well-Known-Text describing the projection
        no_data_value <float>: Value to use for fill
        no_data_locations [<int>]: Locations to be set to nodata
    """

    # Scale the data
    data = ndvi_data * NDVI_SCALE_FACTOR

    # Apply the no data locations, which are populated with CAMEL values
    data[no_data_locations] = no_data_value 

    # Create the ASTER NDVI raster output tile
    logger.debug('Creating an ASTER NDVI tile {}'.format(tile_name))
    Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                             tile_name,
                             data,
                             samps, lines,
                             transform,
                             wkt,
                             no_data_value,
                             gdal.GDT_Float32)
    del data


def generate_tiles(src_info, coefficients, st_data_dir, url, wkt,
                   no_data_value, antimeridian_crossing, first_band,
                   multiple_bands, keep_temporary, camel_filename, band,
                   tirs_only, zenith_angle):
    """Generate tiles for emissivity mean and NDVI from ASTER data

    Args:
        src_info <SourceInfo>: Information about the source data
        coefficients <CoefficientInfo>: coefficients for the math
        st_data_dir <str>: Location of the ST data files
        url <str>: URL to retrieve the file from
        wkt <str>: Well-Known-Text describing the projection
        no_data_value <float>: Value to use for fill
        antimeridian_crossing <boolean>: Flag for scene crossing 180 meridian
        first_band <bool>: Is this the first band to be processed?
        multiple_bands <bool>: Is more than one band being processed?
        keep_temporary <bool>: Keep any temporary products generated
        camel_filename <str>: File with CAMEL tile data
        band <str> : Band to process
        tirs_only <bool>: Is this a TIRS-only scene?
        zenith_angle <float>: Solar zenith angle from metadata

    Returns:
        list(<str>): Mean ASTER GED based emissivity tile names
        list(<str>): Mean ASTER NDVI tile names
        <str>: Mean CAMEL-based emissivity file name
    """

    '''
    Process through the latitude and longitude ASTER tiles which cover
    the Landsat scene we are processing
    - Download them
    - Extract the Emissivity mean bands 13 and 14
    - Extract the NDVI
    - Generate the Landsat EMIS from the 13 and 14 band data
    '''

    # Get list of available ASTER GED tiles and CAMEL equivalents
    (aster_ged_filenames, camel_tile_list) = \
        emis_util.get_aster_ged_tiles_for_src(st_data_dir, src_info,
           antimeridian_crossing)

    # If it's a TIRS-only scene, don't use any ASTER GED tiles that might be
    # available at the location.  TIRS-only doesn't have the Landsat NDVI bands
    # required to adjust ASTER GED tiles to the current scene
    if tirs_only:
        logger.debug('TIRS-only scene, skipping GED')
        aster_ged_filenames = []

    # Don't use ASTER GED if it's around night since NDVI isn't accurate there
    if zenith_angle >= cm_util.NIGHT_THRESHOLD:
        logger.debug('Night or near-night scene with solar zenith angle {0}, skipping GED'.format(zenith_angle))
        aster_ged_filenames = []

    # Initialize filename lists
    aster_emis_mean_filenames = []
    aster_ndvi_mean_filenames = []

    for filename in aster_ged_filenames:
        # Build the output tile names
        root_filename = '.'.join(filename.split('.')[0:5])
        aster_emis_tile_name = ''.join([root_filename, '_emis.tif'])
        aster_ndvi_tile_name = ''.join([root_filename, '_ndvi.tif'])

        # Read the ASTER data
        (aster_b13_data, aster_b14_data, aster_ndvi_data, samps, lines,
         transform, aster_data_available) = (
             extract_aster_data(url=url,
                                filename=filename,
                                first_band=first_band,
                                multiple_bands=multiple_bands,
                                keep_temporary=keep_temporary))

        # Fail if a tile can't be read, but it is in the ASTER GED
        if not aster_data_available:
            raise InaccessibleTileError(
                'Cannot reach tile {} in ASTER GED'.format(filename))

        # Save the no data and low value locations.  Emissivity values under the
        # threshold are unrealistic and likely the result of undetected clouds
        no_data_locations = np.where(((aster_b13_data
                                     * ASTER_EMISSIVITY_MEAN_SCALE_FACTOR)
                                     < emis_util.EMIS_LOW_THRESHOLD) |
                                     (aster_b13_data == no_data_value) |
                                     ((aster_b14_data
                                     * ASTER_EMISSIVITY_MEAN_SCALE_FACTOR)
                                     < emis_util.EMIS_LOW_THRESHOLD) |
                                     (aster_b14_data == no_data_value) |
                                     (aster_ndvi_data == no_data_value))

        # Add the tile names to the list for mosaic building and warping
        aster_emis_mean_filenames.append(aster_emis_tile_name)
        aster_ndvi_mean_filenames.append(aster_ndvi_tile_name)

        generate_estimated_emis_tile(coefficients=coefficients,
                                     tile_name=aster_emis_tile_name,
                                     aster_b13_data=aster_b13_data,
                                     aster_b14_data=aster_b14_data,
                                     samps=samps,
                                     lines=lines,
                                     transform=transform,
                                     wkt=wkt,
                                     no_data_value=no_data_value,
                                     no_data_locations=no_data_locations)

        del aster_b13_data
        del aster_b14_data

        if first_band:
            generate_aster_ndvi_tile(tile_name=aster_ndvi_tile_name,
                                     ndvi_data=aster_ndvi_data,
                                     samps=samps,
                                     lines=lines,
                                     transform=transform,
                                     wkt=wkt,
                                     no_data_value=no_data_value,
                                     no_data_locations=no_data_locations)

        del no_data_locations
        del aster_ndvi_data

    # Build CAMEL file
    camel_emis_mean_filename \
        = cm_util.generate_camel_file(camel_filename, camel_tile_list,
                                      coefficients, wkt, no_data_value,
                                      first_band, "mean", band)

    return (aster_emis_mean_filenames, aster_ndvi_mean_filenames,
            camel_emis_mean_filename)


def build_ls_emis_data(server_name, server_path, st_data_dir, src_info,
                       coefficients, aster_emis_warped_name,
                       aster_ndvi_warped_name, camel_emis_warped_name,
                       no_data_value, keep_temporary, num_threads,
                       use_ias_resampler, band, first_band, multiple_bands,
                       camel_filename, tirs_only, zenith_angle):
    """Build estimated Landsat Emissivity Data

    Args:
        server_name <str>: Name of the ASTER GED server
        server_path <str>: Path on the ASTER GED server
        st_data_dir <str>: Location of the ST data files
        src_info <SourceInfo>: Information about the source data
        coefficients <CoefficientInfo>: coefficients for the math
        aster_emis_warped_name <str>: Path to the warped GED emissivity file
        aster_ndvi_warped_name <str>: Path to the warped ASTER NDVI file
        camel_emis_warped_name <str>: Path to the warped CAMEL emissivity file
        no_data_value <int>: No data (fill) value to use
        keep_temporary <bool>: Keep any temporary products generated
        num_threads <int>: Number of processing threads to use for
                           multithreaded components
        use_ias_resampler <bool>: Do not use gdalwarp for warp_raster,
        band <str> : Band to process
        first_band <bool>: Is this the first band to be processed?
        multiple_bands <bool>: Is more than one band being processed?
        camel_filename <str>: File with CAMEL tile data
        tirs_only <bool>: Is this a TIRS-only scene?
        zenith_angle <float>: Solar zenith angle from metadata

    Returns:
        ged_present <bool>: Is any ASTER GED based data present? 
    """

    # Specify the base URL to use for retrieving the ASTER GED data
    url = ''.join(['http://', server_name, server_path])

    # The ASTER data is in geographic projection so specify that here
    ds_srs = osr.SpatialReference()
    ds_srs.ImportFromEPSG(4326)
    geographic_wkt = ds_srs.ExportToWkt()
    # Save the source wkt string to use during warping
    src_wkt = ds_srs.ExportToWkt()

    # Check for antimeridian crossing
    start_longitude = int(math.floor(src_info.bound.west))
    end_longitude = int(math.floor(src_info.bound.east))

    if start_longitude > 0 and end_longitude < 0:
        antimeridian_crossing = True
    else:
        antimeridian_crossing = False

    (aster_emis_mean_filenames, aster_ndvi_mean_filenames,
        camel_emis_mean_filename) = (
        generate_tiles(src_info=src_info,
                       coefficients=coefficients,
                       st_data_dir=st_data_dir,
                       url=url,
                       wkt=geographic_wkt,
                       no_data_value=no_data_value,
                       antimeridian_crossing=antimeridian_crossing,
                       first_band=first_band,
                       multiple_bands=multiple_bands,
                       keep_temporary=keep_temporary,
                       camel_filename=camel_filename,
                       band=band,
                       tirs_only=tirs_only,
                       zenith_angle=zenith_angle))

    # If there is no ASTER GED data for the scene, skip the mosaic/warp steps
    # for it.
    ged_present = True
    if len(aster_emis_mean_filenames) == 0:
        ged_present = False

    # Define the temporary names
    if multiple_bands:
        aster_emis_mosaic_name = 'aster_emis_' + band + '_mosaic.tif'
        camel_emis_mosaic_name = 'camel_emis_' + band + '_mosaic.tif'
    else:
        aster_emis_mosaic_name = 'aster_emis_mosaic.tif'
        camel_emis_mosaic_name = 'camel_emis_mosaic.tif'

    aster_ndvi_mosaic_name = 'aster_ndvi_mosaic.tif'

    # If the image crosses the 180 meridian, shift the tile longitudes to use
    # the 0..360 range so the mosaic is not confused
    if antimeridian_crossing:
        emis_util.shift_tiles(aster_emis_mean_filenames)

        if first_band:
            emis_util.shift_tiles(aster_ndvi_mean_filenames)

    # Mosaic the estimated Landsat EMIS tiles into the temp EMIS.  The CAMEL
    # mosaic isn't needed since it's 1 file from the beginning now.  Just
    # rename the file
    logger.info('Building mosaics for emissivity for band {0}'.format(band))
    if ged_present:
        Geo.mosaic_tiles_into_one_raster(aster_emis_mean_filenames,
                                         aster_emis_mosaic_name,
                                         no_data_value)
    os.rename(camel_emis_mean_filename, camel_emis_mosaic_name)

    if first_band:
        # Mosaic the ASTER NDVI tiles into the temp NDVI
        if ged_present:
            logger.info('Building mosaic for ASTER NDVI')
            Geo.mosaic_tiles_into_one_raster(aster_ndvi_mean_filenames,
                                             aster_ndvi_mosaic_name,
                                             no_data_value)

    if not keep_temporary:
        # Cleanup the estimated ASTER EMIS tiles
        for emis_filename in aster_emis_mean_filenames:
            if os.path.exists(emis_filename):
                os.unlink(emis_filename)

        # Cleanup the ASTER NDVI tiles
        for ndvi_filename in aster_ndvi_mean_filenames:
            if os.path.exists(ndvi_filename):
                os.unlink(ndvi_filename)

    # Warp estimated Landsat EMIS to match the Landsat data
    logger.info('Warping estimated Landsat EMIS to match Landsat data')
    if use_ias_resampler:
        if ged_present:
            emis_util.create_omf(aster_emis_mosaic_name,
                                 src_info.band_type.thermal.name,
                                 num_threads)
            emis_odl_filename = emis_util.create_odl(aster_emis_warped_name,
                                                 no_data_value)
            emis_util.warp_raster(emis_odl_filename)
        emis_util.create_omf(camel_emis_mosaic_name,
                             src_info.band_type.thermal.name,
                             num_threads)
        emis_odl_filename = emis_util.create_odl(camel_emis_warped_name,
                                                 no_data_value)
        emis_util.warp_raster(emis_odl_filename)
    else:
        if ged_present:
            Geo.warp_raster_using_gdalwarp(src_info, src_wkt, no_data_value,
                                           aster_emis_mosaic_name,
                                           aster_emis_warped_name, 1, 'cubic')
        # Give pixel size for this one so TIRS-only cases don't use red band
        Geo.warp_raster_using_gdalwarp(src_info, src_wkt, no_data_value,
                                       camel_emis_mosaic_name,
                                       camel_emis_warped_name, 1, 'cubic',
                                       x_pixel_size=src_info.band_type.thermal.
                                           pixel_size.x,
                                       y_pixel_size=src_info.band_type.thermal.
                                           pixel_size.y)

    if first_band and ged_present:
        logger.info('Warping ASTER NDVI to match Landsat data')
        if use_ias_resampler:
            emis_util.create_omf(aster_ndvi_mosaic_name,
                                 src_info.band_type.thermal.name,
                                 num_threads)
            aster_odl_filename = emis_util.create_odl(aster_ndvi_warped_name,
                                                      no_data_value)
            emis_util.warp_raster(aster_odl_filename)
        else:
            Geo.warp_raster_using_gdalwarp(src_info, src_wkt, no_data_value,
                                           aster_ndvi_mosaic_name,
                                           aster_ndvi_warped_name, 1,
                                           'cubic')

    if not keep_temporary:
        # Cleanup the temp files
        if os.path.exists(aster_emis_mosaic_name):
            os.unlink(aster_emis_mosaic_name)
        if os.path.exists(camel_emis_mosaic_name):
            os.unlink(camel_emis_mosaic_name)
        if os.path.exists(aster_ndvi_mosaic_name):
            os.unlink(aster_ndvi_mosaic_name)

    return ged_present


def extract_warped_emis_data(aster_emis_warped_name, camel_emis_warped_name,
                             no_data_value, keep_temporary, ged_present):
    """Retrieves the warped emissivity image data

    Args:
        aster_emis_warped_name <str>: Path to the warped ASTER GED emissivity file
        camel_emis_warped_name <str>: Path to the warped CAMEL emissivity file
        no_data_value <float>: Value to use for fill
        keep_temporary <bool>: Keep any temporary products generated
        ged_present <bool>: Is any ASTER GED based data present?

    Returns:
        <numpy.2darray>: ASTER GED emissivity data
        list(<int>): Emissivity locations containing no data (fill) values
        <numpy.2darray>: CAMEL emissivity data
    """

    # Load the warped emissivity source values into memory
    camel_emis_data = Dataset.extract_raster_data(camel_emis_warped_name, 1)

    if ged_present:
        # Load the warped ASTER EMIS into memory
        aster_emis_data = Dataset.extract_raster_data(aster_emis_warped_name, 1)
    else:
        # Initialize ASTER structures to be empty
        aster_emis_data = np.full_like(camel_emis_data, no_data_value)

    # Define nodata locations
    aster_emis_no_data_locations = np.where((aster_emis_data == 0) |
                                         (aster_emis_data == no_data_value))

    if not keep_temporary:
        # Cleanup the temporary files since we have them in memory
        if os.path.exists(aster_emis_warped_name):
            os.unlink(aster_emis_warped_name)
        if os.path.exists(camel_emis_warped_name):
            os.unlink(camel_emis_warped_name)

    # ASTER variables are "ls" on the other side because they will become
    # estimated Landsat structures later
    return (aster_emis_data, aster_emis_no_data_locations, camel_emis_data)


def extract_warped_ndvi_data(aster_ndvi_warped_name, camel_emis_data,
                             no_data_value, ged_present):
                             
    """Retrieve and adjust the warped NDVI data

    Args:
        aster_ndvi_warped_name <str>: Path to the warped ASTER NDVI file
        camel_emis_data <numpy.2darray>: CAMEL emissivity data
        no_data_value <float>: Value to use for fill
        ged_present <bool>: Is any ASTER GED based data present?

    Returns:
        <numpy.2darray>: ASTER NDVI data
    """

    if ged_present:
        # Load the warped ASTER NDVI into memory
        aster_ndvi_data = Dataset.extract_raster_data(aster_ndvi_warped_name, 1)
    else:
        # Initialize ASTER NDVI structure to be empty
        aster_ndvi_data = np.full_like(camel_emis_data, no_data_value)

    # Turn all negative values to zero
    # Use a really small value so that we don't have negative zero (-0.0)
    aster_ndvi_data[aster_ndvi_data < 0.0000001] = 0

    return aster_ndvi_data


def update_water_locations(src_info, ls_emis_final, emis_source_data,
                           samps, lines, transform, wkt, no_data_value,
                           water_emissivity, first_band):
    """
       Update locations identified as water to use the water emissivity
       value if they would otherwise be nodata

    Args:
        src_info <SourceInfo>: Information about the source data
        ls_emis_final <raster>: 2D raster array data
        emis_source_data <raster>: 2D raster array data
        samps <int>: Samples in the data
        lines <int>: Lines in the data
        transform <2x3:float>: GDAL Affine transformation matrix
                               [0] - Map X of upper left corner
                               [1] - Pixel size in X direction
                               [2] - Y rotation
                               [3] - Map Y of upper left corner
                               [4] - X rotation
                               [5] - Pixel size in Y direction
        wkt <str>: Well-Known-Text describing the projection
        no_data_value <int>: No data (fill) value to use
        water_emissivity <float>: Constant value for water emissivity
        first_band <bool>: Is this the first band to be processed?

    Returns:
        <numpy.2darray>: Landsat emissivity
        <numpy.2darray>: Emissivity source
    """
    qa_data = Dataset.extract_raster_data(src_info.pixel_qa, 1)
    qa_water_shifted_mask = np.right_shift(qa_data, PQA_WATER)
    qa_water_mask = np.bitwise_and(qa_water_shifted_mask, PQA_SINGLE_BIT)

    # Replace all QA_PIXEL water locations with predefined water emissivity
    # value.  This reduces artifacts at ASTER GED/CAMEL boundaries
    water_update_locations = np.where(qa_water_mask == 1)
    ls_emis_final[water_update_locations] = water_emissivity

    if first_band:  # Only write out once

        # Update emissivity source where water was updated
        emis_source_data[water_update_locations] = cm_util.WATER_EMIS_SOURCE

    return (ls_emis_final, emis_source_data)


def emis_size(emis_data, ls_emis_no_data_locations, thermal_no_data_locations,
              emis_type, no_data_value):
    """Determines the size of the valid data for the type (ASTER GED or CAMEL).
       The input is for the full scene bounding box, but we only want to count
       the parts in the valid Landsat scene region, not the Landsat fill region.
       For CAMEL, we only want to count the part that would be used, which is
       the ASTER GED locations that are missing.

    Args:
        emis_data <numpy.2darray>: emissivity data (ASTER GED or CAMEL)
        ls_emis_no_data_locations list(<int>): Emissivity locations containing
                                               no data (fill) values
        thermal_no_data_locations list(<int>): Thermal band locations containing
                                               no data(fill) values
        emis_type <str>: Type of data to check ("GED" or "CAMEL")
        no_data_value <int>: No data (fill) value to use

    Returns:
        <int>: Emissivity valid data size
    """

    emis_with_nodata = np.copy(emis_data)
    emis_with_nodata[thermal_no_data_locations] = no_data_value
    if emis_type == "GED":
        valid_size = np.count_nonzero(emis_with_nodata != no_data_value)
    else: # "CAMEL"
        valid_size = np.count_nonzero(
            emis_with_nodata[ls_emis_no_data_locations] != no_data_value)

    return valid_size


def generate_emissivity_data(xml_filename, server_name, server_path,
                             camel_path, st_data_dir, no_data_value,
                             keep_temporary, num_threads, use_ias_resampler,
                             band_list):
    """Provides the main processing algorithm for generating the estimated
       Landsat emissivity product.  It produces the final emissivity product.

    Args:
        xml_filename <str>: Filename for the ESPA Metadata XML
        server_name <str>: Name of the ASTER GED server
        server_path <str>: Path on the ASTER GED server
        camel_path <str>: Path to the CAMEL archive
        st_data_dir <str>: Location of the ST data files
        no_data_value <int>: No data (fill) value to use
        keep_temporary <bool>: Keep any temporary products generated
        num_threads <int>: Number of processing threads to use for
                           multithreaded components
        use_ias_resampler <bool>: Do not use gdalwarp for warp_raster,
                                  use the ias makegeomgrid/geomresampler
        band_list <str>: List of bands to process
    """

    # XML metadata
    espa_metadata = Metadata(xml_filename)
    espa_metadata.parse()
    camel_filename = cm_util.get_camel_filename(espa_metadata, camel_path,
                                                "mean")
    if camel_filename == "CAMEL_file_not_found":
        raise NoTilesError('No matching month of CAMEL data found')
    product_id = espa_metadata.xml_object.global_metadata.product_id.text
    sensor_code = util.get_satellite_sensor_code(product_id)
    tirs_only = False
    if sensor_code in ['LT08', 'LT09']:
        tirs_only = True
    src_info = emis_util.retrieve_metadata_information(espa_metadata, tirs_only)

    # Look up solar zenith angle in metadata
    solar_angles = espa_metadata.xml_object.global_metadata.solar_angles
    zenith_angle = float(solar_angles.get('zenith'))

    # Determine output information
    dataset = gdal.Open(src_info.band_type.thermal.name)
    if dataset is None:
        raise MissingBandError('Missing Level 1 Thermal Band')
    output_srs = osr.SpatialReference()
    output_srs.ImportFromWkt(dataset.GetProjection())
    output_transform = dataset.GetGeoTransform()
    samps = dataset.RasterXSize
    lines = dataset.RasterYSize
    del dataset

    # Read thermal data to identify fill locations
    thermal_data = Dataset.extract_raster_data(
        src_info.band_type.thermal.name, 1)
    thermal_no_data_locations = np.where(thermal_data ==
        src_info.band_type.thermal.fill_value)

    # Build NDVI in memory
    if tirs_only or (zenith_angle >= cm_util.NIGHT_THRESHOLD):
        # TIRS-only has no Landsat NDVI data, and Landsat NDVI isn't accurate
        # at/near night, so make the variable, but empty in those cases
        ls_ndvi_data = np.zeros_like(thermal_data)
    else:
        ls_ndvi_data = generate_landsat_ndvi(src_info, no_data_value)

        if keep_temporary:
            logger.info('Writing Landsat NDVI raster')
            # Write the Landsat NDVI raster
            Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                     'internal_landsat_ndvi.tif',
                                     ls_ndvi_data,
                                     samps,
                                     lines,
                                     output_transform,
                                     output_srs.ExportToWkt(),
                                     no_data_value,
                                     gdal.GDT_Float32)

        # Replace LS NDVI values greater than 1 with 1
        ls_ndvi_data[ls_ndvi_data > 1.0] = 1

    # Determine Snow locations
    snow_locations = get_snow_locations(src_info)

    if len(band_list) > 1:
        multiple_bands = True
    else:
        multiple_bands = False

    logger.info('Multiple bands: ' + str(multiple_bands))
    first_band = True

    for band in band_list:
        logger.info('Processing band: ' + str(band))
        # Initialize coefficients.
        coefficients = sensor_coefficients(espa_metadata.xml_object
                                           .global_metadata.satellite,
                                           band, multiple_bands,
                                           st_data_dir=st_data_dir) 
        logger.debug('Emissivity coefficients for band {0}: '
                     '{1}'.format(str(band), coefficients)) 
                     
        if multiple_bands:
            aster_emis_warped_name = 'aster_emis_' + band + '_warped.tif'
            camel_emis_warped_name = 'camel_emis_' + band + '_warped.tif'
        else:
            aster_emis_warped_name = 'aster_emis_warped.tif'
            camel_emis_warped_name = 'camel_emis_warped.tif'

        aster_ndvi_warped_name = 'aster_ndvi_warped.tif'

        # Build the estimated Landsat EMIS data from the ASTER GED data and
        # warp it to the Landsat scenes projection and image extents
        # For convenience the ASTER NDVI is also extracted and warped to the
        # Landsat scenes projection and image extents
        logger.info('Build thermal emissivity band and retrieve ASTER NDVI')
        ged_present = build_ls_emis_data(server_name=server_name,
                           server_path=server_path,
                           st_data_dir=st_data_dir,
                           src_info=src_info,
                           coefficients=coefficients,
                           aster_emis_warped_name=aster_emis_warped_name,
                           aster_ndvi_warped_name=aster_ndvi_warped_name,
                           camel_emis_warped_name=camel_emis_warped_name,
                           no_data_value=no_data_value,
                           keep_temporary=keep_temporary,
                           num_threads=num_threads,
                           use_ias_resampler=use_ias_resampler,
                           band=band,
                           first_band=first_band,
                           multiple_bands=multiple_bands,
                           camel_filename=camel_filename,
                           tirs_only=tirs_only,
                           zenith_angle=zenith_angle)

        (ls_emis_data, ls_emis_no_data_locations, camel_emis_data) = \
           (extract_warped_emis_data(aster_emis_warped_name=aster_emis_warped_name,
                                     camel_emis_warped_name=camel_emis_warped_name,
                                     no_data_value=no_data_value,
                                     keep_temporary=keep_temporary,
                                     ged_present=ged_present))
        if first_band:
            aster_ndvi_data = \
               (extract_warped_ndvi_data(aster_ndvi_warped_name=aster_ndvi_warped_name,
                                         camel_emis_data=camel_emis_data,
                                         no_data_value=no_data_value,
                                         ged_present=ged_present))

        if first_band:
            ged_size = emis_size(ls_emis_data, ls_emis_no_data_locations,
                                 thermal_no_data_locations, "GED", no_data_value)
            camel_size = emis_size(camel_emis_data, ls_emis_no_data_locations,
                                   thermal_no_data_locations, "CAMEL",
                                   no_data_value)
            ged_fraction = ged_size / (ged_size + camel_size)
            logger.debug('Valid ASTER GED size {0}, valid CAMEL size at ASTER '
                         'GED gaps {1}, ASTER GED fraction {2}, GED threshold '
                         '{3}'.format(ged_size, camel_size, ged_fraction,
                         cm_util.MIX_THRESHOLD))

            if ged_fraction < cm_util.MIX_THRESHOLD:
                skip_decision = "SKIP"
                with open(cm_util.SKIP_GED_FILE, 'w') as skip_ged_fd:
                    skip_ged_fd.write("SKIP")
            else:
                skip_decision = "NO_SKIP"
                with open(cm_util.SKIP_GED_FILE, 'w') as skip_ged_fd:
                    skip_ged_fd.write("NO_SKIP")
        else: # Base the threshold decision for the 2nd band on the first band
            skip_ged_file = cm_util.SKIP_GED_FILE
            skip_decision = "NO_SKIP"
            if os.path.exists(skip_ged_file):
                with open(skip_ged_file, 'r') as skip_ged_fd:
                    skip_decision = skip_ged_fd.readline()
                skip_ged_fd.close()

        if skip_decision == "SKIP":
            # Replace ASTER GED with CAMEL
            ls_emis_data = np.copy(camel_emis_data)
            aster_ndvi_data = np.full_like(ls_emis_data,
                cm_util.CAMEL_DEFAULT_NDVI)
            if first_band:
                emis_source_data = np.full_like(ls_emis_data,
                    cm_util.CAMEL_EMIS_SOURCE)
        else:
            # Fill gaps in ASTER GED with CAMEL
            ls_emis_data[ls_emis_no_data_locations] \
                = camel_emis_data[ls_emis_no_data_locations]
            aster_ndvi_data[ls_emis_no_data_locations] \
                = cm_util.CAMEL_DEFAULT_NDVI
            if first_band:
                emis_source_data = np.zeros_like(ls_emis_data)
                emis_source_data[ls_emis_no_data_locations] \
                    = cm_util.CAMEL_EMIS_SOURCE

        if keep_temporary:
            logger.info('Writing gap_filled raster')
            if multiple_bands:
                gap_filled_emis = 'gap_filled_emis_' + band + '.tif'
            else:
                gap_filled_emis = 'gap_filled_emis.tif'
            Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                     gap_filled_emis,
                                     ls_emis_data,
                                     samps,
                                     lines,
                                     output_transform,
                                     output_srs.ExportToWkt(),
                                     no_data_value,
                                     gdal.GDT_Float32)
        if keep_temporary and first_band:
            Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                     'gap_filled_ndvi.tif',
                                     aster_ndvi_data,
                                     samps,
                                     lines,
                                     output_transform,
                                     output_srs.ExportToWkt(),
                                     no_data_value,
                                     gdal.GDT_Float32)
            Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                     'gap_filled_emis_source.tif',
                                     emis_source_data,
                                     samps,
                                     lines,
                                     output_transform,
                                     output_srs.ExportToWkt(),
                                     util.EMIS_SRC_NO_DATA_VALUE,
                                     gdal.GDT_Byte)

        if first_band:
            # Replace NDVI values greater than 1 with 1
            aster_ndvi_data[aster_ndvi_data > 1.0] = 1

            # Calculate fractional vegetation cover for Landsat
            fv_ls = 1 - ((coefficients.max_ndvi - ls_ndvi_data) / \
                         (coefficients.max_ndvi - coefficients.min_ndvi))

            # Calculate fractional vegetation cover for ASTER
            fv_aster = 1 - ((coefficients.max_ndvi - aster_ndvi_data) /  \
                         (coefficients.max_ndvi - coefficients.min_ndvi))

        # Soil - From prototype code variable name
        logger.info('Calculating bare soil component')

        # Calculate bare soil fraction from ASTER.  The 0.975 is based
        # on ASTER spectral response for bands 13/14 for vegetation
        ls_emis_bare = ((ls_emis_data - 0.975 * fv_aster) / \
                        (1 - fv_aster))

        # Account for instability in (1-fv_aster) denominator when fv_aster is
        # large by fixing bare component to spectral library emissivity of soil
        ls_emis_bare[fv_aster > coefficients.max_ndvi] = \
            coefficients.bare_soil_coeff

        # Smooth bare soil emissivity from ASTER (to minimize artifacts arising
        # from threshold above)
        if keep_temporary and first_band:
            logger.info('Writing non-smoothed ls_emis_bare raster')
            Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                     'ls_emis_bare.tif',
                                     ls_emis_bare,
                                     samps,
                                     lines,
                                     output_transform,
                                     output_srs.ExportToWkt(),
                                     no_data_value,
                                     gdal.GDT_Float32)

        # Build kernel for convolution
        kernel = Box2DKernel(util.KERNEL_SIZE)
        ls_emis_bare = ap_convolve(ls_emis_bare, kernel,
                                   nan_treatment='interpolate',
                                   normalize_kernel=True, boundary='extend')

        if keep_temporary and first_band:
            logger.info('Writing smoothed ls_emis_bare raster')
            Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                     'smoothed_ls_emis_bare.tif',
                                     ls_emis_bare,
                                     samps,
                                     lines,
                                     output_transform,
                                     output_srs.ExportToWkt(),
                                     no_data_value,
                                     gdal.GDT_Float32)

        # Calculate veg adjustment with Landsat
        logger.info('Calculating EMIS Final')

        # Final emissivity adjustment by Landsat NDVI
        logger.info('Adjusting estimated EMIS for vegetation')
        ls_emis_final = coefficients.vegetation_coeff * fv_ls + ls_emis_bare \
            * (1.0 - fv_ls)

        # Memory cleanup
        del ls_emis_bare

        # The above adjustments aren't for CAMEL pixels, which don't have
        # NDVI values
        if skip_decision == "SKIP":
            ls_emis_final = np.copy(camel_emis_data)
        else:
            ls_emis_final[ls_emis_no_data_locations] \
                = camel_emis_data[ls_emis_no_data_locations]

        # Apply default snow unless it's a 100% CAMEL scene.  CAMEL V003 is
        # accurate enough to not need default snow.  Handle snow consistently
        # in mixed GED/CAMEL scenes.
        if skip_decision != "SKIP":
            # Medium snow
            logger.info('Adjusting estimated EMIS for snow')
            if first_band:
                emis_source_data[snow_locations] = cm_util.SNOW_EMIS_SOURCE
            ls_emis_final[snow_locations] = coefficients.snow_emissivity

            if keep_temporary:
                logger.info('Writing ls_emis_final raster after snow update')
                if multiple_bands:
                    emis_final_filename = 'emis_final_after_snow_' + band + '.tif'
                else:
                    emis_final_filename = 'emis_final_after_snow.tif'
                Geo.generate_raster_file(gdal.GetDriverByName('GTiff'),
                                         emis_final_filename,
                                         ls_emis_final,
                                         samps,
                                         lines,
                                         output_transform,
                                         output_srs.ExportToWkt(),
                                         no_data_value,
                                         gdal.GDT_Float32)

        # Final check for emissivity values greater than 1. Reset values greater
        # than 1 to nominal value (should be very few, if any)
        ls_emis_final[np.where(ls_emis_final > 1.0)] = \
            coefficients.nominal_emissivity

        # Reset water values
        if first_band:
            emis_source_data[np.where(ls_emis_data > coefficients.water_emissivity)] = \
                cm_util.WATER_EMIS_SOURCE
        ls_emis_final[np.where(ls_emis_data > coefficients.water_emissivity)] = \
            coefficients.water_emissivity

        # Memory cleanup
        del ls_emis_data

        # Set any emissivity values less than a threshold to CAMEL. This filters
        # out bad ASTER GED emissivity values usually due to undetected cloud
        # that results in emissivity underestimation
        low_locations = np.where((ls_emis_final < emis_util.EMIS_LOW_THRESHOLD) &
                                 (ls_emis_final != no_data_value))
        ls_emis_final[low_locations] = camel_emis_data[low_locations] 
        if first_band:
            emis_source_data[low_locations] = cm_util.CAMEL_EMIS_SOURCE

        # Memory cleanup
        del camel_emis_data

        # Check threshold again in case CAMEL replaced GED but is low there too
        low_locations = np.where((ls_emis_final < emis_util.EMIS_LOW_THRESHOLD) &
                                 (ls_emis_final != no_data_value))
        ls_emis_final[low_locations] = no_data_value

        # Add the fill and scan gaps possibly lost during band math steps
        logger.info('Adding fill and data gaps back into the estimated'
                    ' Landsat emissivity results')
        ls_emis_final[thermal_no_data_locations] = no_data_value

        # We want emis_source_data to also apply to emissivity standard
        # deviation, so don't set nodata based on emissivity-specific values.
        # Landsat NDVI captures the background nodata, and it should be close to
        # emissivity standard deviation which uses red band nodata
        if first_band:
            emis_source_data[thermal_no_data_locations] \
                = util.EMIS_SRC_NO_DATA_VALUE

        # Update known water locations to water emissivity value if they are
        # nodata
        (ls_emis_final, emis_source_data) = \
            update_water_locations(src_info, ls_emis_final, emis_source_data,
                                   samps, lines, output_transform,
                                   output_srs.ExportToWkt(), no_data_value,
                                   coefficients.water_emissivity, first_band)

        # Write emissivity data and metadata
        if multiple_bands:
            ls_emis_img_filename = ''.join([product_id, '_emis_', band, '.img'])
        else:
            ls_emis_img_filename = ''.join([product_id, '_emis', '.img'])

        emis_util.write_emissivity_product(samps=samps,
                                           lines=lines,
                                           transform=output_transform,
                                           wkt=output_srs.ExportToWkt(),
                                           no_data_value=no_data_value,
                                           filename=ls_emis_img_filename,
                                           file_data=ls_emis_final,
                                           data_type=gdal.GDT_Float32)

        emis_util.add_emissivity_band_to_xml(espa_metadata=espa_metadata,
                                             filename=ls_emis_img_filename,
                                             camel_filename=camel_filename,
                                             sensor_code=sensor_code,
                                             no_data_value=no_data_value,
                                             band_type='mean',
                                             band_name=band,
                                             multiple_bands=multiple_bands,
                                             resample_method='cubic convolution')

        if first_band:
            # Write emissivity source data and metadata
            emis_source_img_filename = ''.join([product_id, '_emis_source',
                                               '.img'])
            emis_util.write_emissivity_product(samps=samps,
                                               lines=lines,
                                               transform=output_transform,
                                               wkt=output_srs.ExportToWkt(),
                                               no_data_value=
                                                   util.EMIS_SRC_NO_DATA_VALUE,
                                               filename=
                                                   emis_source_img_filename,
                                               file_data=emis_source_data,
                                               data_type=gdal.GDT_Byte)

            emis_util.add_emissivity_band_to_xml(espa_metadata=espa_metadata,
                                                 filename=
                                                     emis_source_img_filename,
                                                 camel_filename=camel_filename,
                                                 sensor_code=sensor_code,
                                                 no_data_value=
                                                     util.EMIS_SRC_NO_DATA_VALUE,
                                                 band_type='emis_source',
                                                 band_name=band,
                                                 multiple_bands=multiple_bands,
                                                 resample_method='nearest neighbor')

        first_band = False

        # Memory cleanup
        del ls_emis_final
        del ls_emis_no_data_locations

    # Cleanup warped ASTER NDVI file
    if not keep_temporary:
        if os.path.exists(aster_ndvi_warped_name):
            os.unlink(aster_ndvi_warped_name)

    # Memory cleanup
    del snow_locations
    del thermal_no_data_locations
    del ls_ndvi_data
    del fv_ls
    del aster_ndvi_data
    del fv_aster


def main():
    """Generate Landsat EMIS and ASTER NDVI from ASTER GED tiles for the
       specified Landsat scene.
    """

    args = emis_util.retrieve_command_line_arguments()

    # Configure logging
    Sys.setup_logging()

    logger.info('*** Begin Generate Estimated Landsat Emissivity ***')

    try:
        # Set GDAL exceptions
        gdal.UseExceptions()

        # Register all the gdal drivers
        gdal.AllRegister()

        # Get the data directory from the environment
        st_data_dir = Sys.get_env_var('ST_DATA_DIR', None)

        args.bands = [band.strip() for band in args.bands]

        # Call the main processing routine
        generate_emissivity_data(xml_filename=args.xml_filename,
                                 server_name=args.aster_ged_server_name,
                                 server_path=args.aster_ged_server_path,
                                 camel_path=args.camel_path,
                                 st_data_dir=st_data_dir,
                                 no_data_value=util.INTERMEDIATE_NO_DATA_VALUE,
                                 keep_temporary=args.keep_temporary,
                                 num_threads=args.num_threads,
                                 use_ias_resampler=args.use_ias_resampler,
                                 band_list=args.bands)

    except Exception:
        logger.exception('Processing failed')
        sys.exit(1)  # EXIT FAILURE

    logger.info('*** Generate Estimated Landsat Emissivity - Complete ***')


if __name__ == '__main__':
    main()