Landsat系列卫星地表温度批量反演代码（大气校正法）

针对Landsat5、Landsat7、Landsat8的热红外波段反演地表温度的代码，可以批量进行温度的反演，但需要有前期的一些准备，包括大气校正参数的获取、可见光波段的大气校正等，以及文件夹的准备、文件的规范命名等。

1、大气校正法的主要过程参考文档：ENVI下Landsat8大气校正法反演地表温度
2、可见光波段参照以上文档进行大气校正（但是根据覃志豪老师的论文，可见光波段不进行大气校正也可直接用于温度反演，因此在文件命名中有所区别）
3、文件夹准备及文件命名：

（1）近红外波段文件夹
注意如果是Landsat8的数据，开头必须是"LC08"；Landsat7对应"LE07"；Landsat5对应"LT05"
如果近红外波段通过ENVI做了大气校正，则最后的命名必须是"flaash"，如果没有进行大气校正，则结尾不能是"flaash"
注意：如果不做大气校正，则直接把原始图像复制过来即可，不需要用ENVI做辐射定标，代码中会进行辐射定标

（2）红波段文件夹
注意如果是Landsat8的数据，开头必须是"LC08"；Landsat7对应"LE07"；Landsat5对应"LT05"
如果红波段通过ENVI做了大气校正，则最后的命名必须是"flaash"，如果没有进行大气校正，则结尾不能是"flaash"
注意：如果不做大气校正，则直接把原始图像复制过来即可，不需要用ENVI做辐射定标，代码中会进行辐射定标

（3）热红外波段文件夹
从USGS下载的Landsat系列数据，解压后，将热红外波段直接复制过来，不需要重命名

（4）热红外波段大气校正参数的EXCEL文件
注意TIR这一列，需要跟热红外文件夹中热红外波段的命名一致（不需要尾缀）

4、具体代码

import numpy as np
from osgeo import gdal
import matplotlib.ticker as tic
import matplotlib.pyplot as plt
import imageio
import pandas as pd
import os
import copy
import datetime
# 程序开始运行时的时间
start_time = datetime.datetime.now()
# gdal库读取遥感影像数据
def image_open(img):data = gdal.Open(img)if data == 'None':print("图像无法读取")return data
FilePath = r"D:/file/novel coronavirus/Infrared remote sensing/python_temperature"
SavePath = FilePath + '/' + 'Temperature'
if not os.path.exists(SavePath):os.makedirs(SavePath)
# 将三个文件夹中的文件名称存入内存中，并按字符串顺序进行排序
Red_image = sorted(os.listdir(FilePath + "/" + "RED"))  # 红波段影像位置
NIR_image = sorted(os.listdir(FilePath + "/" + "NIR"))  # 近红波段影像位置
TIR_image = sorted(os.listdir(FilePath + "/" + "TIR"))  # 热红外波段影像位置
# 读取大气校正文件
# 将列的名称重命名
# 并用TIR这一列作为索引列
# excel中，TIR列为热红外波段的文件名（不含有.tif）
ATM = pd.read_excel("D:/file/novel coronavirus/Infrared remote sensing/python_temperature/ATM.xlsx",names=['TIR', '大气上行辐射', '大气下行辐射', '大气透过率'], index_col='TIR')
ATM.head()
# 用来验证大气校正参数的读取
# a = ATM.at['LC08_L1TP_123039_20131120_20170428_01_T1_B10', '大气上行辐射']
# print(a)
# 判断语句来验证文件尾是不是tif文件,尽量保证每个文件夹里面除了tif文件之外没有其它文件，保证三个文件夹数据量一致
for i in range(len(TIR_image)):# 将红波段的影像文件名称存入第一列Red_for = os.path.splitext(Red_image[i])[0]# 文件尾缀存入第二列Red_end = os.path.splitext(Red_image[i])[1]NIR_for = os.path.splitext(NIR_image[i])[0]NIR_end = os.path.splitext(NIR_image[i])[1]TIR_for = os.path.splitext(TIR_image[i])[0]TIR_end = os.path.splitext(TIR_image[i])[1]# 获取Landsat数据名landsat_name = TIR_for.split('_')[0]# print(landsat_name)# 获取影像拍摄时间点，用来作为文件保存的命名savename = "/" + landsat_name + "_" + TIR_for.split('_')[3]# print(savename)if (Red_end.lower() == ".tif" and NIR_end.lower() == ".tif" and TIR_end.lower() == ".tif"):  # 注意影像的尾缀是大小写的问题RedImgPath = FilePath + "/" + "RED" + "/" + Red_image[i]NIRImgPath = FilePath + "/" + "NIR" + "/" + NIR_image[i]TIRImgPath = FilePath + "/" + "TIR" + "/" + TIR_image[i]# 读取数据# 先调用image_open函数打开影像# 然后用GetRasterBand()来获取波段# 然后用ReadAsArray()读取像素值，用此函数之前必须先获取波段# 可用ReadAsArray()函数设置读取像素值的数量及位置，将读取数据存放入矩阵中Reddata = image_open(RedImgPath)Redimg = Reddata.GetRasterBand(1).ReadAsArray().astype(np.float32)NIRdata = image_open(NIRImgPath)NIRimg = NIRdata.GetRasterBand(1).ReadAsArray().astype(np.float32)TIRdata = image_open(TIRImgPath)TIRimg = TIRdata.GetRasterBand(1).ReadAsArray().astype(np.float32)# 热红外波段辐射校正# landsat8有两个热红外波段，一般用第10波段，两个波段定标参数一样if landsat_name == 'LC08':# Landsat8红波段辐射定标参数gainRed = 0.010233offsetRed = -51.16713# Landsat8近红波段辐射定标参数gainNIR = 0.0062623offsetNIR = -31.31173# Landsat8热红外波段辐射定标参数gainTIR = 0.0003342offsetTIR = 0.1landsat8_end = TIR_for.split('_')[7]# Landsat8的两个热红外波段有效波长不同if landsat8_end == 'B10':# Landsat8第10波段有效波长lda = 10.904if landsat8_end == 'B11':# Landsat8第11波段有效波长lda = 12.003# Landsat7分为高增益与低增益波段，'1'为低增益波段，用于高温环境；'2'为高增益波段，用于常规条件if landsat_name == 'LE07':# Landsat7红波段辐射定标参数gainRed = 0.62165offsetRed = -5.62165# Landsat7近红波段辐射定标参数gainNIR = 0.63976offsetNIR = -5.73976# Landsat7热红外波段有效波长lda = 11.269# Landsat7热红外波段定标参数# 波段号保存在文件名称的第9个参数中landsat7_end = TIR_for.split('_')[-1]# Landsat7低增益波段辐射定标参数# 低增益波段适用于沙漠等温度较高的环境if landsat7_end == '1':gainTIR = 0.067087offsetTIR = -0.06709# Landsat7高增益波段辐射定标参数# 高增益波段适用于常规环境if landsat7_end == '2':gainTIR = 0.037205offsetTIR = 3.1628# Landsat5只有一个热红外波段if landsat_name == 'LT05':# Landsat5热红外波段有效波长lda = 11.457# Landsat5红波段辐射定标参数gainRed = 1.044offsetRed = -2.21398# Landsat5近红波段辐射定标参数gainNIR = 0.87602offsetNIR = -2.38602# Landsat5热红外波段辐射定标参数gainTIR = 0.055376offsetTIR = 1.18# print(gain)# print(offset)# 红波段辐射定标# 判断红波段是不是经过了ENVI的大气校正，因为红波段与近红外波段的处理过程应该一致，所以只用判断红波段即可red1 = Red_for.split('_')[-1]# 如果红波段与近红波段已经经过了ENVI的大气校正，则直接使用大气校正的结果作为处理的输入值，如果没有进行大气校正，则需要先用公式做辐射定标# 需要注意的是，ENVI大气校正处理结果的命名方式，应该以flaash结尾if ('flaash' in red1):Red_radiance = RedimgNIR_radiance = NIRimgelif not('flaash' in red1):Red_radiance = gainRed * Redimg + offsetRed# 近红波段辐射定标NIR_radiance = gainNIR * NIRimg + offsetNIR# 热红外波段辐射定标TIR_radiance = gainTIR * TIRimg + offsetTIR# 计算NDVI值# 忽略分母为0的情况np.seterr(divide='ignore', invalid='ignore')NDVI = (NIR_radiance - Red_radiance) / (NIR_radiance + Red_radiance)NDVI = NDVI.astype(np.float32)# 计算FVC值# 利用NDVI值来判断FVC值的计算方式# NDVI值大于0.7时认为是纯植被像元# NDVI值小于0.05时判断为裸土像元# NDVI值在0.05与0.7之间时判断为城市地表像元# 需要用深拷贝，不然会改变原来的数组内容FVC = copy.deepcopy(NDVI)FVC = (FVC - 0.05) / (0.7 - 0.05)FVC[NDVI > 0.7] = 1FVC[NDVI < 0.05] = 0# 计算地表比辐射率emissivity = copy.deepcopy(NDVI)# 注意在等号后面的式子中也加上条件，不然会出现维度问题# 同样用NDVI值进行判断# NDVI值大于等于0.7时为纯植被# NDVI值小于等于0时为水体# NDVI值在0到0.7之间时为城市地表emissivity[NDVI >= 0.7] = 0.9625 + 0.0614 * FVC[NDVI >= 0.7] - 0.0461 * (FVC[NDVI >= 0.7]) ** 2emissivity[NDVI < 0.7] = 0.9589 + 0.086 * FVC[NDVI < 0.7] - 0.0671 * (FVC[NDVI < 0.7]) ** 2emissivity[NDVI <= 0] = 0.995# 温度反演# 先做大气校正，从ATM表中提取出相应的参数,计算得到地表亮度温度up = ATM.at[TIR_for, '大气上行辐射']down = ATM.at[TIR_for, '大气下行辐射']trans = ATM.at[TIR_for, '大气透过率']# 计算出地表亮度温度BT = (TIR_radiance - up - trans * (1 - emissivity) * down) / (emissivity * trans)# 求出地表温度，并将单位转化为摄氏度Temperature = 14387.7 / (lda * np.log(119104000 / (BT * (lda ** 5)) + 1)) - 273.15Temperature = Temperature.astype(np.float32)# 输出数据gtiff_driver = gdal.GetDriverByName("GTiff")out_ds = gtiff_driver.Create(SavePath + savename + "T.tif",Temperature.shape[1], Temperature.shape[0], 1, gdal.GDT_Float32)# 设置输出数据坐标投影为原始影像的坐标投影out_ds.SetProjection(TIRdata.GetProjection())out_ds.SetGeoTransform(TIRdata.GetGeoTransform())out_band = out_ds.GetRasterBand(1)out_band.WriteArray(Temperature)out_band.FlushCache()# 输出比辐射率# emissivity_driver = gdal.GetDriverByName("GTiff")# out_e = emissivity_driver.Create(SavePath + savename + "E.tif",#                              emissivity.shape[1], emissivity.shape[0], 1, gdal.GDT_Float32)# out_e.SetProjection(TIRdata.GetProjection())# out_e.SetGeoTransform(TIRdata.GetGeoTransform())# out_band_e = out_e.GetRasterBand(1)# out_band_e.WriteArray(emissivity)# out_band_e.FlushCache()# 输出辐射亮度图# BT_driver = gdal.GetDriverByName("GTiff")# out_BT = BT_driver.Create(SavePath + savename + "BT.tif",#                              BT.shape[1], BT.shape[0], 1, gdal.GDT_Float32)# out_BT.SetProjection(TIRdata.GetProjection())# out_BT.SetGeoTransform(TIRdata.GetGeoTransform())# out_band_BT = out_BT.GetRasterBand(1)# out_band_BT.WriteArray(BT)# out_band_BT.FlushCache()# 处理完一张后显示print('Temperature Retrieval: {}, {}'.format(i+1, landsat_name + "_" + TIR_for.split('_')[3]))
# 程序结束运行时的时间
end_time = datetime.datetime.now()
# # 计算温度反演所用的时间运行的时间
# print('Temperature Retrieval Time: {}'.format(end_time - start_time))

Landsat系列卫星地表温度批量反演代码（大气校正法）相关推荐

Landsat系列卫星介绍及影像下载
1.LandSat系列卫星介绍: 1.Landsat系列卫星概述: 美国NASA的陆地卫星(Landsat)计划从1972年7月23日以来,已发射8颗(第6颗发射失败).目前Landsat1-4均相继 ...
Landsat系列卫星WRS条带号Path Row分布介绍与对照图
WRS,即Worldwide Reference System,是Landsat系列卫星全球影像标记符号系统,用以区分全球各区域对应的Landsat系列卫星影像编号:其用"Path&qu ...
使用大气校正法对landsat-8tirs地表温度进行反演
软件平台 ENVI5.3软件具体方法及原理原理:首先估计大气对地表热辐射的影响, 然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去, 从而得到地表热辐射强度, 再把这一热辐射强度转化为相 ...
Landsat系列卫星：Landsat 9 详解和细节（NASA/USGS）
Landsat 9 是美国国家航空航天局(NASA) 和美国地质调查局 (USGS)之间的合作伙伴,将继续 Landsat 计划在重复全球观测以监测.了解和管理地球自然资源方面的关键作用. 自 197 ...
Landsat系列卫星全球参考系统，指定的PATH和ROW编号详细介绍
全球参考系统(WRS)是陆地卫星数据的一个全球符号系统.它使用户能够通过指定一个由PATH和ROW号码指定的名义场景中心来查询世界上任何部分的卫星图像.事实证明,WRS对于编目.参考和日常使用从Lan ...
基于Google Earth Engine的Landsat单窗算法地表温度（LST）反演
基于Google Earth Engine的Landsat单窗算法地表温度(LST)反演 1 背景知识 2 算法介绍 3 代码 4 效果 1 背景知识基于遥感数据的地表温度(LST)反演目前得到 ...
基于ENVI与ERDAS的Landsat 7 ETM+单窗算法地表温度（LST）反演
基于ENVI与ERDAS的Landsat 7 ETM+单窗算法地表温度(LST)反演 1 原理部分与前期操作准备 1.1 图像预处理 1.2 植被指数反演 1.3 单窗算法原理 2 实际操作部分 2. ...
ENVI5.3.1Landsat 8影像基于单窗算法和辐射传输方程进行地表温度反演
ENVI5.3.1基于Landsat 8影像进行辐射定标和大气校正文章目录一.为什么要进行辐射定标和大气校正? 二.详细步骤 1. 数据获取 2.数据预处理 2.1 辐射定标 2.1.1 多光谱波 ...
【转载】基于ENVI bandmath的地表温度反演
地表温度作为地球环境分析的重要指标,而遥感技术作为现代重要的对地观测手段,使得基于遥感图像的地表温度反演的研究越来越多.主要的地表温度反演方法有:大气校正法,单窗算法,单通道法等等.本文介绍用辐射传输 ...

Landsat系列卫星地表温度批量反演代码（大气校正法）

Landsat系列卫星地表温度批量反演代码（大气校正法）相关推荐

最新文章

热门文章