python对landsat8数据进行辐射校正和大气校正

　　Landsat8 L1 T数据是辐射校正数据使⽤地⾯控制点和数字⾼程模型数据进⾏精确校正后的数据产品，还需要做辐射校正（辐射定标和⼤⽓校正）。

⼀、辐射定标

　　辐射亮度L=DN*Gain+Bias

from osgeo import gdal

from osgeo import gdal_arrayimport numpy as np

from show import TwoPercentLinearfrom matplotlib import pyplot as pltimport 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 globimport osimport sysimport tarfileimport re

from osgeo import gdalimport numpy

from Py6S import *from osgeo import gdalimport pdbimport 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⼤⽓校正处理影像

该代码⼤⽓校正影像

参考链接：