摘要: 本文主要介绍如何使用Python对Landsat8-9影像在本地进行批量预处理(辐射定标及波段合成)。是笔者经过对论文、数据手册和博客等一系列学习资料的一个总结。
文章目录
前言
Landsat系列影像作为遥感领域使用最为广泛、时间序列最长的和影像质量最好的免费遥感数据之一。在科学实践的过程中,往往需要同时使用多张Landsat影像进行科学分析,对于原始数据的预处理过程(辐射定标和波段合成)成为定量分析前必不可少的步骤。
传统本地预处理方法一般采用ENVI或IDL的方法对数据进行预处理,但是:(1)老版本的ENVI对于卫星最新处理结果(Landsat Collection2)的兼容性较差(需要动手调整头文件),且面对众多需要预处理数据时,操作较多耗时较长;(2)IDL是针对遥感数据处理的专用语言,有一定的学习成本。
以下将介绍如何使用Python对Landsat影像在本地进行批量辐射定标和波段合成的预处理工作。
一、概念辨析
一般来讲,辐射定标是指,将遥感图像的数字量化值转化为辐射亮度、反射率或表面温度等物理量的处理过程。其中:
数字量化值:又称DN值,是遥感影像数据的原始值,数值类型一般为整数。(Landsat8-9有效范围为:1-65455)
辐射亮度:Radiance,是某一个面积辐射能量的总和。(有单位,如:
µ
W
/
(
c
m
2
∗
s
r
∗
n
m
)
µW/(cm2*sr*nm)
µW/(cm2∗sr∗nm))
反射率:Reflectance,反射率是物体表面所能反射的辐射量和它所接受的辐射量的比值,一般在[0,1]范围。
表面温度:分为大气表观或地表温度。
TOA:Top Of Atmosphere,代表大气表观(反射率/温度)
关于Landsat系列卫星的分级和命名标准,参见:https://blog.csdn.net/Mrwindccc/article/details/125407400
二、具体实现
1.实现步骤概要
为了实现python自动Landsat影像数据的定标及波段合成工作,将总体工作分解为以下几个部分:
|----预处理及准备
||--------引用python库
||--------解压Landsat数据压缩包
||--------读取MLT文件
||--------读取遥感图像
|----Landsat C2 L1数据定标
||--------读取L1定标数据
||--------定标为TOA辐射亮度
||--------定标为TOA反射率
|----Landsat C2 L2数据定标
||--------读取L2定标数据
||--------定标为地表反射率
|----波段合成
2.预处理及准备
2.1 引用库
本方法引用了4个常见的py库,分别是os库、tarfile库、numpy库和gdal库。已经有大量的教程介绍如何安装这些库,此处不再赘述。
import tarfile
import os
import numpy as np
from osgeo import gdal
2.2 解压Landsat数据压缩包
def untar(tar_file):
'''
解包.tar文件
@ tar_file : .tar文件路径
return 解包文件夹
'''
untar_dir = tar_file.split('.')[0]
with tarfile.open(tar_file, 'r') as t:
t.extractall(untar_dir)
return untar_dir
2.3 读取MLT文件
def read_MTL(untar_dir):
"""
读取MLT文件
@ untar_dir : 解压产生的文件夹路径
"""
data_list = os.listdir(untar_dir)
for i in data_list:
if 'MTL.txt' in i:
MTL = i
break
MTL = os.path.join(untar_dir, MTL)
with open(MTL, 'r', encoding='utf-8') as f:
MTL_data = f.readlines()
return MTL_data
2.4 读取图像
def read_img(img_filename):
"""
读取影像信息
输入
输出
@im_proj : 投影属性
@im_geotrans : 仿射属性
@im_data : 图像数组
"""
dataset = gdal.Open(img_filename) # 打开文件
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
# print("width:",im_width)
# print("height:",im_height)
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵
im_proj = dataset.GetProjection() # 地图投影信息
im_data = dataset.ReadAsArray(0, 0, im_width, im_height) # 将数据写成数组,对应栅格矩阵
del dataset
return im_proj, im_geotrans, im_data
3.Landsat C2 L1数据定标
3.1 读取L1定标数据
根据MLT文件指标解释,以及L1数据所需的定标参数(分别是’RADIANCE_MULT_BAND’, ‘RADIANCE_ADD_BAND’,‘REFLECTANCE_MULT_BAND’,‘REFLECTANCE_ADD_BAND’),对L1数据的MLT文件进行参数提取。
def get_L1_MTL_info(MTL_data):
"""
@ MTL_data : 获取的MTL数据
"""
info_dict = {'FILE_NAME_BAND': {},
'RADIANCE_MULT_BAND':{},
'RADIANCE_ADD_BAND':{},
'REFLECTANCE_MULT_BAND': {},
'REFLECTANCE_ADD_BAND': {},
}
# FILE_NAME_BAND & REFLECTANCE_MULT_BAND & REFLECTANCE_ADD_BAND
for i in MTL_data:
if 'FILE_NAME_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['FILE_NAME_BAND'][key] = value
elif 'RADIANCE_MULT_BAND' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['RADIANCE_MULT_BAND'][key] = float(value)
elif 'RADIANCE_ADD_BAND' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['RADIANCE_ADD_BAND'][key] = float(value)
elif 'REFLECTANCE_MULT_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_MULT_BAND'][key] = float(value)
elif 'REFLECTANCE_ADD_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_ADD_BAND'][key] = float(value)
elif 'SUN_ELEVATION' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict[key] = float(value)
return info_dict
3.2 定标为TOA辐射亮度
Landsat L1 数据可以使用MTL文件中的辐射重新缩放因子通过以下公式转换为TOA辐射亮度:
L
λ
=
M
L
Q
c
a
l
+
A
L
L_{\lambda}=M_{L}Q_{cal}+A_{L}
Lλ=MLQcal+AL
其中:
L
λ
L_{\lambda}
Lλ = TOA辐射亮度
M
L
M_{L}
ML =MLT文件中的乘法重缩放因子(RADIANCE_MULT_BAND_x,其中x是波段号)
A
L
A_{L}
AL =MLT文件中的加法重缩放因子 (RADIANCE_ADD_BAND_x, 其中x是波段号)
Q
c
a
l
Q_{cal}
Qcal =量化和校准的标准产品像素值(DN值)
def L1OneLayerDN2Radiance(info_dict, band, untar_dir):
'''
@ info_dict : 定标信息
@ band : 定标波段
@ untar_dir : 解包文件夹
'''
filename = os.path.join(untar_dir, info_dict['FILE_NAME_BAND']['FILE_NAME_BAND_' + str(band)])
REFLECTANCE_MULT_BAND = info_dict['RADIANCE_MULT_BAND']['RADIANCE_MULT_BAND_' + str(band)]
REFLECTANCE_ADD_BAND = info_dict['RADIANCE_ADD_BAND']['RADIANCE_ADD_BAND_' + str(band)]
im_proj, im_geotrans, im_data = read_img(filename)
im_data = np.where(im_data, im_data, np.nan)
im_data = (REFLECTANCE_MULT_BAND * im_data.astype(np.float32) + REFLECTANCE_ADD_BAND)
return im_proj, im_geotrans, im_data
3.3 定标为TOA反射率
可以使用MTL文件中的重新缩放系数通过以下公式将DN值转换为TOA反射率:
ρ
λ
=
(
M
ρ
⋅
Q
c
a
l
+
A
ρ
)
/
sin
(
θ
S
E
)
{\rho}_{\lambda}=(M_{\rho}{\cdot}Q_{cal}+A_{\rho})/\sin({\theta}_{SE})
ρλ=(Mρ⋅Qcal+Aρ)/sin(θSE)
其中:
ρ
λ
{\rho}_{\lambda}
ρλ=TOA反射率
M
ρ
M_{\rho}
Mρ=MLT文件中的乘法重缩放因子(REFLECTANCE_MULT_BAND_x,其中x是波段号)
A
ρ
A_{\rho}
Aρ=MLT文件中的乘法重缩放因子(REFLECTANCE_ADD_BAND_x,其中x是波段号)
Q
c
a
l
Q_{cal}
Qcal=量化和校准的标准产品像素值(DN值)
θ
S
E
{\theta}_{SE}
θSE=局部太阳仰角。MLT中提供了以度为单位的场景中心太阳仰角(SUN_ELEVATION)
def L1OneLayerDN2Reflectance(info_dict, band, untar_dir):
'''
@ info_dict : 定标信息
@ band : 定标波段
@ untar_dir : 解包文件夹
'''
filename = os.path.join(untar_dir, info_dict['FILE_NAME_BAND']['FILE_NAME_BAND_' + str(band)])
REFLECTANCE_MULT_BAND = info_dict['REFLECTANCE_MULT_BAND']['REFLECTANCE_MULT_BAND_' + str(band)]
REFLECTANCE_ADD_BAND = info_dict['REFLECTANCE_ADD_BAND']['REFLECTANCE_ADD_BAND_' + str(band)]
SUN_ELEVATION = info_dict['SUN_ELEVATION']*np.pi/180.0
im_proj, im_geotrans, im_data = self.read_img(filename)
im_data = np.where(im_data, im_data, np.nan)
im_data = (REFLECTANCE_MULT_BAND * im_data.astype(np.float32) + REFLECTANCE_ADD_BAND)/np.sin(SUN_ELEVATION)
return im_proj, im_geotrans, im_data
4.Landsat C2 L2数据(Science Product)
4.1 读取L2定标数据
根据MLT文件指标解释,以及L2数据所需的定标参数(分别是’REFLECTANCE_MULT_BAND’,‘REFLECTANCE_ADD_BAND’),对L2数据的MLT文件进行参数提取。
def get_L2_MTL_info(MTL_data):
"""
@ level : 处理的数据等级
@ MTL_data : 获取的MTL数据
"""
info_dict = {'FILE_NAME_BAND': {},
'REFLECTANCE_MULT_BAND': {},
'REFLECTANCE_ADD_BAND': {}}
# FILE_NAME_BAND & REFLECTANCE_MULT_BAND & REFLECTANCE_ADD_BAND
for i in MTL_data:
if 'FILE_NAME_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['FILE_NAME_BAND'][key] = value
elif 'REFLECTANCE_MULT_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_MULT_BAND'][key] = float(value)
elif 'REFLECTANCE_ADD_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_ADD_BAND'][key] = float(value)
elif 'LEVEL1' in i:
break
return info_dict
4.2 定标为地表反射率
可以使用MTL文件中的重新缩放系数通过以下公式将DN值转换为地表反射率:
L
λ
=
M
L
Q
c
a
l
+
A
L
L_{\lambda}=M_{L}Q_{cal}+A_{L}
Lλ=MLQcal+AL
其中:
L
λ
L_{\lambda}
Lλ = TOA辐射亮度
M
L
M_{L}
ML =MLT文件中的乘法重缩放因子(RADIANCE_MULT_BAND_x,其中x是波段号)
A
L
A_{L}
AL =MLT文件中的加法重缩放因子 (RADIANCE_ADD_BAND_x, 其中x是波段号)
Q
c
a
l
Q_{cal}
Qcal =量化和校准的标准产品像素值(DN值)
注:在L2数据中
M
L
M_{L}
ML和
A
L
A_{L}
AL为常数分别是0.0000275和-0.2,这里保持统一也从MLT文件中读取,也可以直接设置为0.0000275和-0.2,不读取MLT。
def L2OneLayerDN2Reflectance(info_dict, band, untar_dir):
'''
@ info_dict : 定标信息
@ band : 定标波段
@ untar_dir : 解包文件夹
'''
filename = os.path.join(untar_dir, info_dict['FILE_NAME_BAND']['FILE_NAME_BAND_' + str(band)])
REFLECTANCE_MULT_BAND = info_dict['REFLECTANCE_MULT_BAND']['REFLECTANCE_MULT_BAND_' + str(band)]
REFLECTANCE_ADD_BAND = info_dict['REFLECTANCE_ADD_BAND']['REFLECTANCE_ADD_BAND_' + str(band)]
# print(filename)
im_proj, im_geotrans, im_data = self.read_img(filename)
im_data = np.where(im_data, im_data, np.nan)
im_data = REFLECTANCE_MULT_BAND * im_data.astype(np.float32) + REFLECTANCE_ADD_BAND
# 以上代码也可以写成
# im_data = 0.0000275 * im_data.astype(np.float32) -0.2
return im_proj, im_geotrans, im_data
5.波段合成
此工作无科学内涵,具体实现参考代码汇总中,函数save_result部分。
代码汇总
import tarfile
import os
import numpy as np
from osgeo import gdal
class Radiometric_Calibration:
def __init__(self):
pass
def untar(self,tar_file):
'''
解包.tar文件
@ tar_file : .tar文件路径
return 解包文件夹
'''
untar_dir = tar_file.split('.')[0]
with tarfile.open(tar_file, 'r') as t:
t.extractall(untar_dir)
return untar_dir
def remove_dir(self,untar_dir):
'''
删除解压产成的文件
@ untar_dir : 解压产生的文件夹路径
'''
data_list = os.listdir(untar_dir)
for i in data_list:
os.remove(os.path.join(untar_dir, i))
os.rmdir(untar_dir)
print(untar_dir + "删除成功")
def read_MTL(self,untar_dir):
"""
读取MLT文件
@ untar_dir : 解压产生的文件夹路径
"""
data_list = os.listdir(untar_dir)
for i in data_list:
if 'MTL.txt' in i:
MTL = i
break
MTL = os.path.join(untar_dir, MTL)
with open(MTL, 'r', encoding='utf-8') as f:
MTL_data = f.readlines()
return MTL_data
def mkdir(self,path):
folder = os.path.exists(path)
if not folder:
os.makedirs(path)
print('文件夹创建成功:', path)
else:
print('文件夹已经存在:', path)
def create_filefolder(self,work_place):
"""
@ work_place : 原始文件夹
"""
dirname, oname = os.path.split(work_place)
result_folder = os.path.join(dirname, oname + '_result')
self.mkdir(result_folder)
return result_folder
def get_L2_MTL_info(self,MTL_data):
"""
@ level : 处理的数据等级
@ MTL_data : 获取的MTL数据
"""
info_dict = {'FILE_NAME_BAND': {},
'REFLECTANCE_MULT_BAND': {},
'REFLECTANCE_ADD_BAND': {}}
for i in MTL_data:
if 'FILE_NAME_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['FILE_NAME_BAND'][key] = value
elif 'REFLECTANCE_MULT_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_MULT_BAND'][key] = float(value)
elif 'REFLECTANCE_ADD_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_ADD_BAND'][key] = float(value)
elif 'LEVEL1' in i:
break
return info_dict
def get_L1_MTL_info(self,MTL_data):
"""
@ level : 处理的数据等级
@ MTL_data : 获取的MTL数据
"""
info_dict = {'FILE_NAME_BAND': {},
'RADIANCE_MULT_BAND':{},
'RADIANCE_ADD_BAND':{},
'REFLECTANCE_MULT_BAND': {},
'REFLECTANCE_ADD_BAND': {},
}
for i in MTL_data:
if 'FILE_NAME_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['FILE_NAME_BAND'][key] = value
elif 'RADIANCE_MULT_BAND' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['RADIANCE_MULT_BAND'][key] = float(value)
elif 'RADIANCE_ADD_BAND' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['RADIANCE_ADD_BAND'][key] = float(value)
elif 'REFLECTANCE_MULT_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_MULT_BAND'][key] = float(value)
elif 'REFLECTANCE_ADD_BAND_' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict['REFLECTANCE_ADD_BAND'][key] = float(value)
elif 'SUN_ELEVATION' in i:
key, value = i.split(' = ')
key = key.replace(' ', '')
value = value.replace('"', '').replace('\n', '')
info_dict[key] = float(value)
return info_dict
def read_img(self,img_filename):
"""
读取影像信息
return :
@im_proj : 投影属性
@im_geotrans : 仿射属性
@im_data : 图像数组
"""
dataset = gdal.Open(img_filename) # 打开文件
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
# print("width:",im_width)
# print("height:",im_height)
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵
im_proj = dataset.GetProjection() # 地图投影信息
im_data = dataset.ReadAsArray(0, 0, im_width, im_height) # 将数据写成数组,对应栅格矩阵
del dataset
return im_proj, im_geotrans, im_data
def L2OneLayerDN2Reflectance(self,info_dict, band, untar_dir):
'''
L2数据单波段定标为反射率
@ info_dict : 定标信息
@ band : 定标波段
@ untar_dir : 解包文件夹
'''
filename = os.path.join(untar_dir, info_dict['FILE_NAME_BAND']['FILE_NAME_BAND_' + str(band)])
REFLECTANCE_MULT_BAND = info_dict['REFLECTANCE_MULT_BAND']['REFLECTANCE_MULT_BAND_' + str(band)]
REFLECTANCE_ADD_BAND = info_dict['REFLECTANCE_ADD_BAND']['REFLECTANCE_ADD_BAND_' + str(band)]
# print(filename)
im_proj, im_geotrans, im_data = self.read_img(filename)
im_data = np.where(im_data, im_data, np.nan)
im_data = REFLECTANCE_MULT_BAND * im_data.astype(np.float32) + REFLECTANCE_ADD_BAND
return im_proj, im_geotrans, im_data
def L2DN2Reflectance(self, info_dict, untar_dir, bands=None):
'''
L2数据单景影像定标为反射率
@ info_dict : 定标信息
@ bands : 定标波段数
@ untar_dir : 解包文件夹
'''
if bands is None:
bands = [1, 2, 3, 4, 5, 6, 7]
data = []
for band in bands:
im_proj, im_geotrans, im_data = self.L2OneLayerDN2Reflectance(info_dict, band, untar_dir)
data.append(im_data)
return im_proj, im_geotrans, np.array(data)
def save_result(self,result_folder, tar_file, im_proj, im_geotrans, im_data):
'''
波段合成
'''
filename = os.path.join(result_folder, tar_file.split('.')[0] + '.TIF')
driver = gdal.GetDriverByName('GTiff')
(d, col, row) = im_data.shape
outraster = driver.Create(filename, row, col, d, gdal.GDT_Float32)
outraster.SetGeoTransform(im_geotrans)
outraster.SetProjection(im_proj)
for i in range(d):
outraster.GetRasterBand(i + 1).WriteArray(im_data[i])
outraster.FlushCache()
del outraster
def L2_batch_DN2Reflectance(self, work_place, bands=None):
'''
多景影像批量定标
'''
if bands is None:
bands = [1, 2, 3, 4, 5, 6, 7]
result_folder = self.create_filefolder(work_place) # 创建结果保存文件夹
work_file_dir = os.listdir(work_place)
for num,tar_file in enumerate(work_file_dir):
if '.tar' in tar_file:
untar_dir = self.untar(os.path.join(work_place, tar_file))
MTL_data = self.read_MTL(untar_dir)
info_dict = self.get_L2_MTL_info(MTL_data)
im_proj, im_geotrans, im_data = self.L2DN2Reflectance(info_dict, untar_dir, bands)
self.save_result(result_folder, tar_file, im_proj, im_geotrans, im_data)
self.remove_dir(untar_dir)
print("{}.{} have done!".format(num,tar_file))
def L1OneLayerDN2Reflectance(self,info_dict, band, untar_dir):
'''
L1单波段影像定标为反射率
@ info_dict : 定标信息
@ band : 定标波段
@ untar_dir : 解包文件夹
'''
filename = os.path.join(untar_dir, info_dict['FILE_NAME_BAND']['FILE_NAME_BAND_' + str(band)])
REFLECTANCE_MULT_BAND = info_dict['REFLECTANCE_MULT_BAND']['REFLECTANCE_MULT_BAND_' + str(band)]
REFLECTANCE_ADD_BAND = info_dict['REFLECTANCE_ADD_BAND']['REFLECTANCE_ADD_BAND_' + str(band)]
SUN_ELEVATION = info_dict['SUN_ELEVATION']*np.pi/180.0
im_proj, im_geotrans, im_data = self.read_img(filename)
im_data = np.where(im_data, im_data, np.nan)
im_data = (REFLECTANCE_MULT_BAND * im_data.astype(np.float32) + REFLECTANCE_ADD_BAND)/np.sin(SUN_ELEVATION)
return im_proj, im_geotrans, im_data
def L1DN2Reflectance(self, info_dict, untar_dir, bands=None):
'''
单景影像定标为反射率
@ info_dict : 定标信息
@ bands : 定标波段数
@ untar_dir : 解包文件夹
'''
if bands is None:
bands = [1, 2, 3, 4, 5, 6, 7, 9]
data = []
for band in bands:
im_proj, im_geotrans, im_data = self.L1OneLayerDN2Reflectance(info_dict, band, untar_dir)
data.append(im_data)
return im_proj, im_geotrans, np.array(data)
def L1_batch_DN2Reflectance(self, work_place, bands=None):
'''
多景影像批量定标
'''
if bands is None:
bands = [1, 2, 3, 4, 5, 6, 7, 9]
result_folder = self.create_filefolder(work_place) # 创建结果保存文件夹
work_file_dir = os.listdir(work_place)
for num,tar_file in enumerate(work_file_dir):
if '.tar' in tar_file:
untar_dir = self.untar(os.path.join(work_place, tar_file))
MTL_data = self.read_MTL(untar_dir)
info_dict = self.get_L1_MTL_info(MTL_data)
im_proj, im_geotrans, im_data = self.L1DN2Reflectance(info_dict, untar_dir, bands)
self.save_result(result_folder, tar_file, im_proj, im_geotrans, im_data)
self.remove_dir(untar_dir)
print("{}.{} have done!".format(num,tar_file))
def L1OneLayerDN2Radiance(self,info_dict, band, untar_dir):
'''
单波段定标为辐亮度
@ info_dict : 定标信息
@ band : 定标波段
@ untar_dir : 解包文件夹
'''
filename = os.path.join(untar_dir, info_dict['FILE_NAME_BAND']['FILE_NAME_BAND_' + str(band)])
REFLECTANCE_MULT_BAND = info_dict['RADIANCE_MULT_BAND']['RADIANCE_MULT_BAND_' + str(band)]
REFLECTANCE_ADD_BAND = info_dict['RADIANCE_ADD_BAND']['RADIANCE_ADD_BAND_' + str(band)]
im_proj, im_geotrans, im_data = self.read_img(filename)
im_data = np.where(im_data, im_data, np.nan)
im_data = (REFLECTANCE_MULT_BAND * im_data.astype(np.float32) + REFLECTANCE_ADD_BAND)
return im_proj, im_geotrans, im_data
def L1DN2Radiance(self, info_dict, untar_dir, bands=None):
'''
单景定标为辐亮度
@ info_dict : 定标信息
@ bands : 定标波段数
@ untar_dir : 解包文件夹
'''
if bands is None:
bands = [1, 2, 3, 4, 5, 6, 7, 9]
data = []
for band in bands:
im_proj, im_geotrans, im_data = self.L1OneLayerDN2Radiance(info_dict, band, untar_dir)
data.append(im_data)
return im_proj, im_geotrans, np.array(data)
def L1_batch_DN2Radiance(self, work_place, bands = None):
'''
批量定标为辐亮度
'''
# 1. 准备工作
if bands is None:
bands = [1, 2, 3, 4, 5, 6, 7, 9]
result_folder = self.create_filefolder(work_place) # 创建结果保存文件夹
work_file_dir = os.listdir(work_place)
for num,tar_file in enumerate(work_file_dir):
if '.tar' in tar_file:
untar_dir = self.untar(os.path.join(work_place, tar_file))
MTL_data = self.read_MTL(untar_dir)
info_dict = self.get_L1_MTL_info(MTL_data)
im_proj, im_geotrans, im_data = self.L1DN2Radiance(info_dict, untar_dir, bands)
self.save_result(result_folder, tar_file, im_proj, im_geotrans, im_data)
self.remove_dir(untar_dir)
print("{}.{} have done!".format(num,tar_file))
if __name__ == '__main__' :
RC = Radiometric_Calibration()
nocloud_dir = r''#此处填写需要定标的文件所在文件夹的路径
RC.L1_batch_DN2Reflectance(nocloud_dir,[1, 2, 3, 4, 5])
参考链接
[1] https://www.usgs.gov/landsat-missions/using-usgs-landsat-level-1-data-product