#! /usr/bin/env python ''' FILE: emissivity_utilities.py PURPOSE: Provide a library of routines to be used by ST emissivity applications. 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 urllib.parse import logging import datetime from argparse import ArgumentParser from collections import namedtuple import requests import lxml from lxml import objectify as objectify from osgeo import gdal from espa import Meta, Sys, Geo, Web from st_exceptions import MissingBandError # Import local modules import st_utilities as util logger = logging.getLogger(__name__) # Formats to use for ASTER GED tile names # Filename format is modifiable by command-line argument __aster_ged_filename_format = '' # Format latitude as optional negative sign plus 2 digits (0-padded) ASTER_GED_LAT_FORMAT = '{0: 03d}' # Format longitude as optional negative sign plus 3 digits (0-padded) ASTER_GED_LON_FORMAT = '{0: 04d}' EMIS_LOW_THRESHOLD = 0.6 def data_resolution_and_size(name, x_min, x_max, y_min, y_max): """Calculates dataset resolution and retrieves size information Args: name : Full path dataset name x_min : Minimum longitude value x_max : Maximum longitude value y_min : Minimum latitude value y_max : Maximum latitude value Returns: : Longitude resolution : Latitude resolution : Samples in the data : Lines in the data """ dataset = gdal.Open(name) if dataset is None: raise RuntimeError('GDAL failed to open {0}'.format(name)) return ((x_max - x_min) / float(dataset.RasterXSize), (y_max - y_min) / float(dataset.RasterYSize), dataset.RasterXSize, dataset.RasterYSize) XYInfo = namedtuple('XYInfo', ('x', 'y')) BandInfo = namedtuple('BandInfo', ('name', 'scale_factor', 'add_offset', 'pixel_size', 'fill_value')) BandTypeInfo = namedtuple('BandTypeInfo', ('red', 'nir', 'thermal')) ExtentInfo = namedtuple('ExtentInfo', ('min', 'max')) BoundInfo = namedtuple('BoundInfo', ('north', 'south', 'east', 'west')) SourceInfo = namedtuple('SourceInfo', ('bound', 'extent', 'wkt', 'pixel_qa', 'band_type')) TileInfo = namedtuple('TileInfo', ('h5_file_path', 'downloaded')) def get_band_info(band): """Returns a populated BandInfo Args: band : Current band being processed Returns: : Populated with band information """ if band.get('name') in ['b10', 'b6', 'b61']: return BandInfo(name=str(band.file_name), scale_factor=float(band.radiance.get('gain')), add_offset=float(band.radiance.get('bias')), fill_value=float(band.get('fill_value')), pixel_size=XYInfo(x=float(band.pixel_size.get('x')), y=float(band.pixel_size.get('y')))) else: return BandInfo(name=str(band.file_name), scale_factor=float(band.get('scale_factor')), add_offset=float(band.get('add_offset')), fill_value=float(band.get('fill_value')), pixel_size=XYInfo(x=float(band.pixel_size.get('x')), y=float(band.pixel_size.get('y')))) def retrieve_metadata_information(espa_metadata, tirs_only): """Reads required information from the metadata XML file Args: espa_metadata : XML metadata tirs_only : Is this a TIRS-only scene? Returns: : Populated with source information """ bi_red = None bi_nir = None bi_thermal = None pixel_qa_name = None satellite = espa_metadata.xml_object.global_metadata.satellite # Find the bands to extract information from. Use SR if available for band in espa_metadata.xml_object.bands.band: if satellite == 'LANDSAT_8' or satellite == 'LANDSAT_9': if (band.get('product') == 'sr_refl' and band.get('name') == 'sr_band4'): bi_red = get_band_info(band) if (band.get('product') == 'sr_refl' and band.get('name') == 'sr_band5'): bi_nir = get_band_info(band) if bi_red is None: if (band.get('product') == 'toa_refl' and band.get('name') == 'toa_band4'): bi_red = get_band_info(band) if bi_nir is None: if (band.get('product') == 'toa_refl' and band.get('name') == 'toa_band5'): bi_nir = get_band_info(band) if (band.get('product').startswith('L1') and band.get('name') == 'b10'): bi_thermal = get_band_info(band) else: if (band.get('product') == 'sr_refl' and band.get('name') == 'sr_band3'): bi_red = get_band_info(band) if (band.get('product') == 'sr_refl' and band.get('name') == 'sr_band4'): bi_nir = get_band_info(band) if bi_red is None: if (band.get('product') == 'toa_refl' and band.get('name') == 'toa_band3'): bi_red = get_band_info(band) if bi_nir is None: if (band.get('product') == 'toa_refl' and band.get('name') == 'toa_band4'): bi_nir = get_band_info(band) if satellite == 'LANDSAT_4' or satellite == 'LANDSAT_5': if (band.get('product').startswith('L1') and band.get('name') == 'b6'): bi_thermal = get_band_info(band) elif satellite == 'LANDSAT_7': if (band.get('product').startswith('L1') and band.get('name') == 'b61'): bi_thermal = get_band_info(band) # Get metadata information for level 1 BQA pixel band if (band.get('product').startswith('L1') and band.get('name') == 'qa_pixel'): pixel_qa_name = str(band.file_name) # Error if we didn't find the required SR or TOA bands in the data if not tirs_only: if bi_red is None: raise MissingBandError('Failed to find the SR or TOA RED band' ' in the input data') if bi_nir is None: raise MissingBandError('Failed to find the SR or TOA NIR band' ' in the input data') if bi_thermal is None: raise MissingBandError('Failed to find the Level 1 thermal band' ' in the input data') if pixel_qa_name is None: raise MissingBandError('Failed to find the PIXEL QA band' ' in the input data') # Get the output wkt string wkt = Geo.get_wkt_projection_string(bi_thermal.name) return SourceInfo(bound=Meta.bound_info(espa_metadata), extent=Meta.extent_info(espa_metadata, bi_thermal), wkt=wkt, pixel_qa=pixel_qa_name, band_type=BandTypeInfo(red=bi_red, nir=bi_nir, thermal=bi_thermal)) def get_aster_ged_tiles_for_src(st_data_dir, src_info, antimeridian_crossing): """Gets the names of ASTER GED tiles for the region of the source image Args: st_data_dir : Location of the ST data file src_info : Information about the source data antimeridian_crossing : Flag for scene crossing 180 meridian Returns: )>: List of ASTER GED files to process )>: List of lat/lon values including missing ASTER GED tiles """ # Read the ASTER GED tile list ged_tile_file = 'aster_ged_tile_list.txt' with open(os.path.join(st_data_dir, ged_tile_file)) as ged_file: tiles = [os.path.splitext(line.rstrip('\n'))[0] for line in ged_file] filename_format = get_aster_ged_filename_format() if antimeridian_crossing: lon_range = list(range(int(src_info.bound.west), 180)) \ + list(range(int(src_info.bound.east), -181, -1)) else: lon_range = list( range(int(src_info.bound.west), int(src_info.bound.east)+1)) ged_file_list = [] camel_tile_list = [] for (lat, lon) in [(lat, lon) for lat in range(int(src_info.bound.south), int(src_info.bound.north)+1) for lon in lon_range]: # Build the base filename using the correct format filename = filename_format.format( ASTER_GED_LAT_FORMAT.format(lat).strip(), ASTER_GED_LON_FORMAT.format(lon).strip()) # Skip the tile if it isn't in the ASTER GED tile list # (ignore filename extension) tile_name = '.'.join(filename.split('.')[:5]) if tile_name in tiles: ged_file_list.append(filename) camel_tile_list.append( ASTER_GED_LAT_FORMAT.format(lat).strip() + "." + ASTER_GED_LON_FORMAT.format(lon).strip()) return ged_file_list, camel_tile_list def locate_aster_ged_tile(url, filename): """Locate the specified tile, either on disk or download if needed Args: url : URL to retrieve the file from filename : Tile filename Returns: : tile filename and downloaded flag """ h5_file_path = None downloaded = False # If the file exists in the current working directory, that means it was # previously downloaded to here if os.path.exists(filename): h5_file_path = filename downloaded = True else: # See if the file is accessible locally, but in another directory # If URL is just a path local_h5_file_path = os.path.join(url, filename) if os.path.exists(local_h5_file_path): h5_file_path = local_h5_file_path else: # Try parsing the URL, if url includes file://hostname/path url_parts = urllib.parse.urlparse(local_h5_file_path) local_h5_file_path = os.path.abspath(os.path.join(url_parts.netloc, url_parts.path)) if os.path.exists(local_h5_file_path): h5_file_path = local_h5_file_path else: download_aster_ged_tile(url=url, h5_file_path=filename) h5_file_path = filename downloaded = True return TileInfo(h5_file_path, downloaded) def download_aster_ged_tile(url, h5_file_path): """Retrieves the specified tile from the host Args: url : URL to retrieve the file from h5_file_path : Full path on the remote system Raises: Exception: If issue transfering data """ # Build the complete URL and download the tile url_path = ''.join([url, h5_file_path]) status_code = Web.http_transfer_file(url_path, h5_file_path) # If a tile is requested, it is expected to be there if status_code != requests.codes['ok']: raise Exception('HTTP - Transfer Failed') def warp_raster(odl_file): """Executes makegeomgrid and geomresample using the information from the supplied odl file to warp to a specific location and extent Args: odl_file : Name of the source ODL file """ logger = logging.getLogger(__name__) # Need to first call makegeomgrid output = '' makegeomgridcmd = "makegeomgrid " + odl_file try: logger.info('Executing [{0}]'.format(makegeomgridcmd)) output = Sys.execute_cmd(makegeomgridcmd) except Exception: logger.error('Failed during makegeomgrid') raise finally: if len(output) > 0: logger.info(output) # Now call geomresample output = '' geomresamplecmd = "geomresample " + odl_file try: logger.info('Executing [{0}]'.format(geomresamplecmd)) output = Sys.execute_cmd(geomresamplecmd) except Exception: logger.error('Failed during geomresample') raise finally: if len(output) > 0: logger.info(output) def shift_longitude(tile_name, shifted_tile_name, offset): """Shift the longitude of the tile data and put the results in the requested output file Args: tile_name : Filename of tile to shift shifted_tile_name : Filename of output tile to place results offset : Number of degrees to add to the longitude """ # Set up the base command cmd = ['gdal_translate', '-a_ullr'] # Get the current locations tile_src = gdal.Open(tile_name) ulx, xres, xskew, uly, yskew, yres = tile_src.GetGeoTransform() # Compute the adjusted longitude locations lrx = ulx + (tile_src.RasterXSize * xres) lry = uly + (tile_src.RasterYSize * yres) new_ulx = ulx + offset new_lrx = lrx + offset # Close the dataset tile_src = None # Add updated coordinates to the command cmd.extend([str(new_ulx), str(uly), str(new_lrx), str(lry)]) # Add source and destination files to the command cmd.extend([tile_name, shifted_tile_name]) # Convert to a string for the execution cmd = ' '.join(cmd) output = '' try: logger.info('Executing [{0}]'.format(cmd)) output = Sys.execute_cmd(cmd) finally: if len(output) > 0: logger.info(output) # Shift tile longitudes def shift_tiles(tiles): """Shift the longitude of the tiles that need it to ensure they are in the 0..360 range. This is intended to be used in antimeridian crossing cases. Making the longitude range contiguous enables mosaicking using GDAL tools. Args: tiles )>: List of tiles to shift """ for tile in tiles: longitude = int(tile.split(".")[3]) # Only shift the longitude of the tiles with negative longitude if longitude < 0: # Name the shifted output file shifted_tile = tile + '_shifted' # Shift the longitude values shift_longitude(tile, shifted_tile, 360) # Move destination file back to source file output = '' try: cmd = 'mv {0} {1}'.format(shifted_tile, tile) logger.info('Executing [{0}]'.format(cmd)) output = Sys.execute_cmd(cmd) finally: if len(output) > 0: logger.info(output) def write_emissivity_product(samps, lines, transform, wkt, no_data_value, filename, file_data, data_type): """Creates the emissivity band file Args: samps : Samples in the data lines : 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 : Well-Known-Text describing the projection no_data_value : Value to use for fill filename : Full path for the output file to create file_data : Binary image data to be output data_type: GDAL data type """ logger.info('Creating {0}'.format(filename)) Geo.generate_raster_file(gdal.GetDriverByName('ENVI'), filename, file_data, samps, lines, transform, wkt, no_data_value, data_type) hdr_filename = Sys.replace_extension(filename, '.hdr') logger.info('Updating {0}'.format(hdr_filename)) util.ST_Geo.update_envi_header(hdr_filename, no_data_value) # Remove the *.aux.xml file generated by GDAL aux_filename = Sys.replace_extension(filename, '.img.aux.xml') if os.path.exists(aux_filename): os.unlink(aux_filename) def add_emissivity_band_to_xml(espa_metadata, filename, camel_filename, sensor_code, no_data_value, band_type, band_name, multiple_bands, resample_method): """Adds the emissivity band to the Metadata XML file Args: espa_metadata : XML metadata information filename : Full path for the output file to create camel_filename : Full path for the source file for CAMEL data sensor_code : Name prefix for the sensor no_data_value : Value to use for fill band_type : Emissivity mean, standard deviation, or source band_name : band to process multiple_bands : Is more than one band being processed? resample_method : gdalwarp resample method string """ logger.info('Adding {0} to Metadata XML'.format(filename)) # Create an element maker maker = objectify.ElementMaker(annotate=False, namespace=None, nsmap=None) # Find the SR Band 3 to use for the specific band details sr_product = 'sr_refl' # This really only applies to ASTER GED pixels toa_product = 'toa_refl' # This really only applies to ASTER GED pixels camel_product = 'camel' source_product = None base_band = util.get_thermal_band_metadata(espa_metadata) for band in espa_metadata.xml_object.bands.band: if (band.get('product') == sr_product and band.get('name') == 'sr_band3'): source_product = sr_product break # If SR Band 3 isn't available, find TOA Band 3 if source_product is None: for band in espa_metadata.xml_object.bands.band: if (band.get('product') == toa_product and band.get('name') == 'toa_band3'): source_product = toa_product break # If neither SR nor TOA are available, it's TIRS-only so there is no SR # or TOA NDVI used and the source is CAMEL if source_product is None: source_product = camel_product emis_band = maker.band() emis_band.set('product', 'st_intermediate') emis_band.set('source', source_product) if band_type == 'mean': emis_band.set('name', 'emis') elif band_type == 'stdev': emis_band.set('name', 'emis_stdev') else: # band_type == 'emis_source' emis_band.set('name', 'emis_source') emis_band.set('category', 'image') if band_type == 'mean' or band_type == 'stdev': emis_band.set('data_type', 'FLOAT32') else: # band_type == 'emis_source' emis_band.set('data_type', 'UINT8') emis_band.set('nlines', base_band.attrib['nlines']) emis_band.set('nsamps', base_band.attrib['nsamps']) if band_type == 'mean': emis_band.short_name = maker.element('{0}EMIS'.format(sensor_code)) if multiple_bands: emis_band.set('name', 'emis_' + band_name) emis_band.long_name = \ maker.element('Landsat emissivity estimated from ASTER GED or' ' CAMEL data for {0}'.format(band_name)) else: emis_band.set('name', 'emis') emis_band.long_name = maker.element('Landsat emissivity estimated' ' from ASTER GED or CAMEL data') elif band_type == 'stdev': emis_band.short_name \ = maker.element('{0}EMIS_STDEV'.format(sensor_code)) if multiple_bands: emis_band.set('name', 'emis_stdev_' + band_name) emis_band.long_name = maker.element('Landsat emissivity standard' ' deviation estimated from' ' ASTER GED or CAMEL data for' ' {0}'.format(band_name)) else: emis_band.set('name', 'emis_stdev') emis_band.long_name = maker.element('Landsat emissivity standard' ' deviation estimated from' ' ASTER GED or CAMEL data') else: # band_type == 'emis_source' emis_band.short_name \ = maker.element('{0}EMIS_STDEV'.format(sensor_code)) emis_band.set('name', 'emis_source') emis_band.long_name = maker.element('Landsat emissivity standard' ' deviation estimated from' ' ASTER GED or CAMEL data') emis_band.short_name \ = maker.element('{0}EMIS_SOURCE'.format(sensor_code)) emis_band.long_name = maker.element('Landsat emissivity source used') emis_band.set('fill_value', str(no_data_value)) emis_band.file_name = maker.element(filename) emis_band.pixel_size = base_band.pixel_size emis_band.resample_method = maker.element(resample_method) if band_type == 'mean' or band_type == 'stdev': emis_band.data_units = maker.element('Emissivity Coefficient') else: # band_type == 'emis_source' emis_band.data_units = maker.element('Emissivity source classification') emis_band.valid_range = maker.element() emis_band.valid_range.set('min', '0.0') if band_type == 'mean' or band_type == 'stdev': emis_band.valid_range.set('max', '1.0') else: # band_type == 'emis_source' emis_band.valid_range.set('max', '3.0') # Populate auxiliary source field with the emissivity source information if band_type == 'mean' or band_type == 'stdev': base = os.path.basename(camel_filename) emis_band.auxiliary_source = maker.element() if source_product == camel_product: emis_band.auxiliary_source.set('emissivity', 'CAMEL (' + base + ')') else: emis_band.auxiliary_source.set('emissivity', 'ASTER GED AG100 v003 and CAMEL (' + base + ')') # The emis_source band has a class_values XML structure if band_type == 'emis_source': class_values_str = ('ASTER GEDCAMELWaterSnow') class_values = lxml.objectify.fromstring(class_values_str) emis_band.append(class_values) if multiple_bands: emis_band.app_version = maker.element( util.Version.split_window_app_version()) else: emis_band.app_version = maker.element(util.Version.app_version()) # Get the production date and time in string format # Strip the microseconds and add a Z date_now = ('{0}Z'.format(datetime.datetime.utcnow() .strftime('%Y-%m-%dT%H:%M:%S'))) emis_band.production_date = maker.element(date_now) # Append the band to the XML espa_metadata.xml_object.bands.append(emis_band) # Validate the XML espa_metadata.validate() # Write it to the XML file espa_metadata.write() def retrieve_command_line_arguments(): """Build the command line argument parser with some extra validation Returns: : The command line arguments """ description = ('Estimates Landsat Emissivity from ASTER GED data') parser = ArgumentParser(description=description) parser.add_argument('--version', action='version', version=util.Version.version_text()) parser.add_argument('--xml', action='store', dest='xml_filename', required=False, default=None, help='The XML metadata file to use') parser.add_argument('--aster-ged-server-name', action='store', dest='aster_ged_server_name', required=False, default=None, help='Name of the ASTER GED server') parser.add_argument('--aster-ged-server-path', action='store', dest='aster_ged_server_path', required=False, default=None, help='Path on the ASTER GED server') parser.add_argument('--camel-path', action='store', dest='camel_path', required=False, default=None, help='Path to the CAMEL archive') parser.add_argument('--aster-ged-filename-format', action='store', dest='aster_ged_filename_format', required=False, default='AG100.v003.{0}.{1}.0001.subset.h5', help='ASTER GED filename format, default ' + 'AG100.v003.{0}.{1}.0001.subset.h5') parser.add_argument('--keep-temporary-data', action='store_true', dest='keep_temporary', required=False, default=False, help='Keep any temporary products generated') parser.add_argument('--num_threads', action='store', dest='num_threads', required=False, default=2, help='Number of threads to use for multithreaded ' + 'components') # Use band10, band11, or band10 and band11 for L8. If no bands are # specified, L8 and L9 use band10. This is ignored for L4, L5, and L7, # which automatically use band6. parser.add_argument('--band', action='append', dest='bands', required=False, help='For L8 and L9, names of bands [band10|band11]') parser.add_argument('--use_ias_resampler', action='store_true', dest='use_ias_resampler', required=False, default=False, help='Do not use ESPA gdalwarp (use IAS ' + 'makegeomgrid and geomresample)') 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') if args.aster_ged_server_name is None: raise Exception('--aster-ged-server-name must be specified on the' ' command line') if args.aster_ged_server_name == '': raise Exception('The --aster-ged-server-name provided was empty') if args.aster_ged_server_path is None: raise Exception('--aster-ged-server-path must be specified on the' ' command line') if args.aster_ged_server_path == '': raise Exception('The --aster-ged-server-path provided was empty') global __aster_ged_filename_format __aster_ged_filename_format = args.aster_ged_filename_format # Set default band if not specified if args.bands is None: args.bands = ["band10"] return args def get_aster_ged_filename_format(): return __aster_ged_filename_format def create_omf(source_projection_file, target_projection_file, num_threads): """ Creates OMF for the makegeomgrid and geomgrid applications Args: source_projection_file : File name to use for the input dem target_projection_file : File with the desired projection information """ logger = logging.getLogger(__name__) if source_projection_file.find("aster") != -1: # Creating the aster omf, set up the appropriate variables file_type = "ASTER" else: file_type = "EMIS" omf_name = file_type + "WARP.omf" if num_threads is None: num_threads = 1 logger.info('Creating {0}'.format(omf_name)) with open(omf_name, 'w') as omf: # Write omf parameters omf.write("OBJECT = {0}_WARP_OMF\n".format(file_type)) omf.write(' DEM_FILENAME = \"{0}\"\n'.format(source_projection_file)) omf.write(' GRID_FILENAME_PASS_1 = \"{0}\"\n'. format(target_projection_file)) omf.write(' MAX_PROCESSORS_TO_USE = {0}\n'.format(num_threads)) # Set the flag to buffer the entire output image and then write # at the end so that the TIFs get ordered consistently omf.write(' BUFFER_OUTPUT_IMAGE = 1\n') omf.write("END_OBJECT = {0}_WARP_OMF\n".format(file_type)) omf.write('END\n') def create_odl(output_name, no_data_value): """ Creates ODL for the makegeomgrid and geomgrid applications Args: output_name : Output file name for geomgrid no_data_value : Value to use for fill """ logger = logging.getLogger(__name__) if output_name.find("aster") != -1: file_type = "ASTER" geom_grid_name = "aster_ndvi_mosaic.grd" else: file_type = "EMIS" geom_grid_name = "landsat_emis_mosaic.grd" # makegeomgrid and geomresample use the work order to determine the # omf name, so the work order has to correspond to the OMF name work_order = file_type + "WARP" odl_name = file_type + "WARP.odl" logger.info('Creating {0}'.format(odl_name)) with open(odl_name, 'w') as odl: # Write ODL parameters odl.write("OBJECT = {0}_WARP_ODL\n".format(file_type)) odl.write(' WORK_ORDER_ID = {0}\n'.format(work_order)) odl.write(' PROCESSING_PASS = 1\n') odl.write(' BAND_LIST = 1\n') odl.write(' CELL_LINES = 25\n') odl.write(' CELL_SAMPLES = 25\n') odl.write(' GEOM_GRID_FILENAME = \"{0}\"\n'. format(geom_grid_name)) odl.write(' SOURCE_BAND_NUMBER_LIST = 1\n') odl.write(' TARGET_BAND_NUMBER_LIST = 0\n') odl.write(' OUTPUT_IMAGE_FILENAME = \"{0}\"\n'.format(output_name)) odl.write(' ODTYPE = \"R*4\"\n') odl.write(' RESAMPLE = NN\n') odl.write(' BACKGRND = {0}\n'.format(no_data_value)) odl.write("END_OBJECT = {0}_WARP_ODL\n".format(file_type)) odl.write('END\n') return odl_name