Landsat影像C2（L1/L2）数据批量辐射定标及波段合成

摘要： 本文主要介绍如何使用Python对Landsat8-9影像在本地进行批量预处理（辐射定标及波段合成）。是笔者经过对论文、数据手册和博客等一系列学习资料的一个总结。

---

前言

Landsat系列影像作为遥感领域使用最为广泛、时间序列最长的和影像质量最好的免费遥感数据之一。在科学实践的过程中，往往需要同时使用多张Landsat影像进行科学分析，对于原始数据的预处理过程（辐射定标和波段合成）成为定量分析前必不可少的步骤。
传统本地预处理方法一般采用ENVI或IDL的方法对数据进行预处理，但是：（1）老版本的ENVI对于卫星最新处理结果（Landsat Collection2）的兼容性较差（需要动手调整头文件），且面对众多需要预处理数据时，操作较多耗时较长；（2）IDL是针对遥感数据处理的专用语言，有一定的学习成本。

以下将介绍如何使用Python对Landsat影像在本地进行批量辐射定标和波段合成的预处理工作。

---

一、概念辨析

一般来讲，辐射定标是指，将遥感图像的数字量化值转化为辐射亮度、反射率或表面温度等物理量的处理过程。其中：
数字量化值：又称DN值，是遥感影像数据的原始值，数值类型一般为整数。（Landsat8-9有效范围为：1-65455）
辐射亮度：Radiance，是某一个面积辐射能量的总和。（有单位，如：$\mu W/(cm2\ast sr\ast 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_{\lambda }=M_L{Q}_{cal}+A_L$
其中：
$L_{\lambda }$ = TOA辐射亮度
$M_L$ =MLT文件中的乘法重缩放因子（RADIANCE_MULT_BAND_x，其中x是波段号）
$A_L$ =MLT文件中的加法重缩放因子 （RADIANCE_ADD_BAND_x, 其中x是波段号）
$Q_{cal}$ =量化和校准的标准产品像素值（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反射率：
${\rho }_{\lambda }=(M_{\rho }\cdot Q_{cal}+A_{\rho })/\sin ({\theta }_{SE})$
其中：
${\rho }_{\lambda }$=TOA反射率
$M_{\rho }$=MLT文件中的乘法重缩放因子（REFLECTANCE_MULT_BAND_x，其中x是波段号）
$A_{\rho }$=MLT文件中的乘法重缩放因子（REFLECTANCE_ADD_BAND_x，其中x是波段号）
$Q_{cal}$=量化和校准的标准产品像素值（DN值）
${\theta }_{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_{\lambda }=M_L{Q}_{cal}+A_L$
其中：
$L_{\lambda }$ = TOA辐射亮度
$M_L$ =MLT文件中的乘法重缩放因子（RADIANCE_MULT_BAND_x，其中x是波段号）
$A_L$ =MLT文件中的加法重缩放因子 （RADIANCE_ADD_BAND_x, 其中x是波段号）
$Q_{cal}$ =量化和校准的标准产品像素值（DN值）
注：在L2数据中$M_L$和$A_L$为常数分别是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