Source code for dnppy.landsat.toa_radiance

#standard imports
from grab_meta import grab_meta
from dnppy import core
import arcpy
import os
if arcpy.CheckExtension('Spatial') == 'Available':
    arcpy.CheckOutExtension('Spatial')
    arcpy.env.overwriteOutput = True

__all__=['toa_radiance_8',          # complete
         'toa_radiance_457']        # complete


def toa_radiance_8(band_nums, meta_path, outdir = None):
[docs] """ Top of Atmosphere radiance (in Watts/(square meter x steradians x micrometers)) conversion for landsat 8 data. To be performed on raw Landsat 8 level 1 data. See link below for details: see here http://landsat.usgs.gov/Landsat8_Using_Product.php :param band_nums: A list of desired band numbers such as [3, 4, 5] :param meta_path: The full filepath to the metadata file for those bands :param outdir: Output directory to save converted files. :return output_filelist: List of filepaths created by this function. """ meta_path = os.path.abspath(meta_path) output_filelist = [] #enforce list of band numbers and grab the metadata from the MTL file band_nums = core.enf_list(band_nums) band_nums = map(str, band_nums) meta = grab_meta(meta_path) OLI_bands = ['1','2','3','4','5','6','7','8','9'] #loop through each band for band_num in band_nums: if band_num in OLI_bands: #create the band name band_path = meta_path.replace("MTL.txt","B{0}.tif".format(band_num)) Qcal = arcpy.Raster(band_path) null_raster = arcpy.sa.SetNull(Qcal, Qcal, "VALUE = 0") #scrape the attribute data Ml = getattr(meta,"RADIANCE_MULT_BAND_{0}".format(band_num)) # multiplicative scaling factor Al = getattr(meta,"RADIANCE_ADD_BAND_{0}".format(band_num)) # additive rescaling factor #calculate Top-of-Atmosphere radiance TOA_rad = (null_raster * Ml) + Al del null_raster # create the output name and save the TOA radiance tiff if "\\" in meta_path: name = meta_path.split("\\")[-1] elif "//" in meta_path: name = meta_path.split("//")[-1] rad_name = name.replace("_MTL.txt", "_B{0}".format(band_num)) if outdir is not None: outdir = os.path.abspath(outdir) outname = core.create_outname(outdir, rad_name, "TOA_Rad", "tif") else: folder = os.path.split(meta_path)[0] outname = core.create_outname(folder, rad_name, "TOA_Rad", "tif") TOA_rad.save(outname) output_filelist.append(outname) print("Saved toa_radiance at {0}".format(outname)) #if listed band is not a OLI sensor band, skip it and print message else: print("Can only perform reflectance conversion on OLI sensor bands") print("Skipping band {0}".format(band_num)) return output_filelist def toa_radiance_457(band_nums, meta_path, outdir = None):
[docs] """ Top of Atmosphere radiance (in Watts/(square meter x steradians x micrometers)) conversion for Landsat 4, 5, and 7 data. To be performed on raw Landsat 4, 5, or 7 level 1 data. :param band_nums: A list of desired band numbers such as [3, 4, 5] :param meta_path: The full filepath to the metadata file for those bands :param outdir: Output directory to save converted files. :return output_filelist: List of filepaths created by this function. """ output_filelist = [] meta_path = os.path.abspath(meta_path) band_nums = core.enf_list(band_nums) band_nums = map(str, band_nums) #metadata format was changed August 29, 2012. This tool can process either the new or old format f = open(meta_path) MText = f.read() metadata = grab_meta(meta_path) #the presence of a PRODUCT_CREATION_TIME category is used to identify old metadata #if this is not present, the meta data is considered new. #Band6length refers to the length of the Band 6 name string. In the new metadata this string is longer if "PRODUCT_CREATION_TIME" in MText: Meta = "oldMeta" Band6length = 2 else: Meta = "newMeta" Band6length = 8 #The tilename is located using the newMeta/oldMeta indixes and the date of capture is recorded if Meta == "newMeta": TileName = getattr(metadata, "LANDSAT_SCENE_ID") year = TileName[9:13] jday = TileName[13:16] date = getattr(metadata, "DATE_ACQUIRED") elif Meta == "oldMeta": TileName = getattr(metadata, "BAND1_FILE_NAME") year = TileName[13:17] jday = TileName[17:20] date = getattr(metadata, "ACQUISITION_DATE") #the spacecraft from which the imagery was capture is identified #this info determines the solar exoatmospheric irradiance (ESun) for each band spacecraft = getattr(metadata, "SPACECRAFT_ID") if "7" in spacecraft: ESun = (1969.0, 1840.0, 1551.0, 1044.0, 255.700, 0., 82.07, 1368.00) TM_ETM_bands = ['1','2','3','4','5','7','8'] elif "5" in spacecraft: ESun = (1957.0, 1826.0, 1554.0, 1036.0, 215.0, 0. ,80.67) TM_ETM_bands = ['1','2','3','4','5','7'] elif "4" in spacecraft: ESun = (1957.0, 1825.0, 1557.0, 1033.0, 214.9, 0. ,80.72) TM_ETM_bands = ['1','2','3','4','5','7'] else: arcpy.AddError("This tool only works for Landsat 4, 5, or 7") raise arcpy.ExecuteError() #Calculating values for each band for band_num in band_nums: if band_num in TM_ETM_bands: print("Processing Band {0}".format(band_num)) pathname = meta_path.replace("MTL.txt", "B{0}.tif".format(band_num)) Oraster = arcpy.Raster(pathname) null_raster = arcpy.sa.SetNull(Oraster, Oraster, "VALUE = 0") #using the oldMeta/newMeta indixes to pull the min/max for radiance/Digital numbers if Meta == "newMeta": LMax = getattr(metadata, "RADIANCE_MAXIMUM_BAND_{0}".format(band_num)) LMin = getattr(metadata, "RADIANCE_MINIMUM_BAND_{0}".format(band_num)) QCalMax = getattr(metadata, "QUANTIZE_CAL_MAX_BAND_{0}".format(band_num)) QCalMin = getattr(metadata, "QUANTIZE_CAL_MIN_BAND_{0}".format(band_num)) elif Meta == "oldMeta": LMax = getattr(metadata, "LMAX_BAND{0}".format(band_num)) LMin = getattr(metadata, "LMIN_BAND{0}".format(band_num)) QCalMax = getattr(metadata, "QCALMAX_BAND{0}".format(band_num)) QCalMin = getattr(metadata, "QCALMIN_BAND{0}".format(band_num)) Radraster = (((LMax - LMin)/(QCalMax-QCalMin)) * (null_raster - QCalMin)) + LMin Oraster = 0 del null_raster band_rad = "{0}_B{1}".format(TileName, band_num) #create the output name and save the TOA radiance tiff if outdir is not None: outdir = os.path.abspath(outdir) outname = core.create_outname(outdir, band_rad, "TOA_Rad", "tif") else: folder = os.path.split(meta_path)[0] outname = core.create_outname(folder, band_rad, "TOA_Rad", "tif") Radraster.save(outname) output_filelist.append(outname) del Radraster print("toa radiance saved for Band {0}".format(band_num)) #if listed band is not a TM/ETM+ sensor band, skip it and print message else: print("Can only perform reflectance conversion on TM/ETM+ sensor bands") print("Skipping band {0}".format(band_num)) f.close() return output_filelist