Source code for dnppy.landsat.toa_reflectance

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

__all__=['toa_reflectance_8',       # complete       
         'toa_reflectance_457']     # complete


def toa_reflectance_8(band_nums, meta_path, outdir = None):
[docs] """ Converts Landsat 8 bands to Top-of-Atmosphere reflectance. 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. If left False it will save ouput files in the same directory as input files. :return output_filelist: List of files created by this function """ output_filelist = [] # enforce the list of band numbers and grab metadata from the MTL file band_nums = core.enf_list(band_nums) band_nums = map(str, band_nums) OLI_bands = ['1','2','3','4','5','6','7','8','9'] meta_path = os.path.abspath(meta_path) meta = grab_meta(meta_path) # cycle through each band in the list for calculation, ensuring each is in the list of OLI bands for band_num in band_nums: if band_num in OLI_bands: # scrape data from the given file path and attributes in the MTL file band_path = meta_path.replace("MTL.txt","B{0}.tif".format(band_num)) Qcal = arcpy.Raster(band_path) Mp = getattr(meta,"REFLECTANCE_MULT_BAND_{0}".format(band_num)) # multiplicative scaling factor Ap = getattr(meta,"REFLECTANCE_ADD_BAND_{0}".format(band_num)) # additive rescaling factor SEA = getattr(meta,"SUN_ELEVATION")*(math.pi/180) # sun elevation angle theta_se # get rid of the zero values that show as the black background to avoid skewing values null_raster = arcpy.sa.SetNull(Qcal, Qcal, "VALUE = 0") # calculate top-of-atmosphere reflectance TOA_ref = (((null_raster * Mp) + Ap)/(math.sin(SEA))) # save the data to the automated name if outdir is given or in the parent folder if not if outdir is not None: outdir = os.path.abspath(outdir) outname = core.create_outname(outdir, band_path, "TOA_Ref", "tif") else: folder = os.path.split(meta_path)[0] outname = core.create_outname(folder, band_path, "TOA_Ref", "tif") TOA_ref.save(outname) output_filelist.append(outname) print("Saved output at {0}".format(outname)) # if listed band is not an 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_reflectance_457(band_nums, meta_path, outdir = None):
[docs] """ This function is used to convert Landsat 4, 5, or 7 pixel values from digital numbers to Top-of-Atmosphere Reflectance. To be performed on raw Landsat 4, 5, or 7 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. If left False it will save ouput files in the same directory as input files. :return output_filelist: List of files created by this function """ output_filelist = [] 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() meta_path = os.path.abspath(meta_path) 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() # determing if year is leap year and setting the Days in year accordingly if float(year) % 4 == 0: DIY = 366. else: DIY=365. # using the date to determining the distance from the sun theta = 2 * math.pi * float(jday)/DIY dSun2 = (1.00011 + 0.034221 * math.cos(theta) + 0.001280 * math.sin(theta) + 0.000719 * math.cos(2*theta)+ 0.000077 * math.sin(2 * theta)) SZA = 90. - float(getattr(metadata, "SUN_ELEVATION")) # 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 indices 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 # Calculating temperature for band 6 if present Refraster = (math.pi * Radraster * dSun2) / (ESun[int(band_num[0])-1] * math.cos(SZA*(math.pi/180))) # construc output names for each band based on whether outdir is set (default is False) if outdir is not None: outdir = os.path.abspath(outdir) BandPath = core.create_outname(outdir, pathname, "TOA_Ref", "tif") else: folder = os.path.split(meta_path)[0] BandPath = core.create_outname(folder, pathname, "TOA_Ref", "tif") Refraster.save(BandPath) output_filelist.append(BandPath) del Refraster, Radraster print("Reflectance Calculated 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