使用Python使用大气校正法计算地表温度

使用Python使用大气校正法计算地表温度

前言

也有段时间没有跟新博客了，这次博客就是用新学的python语言来进行一个地表温度的计算，也算是承接了之前的内容吧！

---

一、具体原理及方法

这里不再赘述，有关原理及方法可以查阅我之前的博客文章https://blog.youkuaiyun.com/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=0

if 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+1

ndvi = (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=0

if 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+1

ndvi = (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

---

结果图如上