python對landsat8數據進行輻射校正和大氣校正

　　Landsat8 L1 T數據是輻射校正數據使用地面控制點和數字高程模型數據進行精確校正后的數據產品，還需要做輻射校正（輻射定標和大氣校正）。

一、輻射定標

　　輻射亮度L=DN*Gain+Bias

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"] = point2lon
    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()

python運行完

可以看出，這里只是處理了可見光和全色波段，雲、熱紅外波段都沒有處理。

文件：

 結果顯示

原始影像

 ENVI大氣校正處理影像

 該代碼大氣校正影像

參考鏈接：

https://blog.csdn.net/weixin_40501429/article/details/115671738

https://blog.csdn.net/weixin_40501429/article/details/115856301