4、利用python实现Landsat系列数据(4,5,8)影像波段星上反射率(TOA reflectance)定标。

   前面已经实现了python对Landsat解压后的元数据中各种关于影像参数读取，这里进一步利用这些参数实现波段的星上反射率定标。下面的代码也是在dnppy中相关代码作的一些修改，可酌情参考。

# standard imports
from lst_meta import landsat_metadata
import math
import os
from osgeo import gdal
import numpy as np
def find_last(string,str):
    last_position=-1
    while True:
        position=string.find(str,last_position+1)
        if position==-1:
            return last_position
        last_position=position
def create_outname(outdir, inname, suffix, ext = False):
    """
    Quick way to create unique output file names within iterative functions
    This function is built to simplify the creation of output file names. Function allows
    ``outdir = False`` and will create an outname in the same directory as inname. Function will
    add a the user input suffix, separated by an underscore "_" to generate an output name.
    this is useful when performing a small modification to a file and saving new output with
    a new suffix. Function merely returns an output name, it does not save the file as that name.
    :param outdir:      either the directory of the desired outname or False to create an outname
                        in the same directory as the inname
    :param inname:      the input file from which to generate the output name "outname"
    :param suffix:      suffix to attach to the end of the filename to mark it as output
    :param ext:         specify the file extension of the output filename. Leave blank or False
                        and the outname will inherit the same extension as inname.
    :return outname:    the full filepath at which a new file can be created.
    """
    # isolate the filename from its directory and extension
    if os.path.isfile(inname):
        head, tail = os.path.split(inname)
        noext = tail.split('.')[:-1]
        noext = '.'.join(noext)
    else:
        head = ""
        tail = inname
        if "." in inname:
            noext = tail.split('.')[:-1]
            noext = '.'.join(noext)
        else:
            noext = inname
    # create the suffix
    if ext:
        suffix = "_{0}.{1}".format(suffix, ext)
    else:
        ext = tail.split('.')[-1:]
        suffix = "_{0}.{1}".format(suffix, ''.join(ext))
    if outdir:
        outname = os.path.join(outdir, noext + suffix)
        return outname
    else:
        outname = os.path.join(head, noext + suffix)
        return outname
    
    
def enf_list(item):
    """
    When a list is expected, this function can be used to ensure
    non-list data types are placed inside of a single entry list.
    :param item:    any datatype
    :return list:   a list type
    """
    if not isinstance(item,list) and item:
        return [item]
    else:
        return item
def toa_reflectance_8(band_num, meta_path, ds_name = None, outdir = None):
    """
    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_num:    A desired band numbers such as 3
    :param ds_name:     The name of input image file that used for calibrating TOA reflectance, such as 'd:/band3.tif'
    :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
    """
    band_num = str(band_num)
    OLI_bands = ['1','2','3','4','5','6','7','8','9']
    meta_path = os.path.abspath(meta_path)
    meta = landsat_metadata(meta_path)
    if band_num in OLI_bands:
    # scrape data from the given file path and attributes in the MTL file
        if ds_name == None:
            band_path = meta_path.replace("MTL.txt","B{0}.tif".format(band_num))
        else:
    #        pos = max(find_last(meta_path, '\\'), find_last(meta_path, '/'))
            band_path = ds_name
        inds = gdal.Open(band_path)
        Qcal = inds.ReadAsArray()
        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 = np.where(Qcal!=0, Qcal, np.NaN)
        # calculate top-of-atmosphere reflectance
        TOA_ref = (((null_raster * Mp) + Ap)/(math.sin(SEA)))*10000
        # 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 = create_outname(outdir, band_path, "TOA_Ref", "tif")
        else:
            folder = os.path.split(meta_path)[0]
            outname = create_outname(folder, band_path, "TOA_Ref", "tif")
        # save output
        driver = inds.GetDriver()
        outimg = driver.Create(outname, TOA_ref.shape[1], TOA_ref.shape[0], 1, gdal.GDT_Int16)
        outimg.SetGeoTransform(inds.GetGeoTransform())
        outimg.SetProjection(inds.GetProjection())
        outBand = outimg.GetRasterBand(1)
        outBand.WriteArray(TOA_ref)
        outBand.SetNoDataValue(np.NaN)
        outBand.FlushCache
        outimg = None
        
        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))
    print outname
    return outname
def toa_reflectance_457(band_num, meta_path, ds_name = None, outdir = None):
    """
    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_num:    A desired band numbers such as 3
    :param ds_name:     The name of input image file that used for calibrating TOA reflectance, such as 'd:/band3.tif'
    :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
    """
   
    band_num = str(band_num)
    # 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)
#    print meta_path
    metadata = landsat_metadata(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:
        print "This tool only works for Landsat 4, 5, or 7"
        raise
    # 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
    if band_num in TM_ETM_bands:
        print("Processing Band {0}".format(band_num))
        if ds_name == None:
            pathname = meta_path.replace("MTL.txt","B{0}.tif".format(band_num))
        else:
    #       pos = max(find_last(meta_path, '\\'), find_last(meta_path, '/'))
            pathname = ds_name
        inds = gdal.Open(pathname)
        Oraster = inds.ReadAsArray()
     
        null_raster = np.where(Oraster!=0, Oraster, np.NaN)
        # 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 = 10000*(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 = create_outname(outdir, pathname, "TOA_Ref", "tif")
        else:
            folder = os.path.split(meta_path)[0]
            BandPath = create_outname(folder, pathname, "TOA_Ref", "tif")
         # save output
        driver = inds.GetDriver()
        outimg = driver.Create(BandPath, Refraster.shape[1], Refraster.shape[0], 1, gdal.GDT_Int16)
        outimg.SetGeoTransform(inds.GetGeoTransform())
        outimg.SetProjection(inds.GetProjection())
        outBand = outimg.GetRasterBand(1)
        outBand.WriteArray(Refraster)
        outBand.SetNoDataValue(np.NaN)
        outBand.FlushCache
        outimg = None
        
        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()
    print BandPath
    return BandPath

5、星上反射率定标示例

import lst_toa_ref   #假如上面的代码保存为"lst_toa_ref.py"

meta_path = 'd:/test/LT51240392010355BKT00_MTL.txt'

lst_toa_ref.toa_reflectance_457(4,meta_path)  #对landsat5近红外波段的星上反射率定标（*10000）

运行后，结果展示：