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[, one rror=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