#! /usr/bin/env python ''' FILE: build_st_data.py PURPOSE: Builds the ST product from the intermediate data that was generated. 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 logging import datetime from argparse import ArgumentParser from lxml import objectify as objectify from osgeo import gdal, osr import numpy as np from espa import Metadata, Sys, Geo # Import local modules import st_utilities as util logger = logging.getLogger(__name__) class BuildSTData(object): ''' Description: Defines the processor for generating the Surface Temperature product. ''' def __init__(self, xml_filename, scale, offset): super(BuildSTData, self).__init__() # Keep local copies of this self.xml_filename = xml_filename # Define the scale, offset, and no data values self.scale = scale self.offset = offset self.no_data_value = util.ST_NO_DATA_VALUE self.i_no_data_value = util.INTERMEDIATE_NO_DATA_VALUE self.st_data_dir = '' # Grab the data directory from the environment if 'ST_DATA_DIR' not in os.environ: raise Exception('Environment variable ST_DATA_DIR is' ' not defined') else: self.st_data_dir = os.environ.get('ST_DATA_DIR') # Setup names self.thermal_name = '' self.thermal_input_name = '' self.thermal_input_name_2 = '' self.transmittance_name = '' self.upwelled_name = '' self.downwelled_name = '' self.emissivity_name = '' self.satellite = '' self.product_id = '' def retrieve_metadata_information(self): ''' Description: Loads and reads required information from the metadata XML file. ''' # Read the XML metadata metadata = Metadata(xml_filename=self.xml_filename) # Get satellite and product ID self.satellite = metadata.xml_object.global_metadata.satellite self.product_id = metadata.xml_object.global_metadata.product_id.text self.thermal_name = '' self.thermal_input_name = '' self.thermal_input_name_2 = '' self.transmittance_name = '' self.upwelled_name = '' self.downwelled_name = '' self.emissivity_name = '' # Find the intermediate bands to extract information from for band in metadata.xml_object.bands.band: if (band.get("product") == 'st_intermediate' and band.get("name") == 'st_thermal_radiance'): self.thermal_name = str(band.file_name) if self.satellite == 'LANDSAT_4' or self.satellite == 'LANDSAT_5': if (band.get('product') in ('L1TP', 'L1GT', 'L1GS') and band.get("name") == 'b6'): self.thermal_input_name = str(band.file_name) elif self.satellite == 'LANDSAT_7': if (band.get('product') in ('L1TP', 'L1GT', 'L1GS') and band.get("name") == 'b61'): self.thermal_input_name = str(band.file_name) if (band.get('product') in ('L1TP', 'L1GT', 'L1GS') and band.get("name") == 'b62'): self.thermal_input_name_2 = str(band.file_name) else: # self.satellite == 'LANDSAT_8' or self.satellite == 'LANDSAT_9' if (band.get('product') in ('L1TP', 'L1GT', 'L1GS') and band.get("name") == 'b10'): self.thermal_input_name = str(band.file_name) if (band.get("product") == 'st_intermediate' and band.get("name") == 'st_atmospheric_transmittance'): self.transmittance_name = str(band.file_name) if (band.get("product") == 'st_intermediate' and band.get("name") == 'st_upwelled_radiance'): self.upwelled_name = str(band.file_name) if (band.get("product") == 'st_intermediate' and band.get("name") == 'st_downwelled_radiance'): self.downwelled_name = str(band.file_name) if (band.get("product") == 'st_intermediate' and band.get("name") == 'emis'): self.emissivity_name = str(band.file_name) # Error if we didn't find the required bands in the data if len(self.thermal_name) <= 0: raise Exception('Failed to find the st_thermal_radiance band' ' in the input data') if len(self.thermal_input_name) <= 0: raise Exception('Failed to find the input thermal band (b6, b61, or b10)' ' in the input data') if self.satellite == 'LANDSAT_7': if len(self.thermal_input_name_2) <= 0: raise Exception('Failed to find the second ETM input thermal band b62' ' in the input data') if len(self.transmittance_name) <= 0: raise Exception('Failed to find the st_atmospheric_transmittance' ' in the input data') if len(self.upwelled_name) <= 0: raise Exception('Failed to find the st_upwelled_radiance' ' in the input data') if len(self.downwelled_name) <= 0: raise Exception('Failed to find the st_downwelled_radiance' ' in the input data') if len(self.emissivity_name) <= 0: raise Exception('Failed to find the emis' ' in the input data') del metadata def generate_data(self): ''' Description: Provides the main processing algorithm for building the Surface Temperature product. It produces the final ST product. ''' try: self.retrieve_metadata_information() except Exception: logger.exception('Failed reading input XML metadata file') raise # Set GDAL exceptions gdal.UseExceptions() # Register all the gdal drivers and choose the ENVI for our output gdal.AllRegister() envi_driver = gdal.GetDriverByName('ENVI') # Read the bands into memory # Landsat Radiance at sensor for thermal band logger.info('Loading intermediate thermal band data [%s]', self.thermal_name) dataset = gdal.Open(self.thermal_name) x_dim = dataset.RasterXSize # They are all the same size y_dim = dataset.RasterYSize thermal_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) # Level 1 thermal band logger.info('Loading level 1 thermal band data [%s]', self.thermal_input_name) dataset = gdal.Open(self.thermal_input_name) x_dim = dataset.RasterXSize # They are all the same size y_dim = dataset.RasterYSize thermal_input_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) if self.satellite == 'LANDSAT_7': logger.info('Loading second level 1 thermal band data [%s]', self.thermal_input_name_2) dataset = gdal.Open(self.thermal_input_name_2) x_dim = dataset.RasterXSize # They are all the same size y_dim = dataset.RasterYSize thermal_input_2_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) # Atmospheric transmittance logger.info('Loading intermediate transmittance band data [%s]', self.transmittance_name) dataset = gdal.Open(self.transmittance_name) trans_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) # Atmospheric path radiance - upwelled radiance logger.info('Loading intermediate upwelled band data [%s]', self.upwelled_name) dataset = gdal.Open(self.upwelled_name) upwelled_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) logger.info('Calculating surface radiance') # Surface radiance with np.errstate(invalid='ignore'): surface_radiance = (thermal_data - upwelled_data) / trans_data # Fix the no data locations no_data_locations = np.where(thermal_data == self.i_no_data_value) surface_radiance[no_data_locations] = self.i_no_data_value no_data_locations = np.where(trans_data == self.i_no_data_value) surface_radiance[no_data_locations] = self.i_no_data_value no_data_locations = np.where(upwelled_data == self.i_no_data_value) surface_radiance[no_data_locations] = self.i_no_data_value # Memory cleanup del thermal_data del trans_data del upwelled_data del no_data_locations # Downwelling sky irradiance logger.info('Loading intermediate downwelled band data [%s]', self.downwelled_name) dataset = gdal.Open(self.downwelled_name) downwelled_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) # Landsat emissivity estimated from ASTER GED and/or CAMEL data logger.info('Loading intermediate emissivity band data [%s]', self.emissivity_name) dataset = gdal.Open(self.emissivity_name) emissivity_data = dataset.GetRasterBand(1).ReadAsArray(0, 0, x_dim, y_dim) # Save for the output product ds_srs = osr.SpatialReference() ds_srs.ImportFromWkt(dataset.GetProjection()) ds_transform = dataset.GetGeoTransform() # Memory cleanup del dataset # Estimate Earth-emitted radiance by subtracting off the reflected # downwelling component radiance = (surface_radiance - (1.0 - emissivity_data) * downwelled_data) # Account for surface emissivity to get Plank emitted radiance logger.info('Calculating Plank emitted radiance') with np.errstate(invalid='ignore'): radiance_emitted = radiance / emissivity_data # Fix the no data locations no_data_locations = np.where(surface_radiance == self.i_no_data_value) radiance_emitted[no_data_locations] = self.i_no_data_value no_data_locations = np.where(downwelled_data == self.i_no_data_value) radiance_emitted[no_data_locations] = self.i_no_data_value no_data_locations = np.where(emissivity_data == self.i_no_data_value) radiance_emitted[no_data_locations] = self.i_no_data_value # Memory cleanup del downwelled_data del emissivity_data del surface_radiance del radiance del no_data_locations # Look up radiometric conversion filename for the satellite we're # processing. conversion_lut_name = util.get_radiometric_conversion_filename( self.satellite) # Use Radiometric Conversion LUT to get skin temperature. # For L4-7, the files have the following columns: # # 1. Temperature # 2. Radiance in RTTOV units (wavenumber space) for band 6 # 3. Radiance in microns for band 6 # # For L8-9, the columns are: # # 1. Temperature # 2. Radiance in RTTOV units (wavenumber space) for band 10 # 3. Radiance in RTTOV units (wavenumber space) for band 11 # 4. Radiance in microns for band 10 # 5. Radiance in microns for band 11 # # In this phase we only use the temperature and radiance in microns for # band 10. if self.satellite == 'LANDSAT_4' or self.satellite == 'LANDSAT_5' \ or self.satellite == 'LANDSAT_7': radiance_col = 2 else: # self.satellite == 'LANDSAT_8' or self.satellite == 'LANDSAT_9' radiance_col = 3 conversion_data = np.loadtxt(os.path.join(self.st_data_dir, conversion_lut_name), dtype=float, delimiter=' ', usecols=(0,radiance_col)) temp_lut = conversion_data[:, 0] radiance_lut = conversion_data[:, 1] # B10 radiance in microns logger.info('Generating ST results') st_data = np.interp(radiance_emitted, radiance_lut, temp_lut) # Scale the result st_data = (st_data - self.offset) / self.scale # Mark the locations that are saturated in the level 1 thermal input if self.satellite == 'LANDSAT_4' or self.satellite == 'LANDSAT_5' \ or self.satellite == 'LANDSAT_7': saturation = util.TM_ETM_HIGH_SATURATION else: # self.satellite == 'LANDSAT_8' or self.satellite == 'LANDSAT_9' saturation = util.OLI_HIGH_SATURATION if self.satellite == 'LANDSAT_7': saturation_locations = np.where((thermal_input_data == saturation) & (thermal_input_2_data == saturation)) else: saturation_locations = np.where(thermal_input_data == saturation) # Remove saturated locations st_data[saturation_locations] = self.no_data_value scaled_min = (util.ST_MIN - self.offset) / self.scale scaled_max = (util.ST_MAX - self.offset) / self.scale # Add the fill and scan gaps back into the results, since they may # have been lost logger.info('Adding fill and data gaps back into the Surface' ' Temperature results') # Fix the no data locations no_data_locations = np.where(radiance_emitted == self.i_no_data_value) st_data[no_data_locations] = self.no_data_value # Memory cleanup del radiance_emitted del no_data_locations st_img_filename = ''.join([self.product_id, '_st', '.img']) st_hdr_filename = ''.join([self.product_id, '_st', '.hdr']) st_aux_filename = ''.join([st_img_filename, '.aux', '.xml']) logger.info('Creating %s', st_img_filename) Geo.generate_raster_file(envi_driver, st_img_filename, st_data, x_dim, y_dim, ds_transform, ds_srs.ExportToWkt(), self.no_data_value, gdal.GDT_UInt16) logger.info('Updating %s', st_hdr_filename) util.ST_Geo.update_envi_header(st_hdr_filename, self.no_data_value) # Memory cleanup del ds_srs del ds_transform # Remove the *.aux.xml file generated by GDAL if os.path.exists(st_aux_filename): os.unlink(st_aux_filename) logger.info('Adding %s to %s', st_img_filename, self.xml_filename) # Add the estimated Surface Temperature product to the metadata metadata = Metadata(xml_filename=self.xml_filename) # Create an element maker and band element em = objectify.ElementMaker(annotate=False, namespace=None, nsmap=None) st_band = em.band() # Find the level 1 thermal band to use for the specific band details base_band = util.get_thermal_band_metadata(metadata) # Set attributes for the band element st_band.set('product', 'st') st_band.set('source', 'level_1_thermal') st_band.set('name', 'surface_temperature') st_band.set('category', 'image') st_band.set('data_type', 'UINT16') st_band.set('scale_factor', "%g" % self.scale) st_band.set('add_offset', "%g" % self.offset) st_band.set('nlines', str(int(base_band.get("nlines")))) st_band.set('nsamps', str(int(base_band.get("nsamps")))) st_band.set('fill_value', str(self.no_data_value)) # Add elements to the band object st_band.short_name = em.element('{0}ST'.format(self.product_id[0:4])) st_band.long_name = em.element('Surface Temperature') st_band.file_name = em.element(str(st_img_filename)) # Create a pixel size element and add attributes to it st_band.pixel_size = em.element() st_band.pixel_size.set('x', str(base_band.pixel_size.get('x'))) st_band.pixel_size.set('y', str(base_band.pixel_size.get('x'))) st_band.pixel_size.set('units', str(base_band.pixel_size.get('units'))) st_band.resample_method = em.element('none') st_band.data_units = em.element('temperature (kelvin)') # Create a valid range element and add attributes to it st_band.valid_range = em.element() st_band.valid_range.set('min', "%g" % scaled_min) st_band.valid_range.set('max', "%g" % scaled_max) # Populate auxiliary_source information st_band.auxiliary_source = em.element() st_band.auxiliary_source.set('reanalysis', util.REANALYSIS.get_reanalysis()) # Per discussion, just use GLS DEM now, not what's in elevation_source # Entries like RAMP should be considered part of Collection 2 GLS-DEM st_band.auxiliary_source.set('dem', 'GLS DEM') st_band.app_version = em.element(str(util.Version.app_version())) # Set the date, but first clean the microseconds off of it production_date = ('{0}Z'.format(datetime.datetime.utcnow() .strftime('%Y-%m-%dT%H:%M:%S'))) st_band.production_date = em.element(str(production_date)) # Add the new band to the XML, validate it, and write it metadata.xml_object.bands.append(st_band) metadata.validate() metadata.write(xml_filename=self.xml_filename) # Memory cleanup del metadata del st_band def main(): ''' Description: Create Surface Temperature band using the thermal, emissivity, and atmospheric parameter bands. ''' # Build the command line argument parser parser = ArgumentParser(description='Reads intermediate data generated' ' previously and combines them into' ' the Surface Temperature product') # ---- Add parameters ---- # Required parameters parser.add_argument('--xml', action='store', dest='xml_filename', required=False, default=None, help='The XML metadata file to use') parser.add_argument('--scale', action='store', dest='scale', required=False, default=util.DEFAULT_SCALE, help='The scale factor for ST products') parser.add_argument('--offset', action='store', dest='offset', required=False, default=util.DEFAULT_OFFSET, help='The offset value for ST products') # Optional parameters parser.add_argument('--version', action='version', version=util.Version.version_text()) # Parse the command line arguments args = parser.parse_args() # Command line arguments are required so print the help if none were # provided if len(sys.argv) == 1: parser.print_help() sys.exit(1) # EXIT FAILURE # Verify that the --xml parameter was specified if args.xml_filename is None: raise Exception('--xml must be specified on the command line') # Configure logging Sys.setup_logging() try: build_st_data = BuildSTData(args.xml_filename, args.scale, args.offset) # Call the main processing routine build_st_data.generate_data() except Exception: logger.exception('Processing failed') sys.exit(1) # EXIT FAILURE sys.exit(0) # EXIT SUCCESS if __name__ == '__main__': main()