from osgeo import gdal
from osgeo import gdal_array
import numpy as np
from show import TwoPercentLinear
from matplotlib import pyplot as plt
import cv2 as cv
 
 
class Landsat8Reader(object):
 
    def __init__(self):
        self.base_path = r"D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\LC08_L1TP_134036_20170808_20170813_01_T1"
        self.bands = 7
        self.band_file_name = []
        self.nan_position = []
 
    def read(self):
        for band in range(self.bands):
            band_name = self.base_path + "_B" + str(band + 1) + ".tif"
            self.band_file_name.append(band_name)
 
        ds = gdal.Open(self.band_file_name[0])
 
        image_dt = ds.GetRasterBand(1).DataType
        image = np.zeros((ds.RasterYSize, ds.RasterXSize, self.bands),
                         dtype=np.float)
        for band in range(self.bands):
            ds = gdal.Open(self.band_file_name[band])
            band_image = ds.GetRasterBand(1)
            image[:, :, band] = band_image.ReadAsArray()
 
        self.nan_position = np.where(image == 0)
        image[self.nan_position] = np.nan
 
        del ds
 
        return image
 
    def write(self, image, file_name, bands, format='GTIFF'):
 
        ds = gdal.Open(self.band_file_name[0])
        projection = ds.GetProjection()
        geotransform = ds.GetGeoTransform()
        x_size = ds.RasterXSize
        y_size = ds.RasterYSize
        del ds
        band_count = bands
        dtype = gdal.GDT_Float32
 
        driver = gdal.GetDriverByName(format)
        new_ds = driver.Create(file_name, x_size, y_size, band_count, dtype)
        new_ds.SetGeoTransform(geotransform)
        new_ds.SetProjection(projection)
 
        for band in range(self.bands):
            new_ds.GetRasterBand(band + 1).WriteArray(image[:, :, band])
            new_ds.FlushCache()
 
        del new_ds
 
    def show_true_color(self, image):
 
        index = np.array([3, 2, 1])
        true_color_image = image[:, :, index].astype(np.float32)
        strech_image = TwoPercentLinear(true_color_image)
        plt.imshow(strech_image)
 
    def show_CIR_color(self, image):
 
        index = np.array([4, 3, 2])
        true_color_image = image[:, :, index].astype(np.float32)
        strech_image = TwoPercentLinear(true_color_image)
        plt.imshow(strech_image)
 
    def radiometric_calibration(self):
        # 辐射定标
        image = self.read()
 
        def get_calibration_parameters():
            filename = self.base_path + "_MTL" + ".txt"
            f = open(filename, 'r')
            # f = open(r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\LC08_L1TP_134036_20170808_20170813_01_T1_MTL.txt', 'r')
            # readlines()方法读取整个文件所有行，保存在一个列表(list)变量中，每行作为一个元素，但读取大文件会比较占内存。
            metadata = f.readlines()
            f.close()
            multi_parameters = []
            add_parameters = []
            parameter_start_line = 0
            for lines in metadata:
                # split()方法通过指定分隔符对字符串进行切片，如果参数 num 有指定值，则分隔 num+1 个子字符串。
                # 无指定值，默认为 -1, 即分隔所有符合要求的。
                test_line = lines.split('=')
                if test_line[0] == '    RADIANCE_MULT_BAND_1 ':
                    break
                else:
                    parameter_start_line = parameter_start_line + 1
 
            for lines in range(parameter_start_line, parameter_start_line + 11):
                parameter = float(metadata[lines].split('=')[1])
                multi_parameters.append(parameter)
 
            # 由于RADIANCE_MULT_BAND参数和RADIANCE_ADD_BAND参数是挨着的，所以直接+11个参数即可
            for lines in range(parameter_start_line + 11, parameter_start_line + 22):
                parameter = float(metadata[lines].split('=')[1])
                add_parameters.append(parameter)
 
            return multi_parameters, add_parameters
 
        multi_parameters, add_parameters = get_calibration_parameters()
        cali_image = np.zeros_like(image, dtype=np.float32)
 
        for band in range(self.bands):
            gain = multi_parameters[band]
            offset = add_parameters[band]
            # 辐射定标像元 = DN * 增益 + 偏置，线性关系。
            cali_image[:, :, band] = image[:, :, band] * gain + offset
 
        del image
        return cali_image
 
 
if __name__ == "__main__":
    reader = Landsat8Reader()
    image = reader.read()
    cali_image = reader.radiometric_calibration()
    file_path = r'D:\ProfessionalProfile\LandsatImage\1_RadiationCalibration0414\LC08.tif'
    reader.write(cali_image, file_path, reader.bands)

结果

二、大气校正

　　表观反射率ρ=π*L*D²/（ESUN*cosθ）式中，ρ为大气层顶( TOA) 表观反射率( 无量纲)，π为常量( 球面度sr)，L 为大气层顶进人卫星传感器的光谱辐射亮度（ W*1/（m²*sr*um）），D为日地之间距离( 天文单位)， ESUN为大气层顶的平均太阳光谱辐照度（ W*1/（m²*sr*um））θ为太阳的天顶角。

代码

# glob模块参考： https://blog.csdn.net/gufenchen/article/details/90723418
# glob是python自带的一个操作文件的相关模。用它可以查找符合特定规则的文件路径名。使用该模块查找文件，只需要用到： “*”, “?”, “[]”这三个匹配符;
import glob
import os
import sys
import tarfile
import re
from osgeo import gdal
import numpy
from Py6S import *
from osgeo import gdal
import pdb
import shutil
# argparse是一个Python模块：命令行选项、参数和子命令解析器。
import argparse
# 从base.py文件导入MeanDEM函数
from base import MeanDEM
 
 
def parse_arguments(argv):
 
    # 此部分参考： https://blog.csdn.net/qq_36653505/article/details/83788460
    # 此部分参考： https://blog.csdn.net/Rocky6688/article/details/104295574
    # 使用argparse的第一步是创建一个ArgumentParser对象。ArgumentParser对象包含将命令行解析成Python数据类型所需的全部信息。
    parser = argparse.ArgumentParser()    # 创建一个解析对象
 
    # parser.add_argument() 向该对象中添加你要关注的命令行参数和选项
    # name or flags - 一个命名或者一个选项字符串的列表；type - 命令行参数应当被转换成的类型；
    # help - 一个此选项作用的简单描述；default - 当参数未在命令行中出现时使用的值。
    # parser.add_argument('--Input_dir', type=str, help=r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\.', default=None)
    # parser.add_argument('--Output_dir', type=str, help=r'D:\ProfessionalProfile\LandsatImage\1_OutputPy-RadiationCalibration0414\.', default=None)
 
    parser.add_argument('--Input_dir', type=str, help='Input dir', default=r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1')
    parser.add_argument('--Output_dir', type=str, help='Output dir', default=r'D:\ProfessionalProfile\LandsatImage\1_OutputPy-RadiationCalibration0414')
 
    # parser.add_argument('--Input_dir', type=str, help='Input dir', default=None)
    # parser.add_argument('--Output_dir', type=str, help='Output dir', default=None)
 
    # parse_args(args=None, nampespace=None)，args 参数名称，namespace 赋值。
    return parser.parse_args(argv)    # 进行解析
 
 
# 逐波段辐射定标
def RadiometricCalibration(BandId):
    # LandSat8 TM辐射定标参数
    global data2,ImgRasterData
    parameter_OLI = numpy.zeros((9,2))
 
    #计算辐射亮度参数
    parameter_OLI[0,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_1.+',data2)[0]).split("=")[1])
    parameter_OLI[1,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_2.+',data2)).split("=")[1])
    parameter_OLI[2,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_3.+',data2)).split("=")[1])
    parameter_OLI[3,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_4.+',data2)).split("=")[1])
    parameter_OLI[4,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_5.+',data2)).split("=")[1])
    parameter_OLI[5,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_6.+',data2)).split("=")[1])
    parameter_OLI[6,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_7.+',data2)).split("=")[1])
    parameter_OLI[7,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_8.+',data2)).split("=")[1])
    parameter_OLI[8,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_9.+',data2)).split("=")[1])
 
    parameter_OLI[0,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_1.+',data2)[0]).split("=")[1])
    parameter_OLI[1,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_2.+',data2)).split("=")[1])
    parameter_OLI[2,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_3.+',data2)).split("=")[1])
    parameter_OLI[3,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_4.+',data2)).split("=")[1])
    parameter_OLI[4,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_5.+',data2)).split("=")[1])
    parameter_OLI[5,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_6.+',data2)).split("=")[1])
    parameter_OLI[6,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_7.+',data2)).split("=")[1])
    parameter_OLI[7,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_8.+',data2)).split("=")[1])
    parameter_OLI[8,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_9.+',data2)).split("=")[1])
 
    Gain = parameter_OLI[int(BandId) - 1,0]
    Bias = parameter_OLI[int(BandId) - 1,1]
 
    RaCal = numpy.where(ImgRasterData>0 ,Gain * ImgRasterData + Bias,-9999)
    return (RaCal)
 
# 6s大气校正
def AtmosphericCorrection(BandId):
 
    global data
    # 6S模型
    s = SixS()
 
    s.geometry = Geometry.User()
    s.geometry.solar_z = 90-float(''.join(re.findall('SUN_ELEVATION.+',data2)).split("=")[1])
    s.geometry.solar_a = float(''.join(re.findall('SUN_AZIMUTH.+',data2)).split("=")[1])
    s.geometry.view_z = 0
    s.geometry.view_a = 0
 
 
    # 日期
    Dateparm = ''.join(re.findall('DATE_ACQUIRED.+',data2)).split("=")
    Date = Dateparm[1].split('-')
 
    s.geometry.month = int(Date[1])
    s.geometry.day = int(Date[2])
 
 
    # 中心经纬度
    point1lat = float(''.join(re.findall('CORNER_UL_LAT_PRODUCT.+',data2)).split("=")[1])
    point1lon = float(''.join(re.findall('CORNER_UL_LON_PRODUCT.+',data2)).split("=")[1])
    point2lat = float(''.join(re.findall('CORNER_UR_LAT_PRODUCT.+',data2)).split("=")[1])
    point2lon = float(''.join(re.findall('CORNER_UR_LON_PRODUCT.+',data2)).split("=")[1])
    point3lat = float(''.join(re.findall('CORNER_LL_LAT_PRODUCT.+',data2)).split("=")[1])
    point3lon = float(''.join(re.findall('CORNER_LL_LON_PRODUCT.+',data2)).split("=")[1])
    point4lat = float(''.join(re.findall('CORNER_LR_LAT_PRODUCT.+',data2)).split("=")[1])
    point4lon = float(''.join(re.findall('CORNER_LR_LON_PRODUCT.+',data2)).split("=")[1])
 
    sLongitude = (point1lon + point2lon + point3lon + point4lon) / 4
    sLatitude = (point1lat + point2lat + point3lat + point4lat) / 4
 
    # 大气模式类型
    if sLatitude > -15 and sLatitude <= 15:
        s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)
 
    if sLatitude > 15 and sLatitude <= 45:
        if s.geometry.month > 4 and s.geometry.month <= 9:
            s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
        else:
            s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)
 
    if sLatitude > 45 and sLatitude <= 60:
        if s.geometry.month > 4 and s.geometry.month <= 9:
            s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)
        else:
            s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)
 
    # 气溶胶类型大陆
    s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental)
 
    # 目标地物？？？？？？
    s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36)
 
    # 550nm气溶胶光学厚度,根据日期从MODIS处获取。
    #s.visibility=40.0
    s.aot550 = 0.14497
 
    # 通过研究去区的范围去求DEM高度。
    pointUL = dict()
    pointDR = dict()
    pointUL["lat"] = point1lat
    pointUL["lon"] = point1lon
    pointDR["lat"] = point4lat
    pointDR["lon"] = point4lon
    meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001
 
    # 研究区海拔、卫星传感器轨道高度
    s.altitudes = Altitudes()
    s.altitudes.set_target_custom_altitude(meanDEM)
    s.altitudes.set_sensor_satellite_level()
 
    # 校正波段（根据波段名称）
    if BandId == '1':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)
 
    elif BandId == '2':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B2)
 
    elif BandId == '3':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B3)
 
    elif BandId == '4':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B4)
 
    elif BandId == '5':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B5)
 
    elif BandId == '6':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B6)
 
    elif BandId == '7':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B7)
 
    elif BandId == '8':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B8)
 
    elif BandId == '9':
        s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B9)
 
    # 下垫面非均一、朗伯体
    s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1)
 
    # 运行6s大气模型
    s.run()
 
    xa = s.outputs.coef_xa
    xb = s.outputs.coef_xb
    xc = s.outputs.coef_xc
    x = s.outputs.values
    return (xa, xb, xc)
 
 
if __name__ == '__main__':
 
    #输入数据路径
    # sys.argv参考：https://blog.csdn.net/u013044310/article/details/79499677
    # sys.argv[]说白了就是一个从程序外部获取参数的桥梁，这个“外部”很关键，所以那些试图从代码来说明它作用的解释一直没看明白。因为我们从外部取得的参数可以是多个，所以获得的是一个列表（list)，也就是说sys.argv其实可以看作是一个列表，所以才能用[]提取其中的元素。其第一个元素是程序本身，随后才依次是外部给予的参数。
    RootInputPath = parse_arguments(sys.argv[1:]).Input_dir
    RootOutName = parse_arguments(sys.argv[2:]).Output_dir
 
    #创建日志文件
    LogFile = open(os.path.join(RootOutName,'log.txt'),'w')
 
    # Python os.walk() 方法  https://www.runoob.com/python/os-walk.html
    # os.walk(top[, topdown=True[, onerror=None[, followlinks=False]]])
    # top -- 是你所要遍历的目录的地址, 返回的是一个三元组(root,dirs,files)。
    # root 所指的是当前正在遍历的这个文件夹的本身的地址
    # dirs 是一个 list ，内容是该文件夹中所有的目录的名字(不包括子目录)
    # files 同样是 list , 内容是该文件夹中所有的文件(不包括子目录)
    for root,dirs,RSFiles in os.walk(RootInputPath):
 
        #判断是否进入最底层
        if len(dirs)==0:
            #根据输入输出路径建立生成新文件的路径
            RootInputPathList = RootInputPath.split(os.path.sep)
            RootList = root.split(os.path.sep)
            StartList = len(RootInputPathList)
            EndList = len(RootList)
            outname = RootOutName
            for i in range(StartList,EndList):
                if os.path.exists(os.path.join(outname,RootList[i]))==False:
                    os.makedirs(os.path.join(outname,RootList[i]))
                    outname=os.path.join(outname,RootList[i])
                else:
                    outname=os.path.join(outname,RootList[i])
 
            MeteDatas = glob.glob(os.path.join(root,'*MTL.txt'))
 
            for MeteData in MeteDatas:
                pass
 
            f = open(MeteData)
            data = f.readlines()
            data2 =' '.join(data)
            
            shutil.copyfile(MeteData,os.path.join(outname,os.path.basename(MeteData)))
 
            if len(os.path.basename(MeteData))<10:
 
                RSbands = glob.glob(os.path.join(root,"B0[1-8].tiff"))
            else:
                RSbands = glob.glob(os.path.join(root,"*B[1-8].TIF"))
            print('影像'+root+'开始大气校正')
            print(RSbands)
            for tifFile in RSbands:
 
                BandId = (os.path.basename(tifFile).split('.')[0])[-1]
 
                #捕捉打开数据出错异常
                try:
                    IDataSet = gdal.Open(tifFile)
                except Exception as e:
                    print("文件%S打开失败" % tifFile)
                    LogFile.write('\n'+tifFile+'数据打开失败')
 
                if IDataSet == None:
                    LogFile.write('\n'+tifFile+'数据集读取为空')
                    continue
                else:
                    #获取行列号
                    cols = IDataSet.RasterXSize
                    rows = IDataSet.RasterYSize
                    ImgBand = IDataSet.GetRasterBand(1)
                    ImgRasterData = ImgBand.ReadAsArray(0, 0, cols, rows)
 
                    if ImgRasterData is None:
                        LogFile.write('\n'+tifFile+'栅格数据为空')
                        continue
                    else:
                        #设置输出文件路径
                        outFilename=os.path.join(outname,os.path.basename(tifFile))
 
                        #如果文件存在就跳过，进行下一波段操作
                        if os.path.isfile(outFilename):
                            print("%s已经完成" % outFilename)
                            continue
                        else:
                            # #辐射校正
                            RaCalRaster = RadiometricCalibration(BandId)
                            #大气校正
                            a, b, c = AtmosphericCorrection(BandId)
                            y = numpy.where(RaCalRaster!=-9999,a * RaCalRaster - b,-9999)
                            atc = numpy.where(y!=-9999,(y / (1 + y * c))*10000,-9999)
                            
                            driver = IDataSet.GetDriver()
                            #输出栅格数据集
                            outDataset = driver.Create(outFilename, cols, rows, 1, gdal.GDT_Int16)
 
                            # 设置投影信息，与原数据一样
                            geoTransform = IDataSet.GetGeoTransform()
                            outDataset.SetGeoTransform(geoTransform)
                            proj = IDataSet.GetProjection()
                            outDataset.SetProjection(proj)
 
                            outband = outDataset.GetRasterBand(1)
                            outband.SetNoDataValue(-9999)
                            outband.WriteArray(atc, 0, 0)
                print('第%s波段计算完成'%BandId)
            # print(root+'计算完成')
            print('\n')
    #关闭日志文件
    LogFile.close()