#! /usr/bin/env python

'''
    File: st_build_elevation.py

    Purpose: Build elevation band that covers the reanalysis points to be run.
             The grid points must be defined in the binary grid point file.

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

    License: NASA Open Source Agreement 1.3

'''

import os
import sys
import logging
from argparse import ArgumentParser

from st_exceptions import MissingBandError
from st_grid_points import read_grid_points
from espa import Metadata, Sys
import st_utilities as util


def retrieve_command_line_arguments():
    """Read arguments from the command line

    Returns:
        args <arguments>: The arguments read from the command line
    """

    parser = ArgumentParser(description='Convert intermediate bands to int16')

    parser.add_argument('--xml',
                        action='store', dest='xml_filename',
                        required=True, default=None,
                        help='The XML metadata file to use')

    args = parser.parse_args()

    # Verify that the --xml parameter was specified
    if args.xml_filename is None:
        raise Exception('--xml must be specified on the command line')

    return args


def retrieve_metadata_information(espa_metadata):
    """Reads required information from a reference band from the metadata XML
       file

    Args:
        espa_metadata <espa.Metadata>: XML metadata

    Returns:
        pixel_size_x, pixel_size_y
    """

    pixel_size_x = 0
    pixel_size_y = 0

    # Use the input thermal band as the source for the pixel size information
    band = util.get_thermal_band_metadata(espa_metadata)
    pixel_size_x = float(band.pixel_size.get('x'))
    pixel_size_y = float(band.pixel_size.get('y'))

    # Error if we didn't find the required intermediate band in the data
    if pixel_size_x == 0 or pixel_size_y == 0:
        raise MissingBandError('Failed to find the reference band in the'
                               ' input data')

    return pixel_size_x, pixel_size_y


def next_map_pixel(map_value, pixel_size):
    """Adjust map value in the expanding direction to be divisible by pixel 
       size

    Args:
        map_value <float>: Map coordinate value 
        pixel_size <float>: Size of a pixel in map coordinate value's direction 
    """

    quotient = int(map_value / pixel_size)

    if map_value > 0:
        return pixel_size * (quotient + 1)
    else:
        return pixel_size * quotient


def build_elevation(xml_filename):
    """Build an elevation file that includes data for all reanalysis points

    Args:
        xml_filename <str>: XML metadata filename
    """

    logger = logging.getLogger(__name__)

    # Initialize boundary values.
    min_lon = sys.float_info.max
    min_lat = sys.float_info.max
    max_lon = -sys.float_info.max
    max_lat = -sys.float_info.max
    min_x = sys.float_info.max
    min_y = sys.float_info.max
    max_x = -sys.float_info.max
    max_y = -sys.float_info.max

    # XML Metadata
    espa_metadata = Metadata(xml_filename)
    espa_metadata.parse()

    # Look up the pixel size from the metadata using a reference band.
    (pixel_size_x, pixel_size_y) = retrieve_metadata_information(espa_metadata)

    # Lookup up the longitude boundaries from the metadata.
    east_lon = float(espa_metadata.xml_object.
                     global_metadata.bounding_coordinates.east)
    west_lon = float(espa_metadata.xml_object.
                     global_metadata.bounding_coordinates.west)

    # Check if it's an antimeridian-crossing scene.
    if east_lon > west_lon:
        cross_antimeridian = False
    else:
        cross_antimeridian = True

    # Load the grid information.
    (grid_points, dummy1, dummy2) = read_grid_points()

    # There are cases where it's an antimeridian crossing scene but there are
    # no reanalysis points on one of the sides.  Treat those like non-
    # crossers.
    no_negative_longitude = True
    no_positive_longitude = True
    for point in grid_points:

        # Ignore reanalysis points that will not be processed.
        if not point.run_model:
            continue
        if point.lon < 0:
            no_negative_longitude = False
        else:
            no_positive_longitude = False
    if no_negative_longitude or no_positive_longitude:
        cross_antimeridian = False

    for point in grid_points:

        # Ignore reanalysis points that will not be processed.
        if not point.run_model:
            continue

        # Find minimum and maximum lat and lon values from the reanalysis points
        # that will be used.
        if point.lat < min_lat:
            min_lat = point.lat
        if point.lat > max_lat:
            max_lat = point.lat
        if cross_antimeridian:
            # Handle east: find the largest negative
            if point.lon < 0:
                if point.lon > max_lon:
                    max_lon = point.lon
            # Handle west: find the smallest positive
            else:
                if point.lon < min_lon:
                    min_lon = point.lon
        else:
            if point.lon < min_lon:
                min_lon = point.lon
            if point.lon > max_lon:
                max_lon = point.lon

        # Find minimum and maximum x and y values from the reanalysis points
        # that will be used.
        if point.map_x < min_x:
            min_x = point.map_x
        if point.map_y < min_y:
            min_y = point.map_y
        if point.map_x > max_x:
            max_x = point.map_x
        if point.map_y > max_y:
            max_y = point.map_y

    # This gives us the exact boundaries for the points, but we want to
    # adjust the values to be divisible by the pixel size.
    min_x = next_map_pixel(min_x, pixel_size_x)
    min_y = next_map_pixel(min_y, pixel_size_y)
    max_x = next_map_pixel(max_x, pixel_size_x)
    max_y = next_map_pixel(max_y, pixel_size_y)

    # Adjust the center of pixel coordinates to get the actual extent (0.5).
    # Also adjust the map coordinates by 1 pixel so small differences don't
    # cause us to miss the reanalysis points on the edge.
    min_x -= pixel_size_x * 1.5
    min_y -= pixel_size_y * 1.5
    max_x += pixel_size_x * 1.5
    max_y += pixel_size_y * 1.5

    # Build ST elevation filename.  This should be a separate file from the
    # regular elevation file in case both are used in the processing flow.
    elevation_filename = ''.join([xml_filename.split('.xml')[0],
                                  '_st_elevation', '.img'])

    # Build the elevation file.  Use the minimum and maximum boundary values
    # defined by the reanalysis point list to define the extent of the file.
    output = ''
    try:
        cmd = ['build_elevation_band.py',
               '--xml', xml_filename,
               '--elevation', elevation_filename,
               '--extent-minx', str(min_x),
               '--extent-maxx', str(max_x),
               '--extent-miny', str(min_y),
               '--extent-maxy', str(max_y),
               '--nbound-lat', str(max_lat),
               '--sbound-lat', str(min_lat),
               '--wbound-lon', str(min_lon),
               '--ebound-lon', str(max_lon)]

        output = Sys.execute_cmd(' '.join(cmd))
    finally:
        if len(output) > 0:
            logger = logging.getLogger(__name__)
            logger.info(output)


def main():
    """Main processing for building elevation band
    """

    # Command Line Arguments
    args = retrieve_command_line_arguments()

    # Check logging level
    logging_level = logging.INFO

    # Setup the default logger format and level.  Log to STDOUT.
    logging.basicConfig(format=('%(asctime)s.%(msecs)03d %(process)d'
                                ' %(levelname)-8s'
                                ' %(filename)s:%(lineno)d:'
                                '%(funcName)s -- %(message)s'),
                        datefmt='%Y-%m-%d %H:%M:%S',
                        level=logging_level,
                        stream=sys.stdout)
    logger = logging.getLogger(__name__)

    logger.info('*** Begin ST Build Elevation Band ***')

    try:
        # Call the main processing routine
        build_elevation(xml_filename=args.xml_filename)

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

    logger.info('*** ST Build Elevation Band - Complete ***')

if __name__ == '__main__':
    main()