使用Python使用大气校正法计算地表温度
2024-04-25 05:06:20
使用Python使用大气校正法计算地表温度
前言
也有段时间没有跟新博客了，这次博客就是用新学的python语言来进行一个地表温度的计算，也算是承接了之前的内容吧！
一、具体原理及方法
这里不再赘述，有关原理及方法可以查阅我之前的博客文章https://blog.csdn.net/qq_44935981/article/details/90319487?spm=1001.2014.3001.5501
二、代码及步骤如下
1.引入库
代码如下：
import os
import numpy as np
from osgeo import gdal
import glob
import math
import matplotlib.pyplot as plt
#author sansuixueshi
2.读入数据及处理
代码如下：
input_path = 'E:/study/2018.4.09/2018.4.9/LC81300372018099LGN00/'
output_path = 'E:/study/2018.4.09/2018.4.9/LC81300372018099LGN00/out_data/'
file_all = os.listdir(input_path)
file_tail='.TIF'
radiance_scale=[1.2519E-02,1.2819E-02,1.1813E-02,9.9612E-03,6.0957E-03,1.5160E-03,5.1096E-04,1.1273E-02,2.3824E-03,3.3420E-04,3.3420E-04]
radiance_add=[-62.59278,-64.09577,-59.06370,-49.80584,-30.47869,-7.57977,-2.55479,-56.36651,-11.91176,0.10000,0.10000]
i=0if not os.path.exists(output_path):os.mkdir(output_path)for file_i in file_all:if file_i.endswith(file_tail):file_name = input_path + file_i(prename,suffix) = os.path.split(file_i)if i == 5:band = gdal.Open(file_name)red = band.ReadAsArray()*radiance_scale[3]+radiance_add[3]red = red.astype(np.float32)if i==6:band = gdal.Open(file_name)nir = band.ReadAsArray()*radiance_scale[4]+radiance_add[4]nir = nir.astype(np.float32)if i==1:band = gdal.Open(file_name)thermal = band.ReadAsArray()*radiance_scale[9]+radiance_add[9]thermal = thermal.astype(np.float32)#输出第十波段影像gtiff_driver=gdal.GetDriverByName('GTiff')out_file = gtiff_driver.Create(output_path + 'band_10.tif',thermal.shape[1],thermal.shape[0],1,6)out_file.SetProjection(out_file.GetProjection())out_file.SetGeoTransform(out_file.GetGeoTransform())out_band=out_file.GetRasterBand(1)out_band.WriteArray(thermal)out_band.FlushCache()i=i+1ndvi = (nir - red)/(nir + red)
nan_index=np.isnan(ndvi)
ndvi[nan_index]=0
#将输出影像转为geotiff格式
gtiff_driver=gdal.GetDriverByName('GTiff')
out_file = gtiff_driver.Create(output_path + 'ndvi.tif',ndvi.shape[1],ndvi.shape[0],1,6
)
out_file.SetProjection(out_file.GetProjection())
out_file.SetGeoTransform(out_file.GetGeoTransform())
#将影像展示，这一步可以省略
out_band=out_file.GetRasterBand(1)
out_band.WriteArray(ndvi)
out_band.FlushCache()
plt.imshow(ndvi)
plt.axis('off')
plt.show()#计算植被覆盖度
one_index = np.where(ndvi > 0.7)
zero_index = np.where(ndvi < 0.05)
ndvi[one_index] = 1.0
ndvi[zero_index] = 0.0
pv = ndvi#计算地表比辐射率
bfsl = pv*0.004 + 0.986#查询大气剖面参数，网站：http://atmcorr.gsfc.nasa.gov/
t=0.96
lu=0.26
ld=0.46
black_lightness=(thermal - lu - t*(1 - bfsl)*ld)/(t*bfsl)#计算地表温度
num = 774.89/black_lightness + 1
LST=(1321.08)/np.log(num)#减去273即为摄氏温度#输出LST影像
gtiff_driver = gdal.GetDriverByName('GTiff')
out_file = gtiff_driver.Create(output_path + 'LST.tif',LST.shape[1],LST.shape[0],1,6
)
out_file.SetProjection(out_file.GetProjection())
out_file.SetGeoTransform(out_file.GetGeoTransform())
out_band=out_file.GetRasterBand(1)
out_band.WriteArray(LST)
out_band.FlushCache()#author sansuixueshi
处理结果与全部代码
代码如下：
import os
import numpy as np
from osgeo import gdal
import glob
import math
import matplotlib.pyplot as plt#读取影像数据对应波段并进行定标操作
input_path = 'E:/study/2018.4.09/2018.4.9/LC81300372018099LGN00/'
output_path = 'E:/study/2018.4.09/2018.4.9/LC81300372018099LGN00/out_data/'
file_all = os.listdir(input_path)
file_tail='.TIF'
radiance_scale=[1.2519E-02,1.2819E-02,1.1813E-02,9.9612E-03,6.0957E-03,1.5160E-03,5.1096E-04,1.1273E-02,2.3824E-03,3.3420E-04,3.3420E-04]
radiance_add=[-62.59278,-64.09577,-59.06370,-49.80584,-30.47869,-7.57977,-2.55479,-56.36651,-11.91176,0.10000,0.10000]
i=0if not os.path.exists(output_path):os.mkdir(output_path)for file_i in file_all:if file_i.endswith(file_tail):file_name = input_path + file_i(prename,suffix) = os.path.split(file_i)if i == 5:band = gdal.Open(file_name)red = band.ReadAsArray()*radiance_scale[3]+radiance_add[3]red = red.astype(np.float32)if i==6:band = gdal.Open(file_name)nir = band.ReadAsArray()*radiance_scale[4]+radiance_add[4]nir = nir.astype(np.float32)if i==1:band = gdal.Open(file_name)thermal = band.ReadAsArray()*radiance_scale[9]+radiance_add[9]thermal = thermal.astype(np.float32)#输出第十波段影像gtiff_driver=gdal.GetDriverByName('GTiff')out_file = gtiff_driver.Create(output_path + 'band_10.tif',thermal.shape[1],thermal.shape[0],1,6)out_file.SetProjection(out_file.GetProjection())out_file.SetGeoTransform(out_file.GetGeoTransform())out_band=out_file.GetRasterBand(1)out_band.WriteArray(thermal)out_band.FlushCache()i=i+1ndvi = (nir - red)/(nir + red)
nan_index=np.isnan(ndvi)
ndvi[nan_index]=0
#将输出影像转为geotiff格式
gtiff_driver=gdal.GetDriverByName('GTiff')
out_file = gtiff_driver.Create(output_path + 'ndvi.tif',ndvi.shape[1],ndvi.shape[0],1,6
)
out_file.SetProjection(out_file.GetProjection())
out_file.SetGeoTransform(out_file.GetGeoTransform())
#将影像展示，这一步可以省略
out_band=out_file.GetRasterBand(1)
out_band.WriteArray(ndvi)
out_band.FlushCache()
plt.imshow(ndvi)
plt.axis('off')
plt.show()#计算植被覆盖度
one_index = np.where(ndvi > 0.7)
zero_index = np.where(ndvi < 0.05)
ndvi[one_index] = 1.0
ndvi[zero_index] = 0.0
pv = ndvi#计算地表比辐射率
bfsl = pv*0.004 + 0.986#查询大气剖面参数，网站：http://atmcorr.gsfc.nasa.gov/
t=0.96
lu=0.26
ld=0.46
black_lightness=(thermal - lu - t*(1 - bfsl)*ld)/(t*bfsl)#计算地表温度
num = 774.89/black_lightness + 1
LST=(1321.08)/np.log(num)#输出LST影像
gtiff_driver = gdal.GetDriverByName('GTiff')
out_file = gtiff_driver.Create(output_path + 'LST.tif',LST.shape[1],LST.shape[0],1,6
)
out_file.SetProjection(out_file.GetProjection())
out_file.SetGeoTransform(out_file.GetGeoTransform())
out_band=out_file.GetRasterBand(1)
out_band.WriteArray(LST)
out_band.FlushCache()#author sansuixueshi
结果图如上
使用Python使用大气校正法计算地表温度相关推荐
- ENVI下Landsat8大气校正法反演地表温度
尝试过ETM+热红外的大气校正法.单窗算法:ASTER的单窗算法.劈窗算法:Landsat8的大气校正法.单窗算法.劈窗算法反演地表温度.这里首先想记录Landsat8的大气校正法反演地表温度. 我认 ...
- 如何使用 Landsat 8 卫星影像计算地表温度
首先当然是要下载Landsat 8卫星影像.打开以下网址 https://eos.com/landviewer/或者https://earthexplorer.usgs.gov/或者 访问之前数据资源 ...
- GEEer成长日记十三：Landsat_SR计算地表温度时间序列
更多精彩内容请关注微信公众号:GEEer成长日记 上期我们介绍了Modis_LST产品MODIS/006/MOD11A1的时间序列,因为这款产品是官方已经经过各种矫正和处理的产品,精度较高,且范围广, ...
- Google Earth Engine（GEE）——python法国里昂地区的地表温度和地面高程的静态制图
获取静态地图 现在,我们想要获得感兴趣区域周围地表温度和地面高程的静态地图.我们使用法国里昂周围 1000 公里的缓冲区来定义这个感兴趣的区域. # 用里昂周围 1000 公里的缓冲区定义一个感兴趣的 ...
- Landsat系列卫星地表温度批量反演代码（大气校正法）
针对Landsat5.Landsat7.Landsat8的热红外波段反演地表温度的代码,可以批量进行温度的反演,但需要有前期的一些准备,包括大气校正参数的获取.可见光波段的大气校正等,以及文件夹的准备 ...
- [ENVI] 定量遥感实验-地表温度反演与地表温度测定 (超详细步骤)
实验目的 得到地表气温专题图反演结果之间的散点图 实验内容及实验步骤 1.MODIS地表温度产品的使用 软件环境:ENVI及MRT(已提供,请提前安装好)或MCTK扩展工具 实验数据:地表温度与发射率 ...
- （五）Landat_5 TM 遥感影像计算NDVI、MNDWI、NDBI以及地表温度反演
流程概况图 一.植被覆盖指数(NDVI) 所谓植被指数,就是利用多波段遥感图像的可见光波段以及近红外波段相组合,形成能够反映突出植被分布状况和强度的指数.植被指数的求解方法不唯一,定义的植被指数已经多 ...
- 地表温度反演操作总结
地表温度反演操作总结 由于博主也是自己边做实验便总结,难免会有混乱出错的地方,欢迎大家指正~ 整体实验流程: 1 数据预处理 1.1 数据下载 本实验使用的是Landsat8 OLI_TIRS数据,在 ...
- 基于ENVI与ERDAS的Landsat 7 ETM+单窗算法地表温度（LST）反演
基于ENVI与ERDAS的Landsat 7 ETM+单窗算法地表温度(LST)反演 1 原理部分与前期操作准备 1.1 图像预处理 1.2 植被指数反演 1.3 单窗算法原理 2 实际操作部分 2. ...
最新文章
- latex教程详细笔记
- 2008-2013年写的10个小软件
- 在Qt中使用C++代码创建界面
- sql：CallableStatement执行存储过程
- power bi 中计算_Power BI中的期间比较
- oracle+资料类型不一致吗,oracle数据库中，字段类型不一致，导致查询慢
- iOS 开发中，单款应用程序的最大可用内存是多少？
- 一个农民父亲令人震撼的力量
- vue样式控制的方式
- c语言字符合法,C语言字符数据的合法形式
- VS系列编译器基本调试快捷键的使用
- 超棒工具8个高质量图标搜索引擎
- Molecular Plant: 王二涛组及合作者揭示丛枝菌根共生与根瘤共生的协同进化机制...
- 青龙面板-喜马L雅签到
- UG NX 12 点集
- 六大类二叉树面试题汇总解答
- BIM建模助手上线一周，有哪些BUG被用户找到？
- 高德地图自定义点标记大小_高德地图 自定义点标记 图标大小
- 什么是微信不死域名？
- ANSYS中按照X坐标提取节点应力值