python对landsat8数据进行辐射校正
Landsat8 L1 T數據是輻射校正數據使用地面控制點和數字高程模型數據進行精確校正后的數據產品,還需要做輻射校正(輻射定標和大氣校正)。
一、python
1.輻射定標
輻射亮度L=DN*Gain+Bias
(1)代碼
傾向于TIF格式的影像,所以output format選擇的是“GTIFF”格式。
from osgeo import gdal from osgeo import gdal_array import numpy as np from show import TwoPercentLinear from matplotlib import pyplot as plt import cv2 as cvclass Landsat8Reader(object):def __init__(self):self.base_path = r"D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\LC08_L1TP_134036_20170808_20170813_01_T1"self.bands = 7self.band_file_name = []self.nan_position = []def read(self):for band in range(self.bands):band_name = self.base_path + "_B" + str(band + 1) + ".tif"self.band_file_name.append(band_name)ds = gdal.Open(self.band_file_name[0])image_dt = ds.GetRasterBand(1).DataTypeimage = np.zeros((ds.RasterYSize, ds.RasterXSize, self.bands),dtype=np.float)for band in range(self.bands):ds = gdal.Open(self.band_file_name[band])band_image = ds.GetRasterBand(1)image[:, :, band] = band_image.ReadAsArray()self.nan_position = np.where(image == 0)image[self.nan_position] = np.nandel dsreturn imagedef write(self, image, file_name, bands, format='GTIFF'):ds = gdal.Open(self.band_file_name[0])projection = ds.GetProjection()geotransform = ds.GetGeoTransform()x_size = ds.RasterXSizey_size = ds.RasterYSizedel dsband_count = bandsdtype = gdal.GDT_Float32driver = gdal.GetDriverByName(format)new_ds = driver.Create(file_name, x_size, y_size, band_count, dtype)new_ds.SetGeoTransform(geotransform)new_ds.SetProjection(projection)for band in range(self.bands):new_ds.GetRasterBand(band + 1).WriteArray(image[:, :, band])new_ds.FlushCache()del new_dsdef show_true_color(self, image):index = np.array([3, 2, 1])true_color_image = image[:, :, index].astype(np.float32)strech_image = TwoPercentLinear(true_color_image)plt.imshow(strech_image)def show_CIR_color(self, image):index = np.array([4, 3, 2])true_color_image = image[:, :, index].astype(np.float32)strech_image = TwoPercentLinear(true_color_image)plt.imshow(strech_image)def radiometric_calibration(self):# 輻射定標image = self.read()def get_calibration_parameters():filename = self.base_path + "_MTL" + ".txt"f = open(filename, 'r')# f = open(r'D:\ProfessionalProfile\LandsatImage\LC08_L1TP_134036_20170808_20170813_01_T1\LC08_L1TP_134036_20170808_20170813_01_T1_MTL.txt', 'r')# readlines()方法讀取整個文件所有行,保存在一個列表(list)變量中,每行作為一個元素,但讀取大文件會比較占內存。metadata = f.readlines()f.close()multi_parameters = []add_parameters = []parameter_start_line = 0for lines in metadata:# split()方法通過指定分隔符對字符串進行切片,如果參數 num 有指定值,則分隔 num+1 個子字符串。# 無指定值,默認為 -1, 即分隔所有符合要求的。test_line = lines.split('=')if test_line[0] == ' RADIANCE_MULT_BAND_1 ':breakelse:parameter_start_line = parameter_start_line + 1for lines in range(parameter_start_line, parameter_start_line + 11):parameter = float(metadata[lines].split('=')[1])multi_parameters.append(parameter)# 由于RADIANCE_MULT_BAND參數和RADIANCE_ADD_BAND參數是挨著的,所以直接+11個參數即可for lines in range(parameter_start_line + 11, parameter_start_line + 22):parameter = float(metadata[lines].split('=')[1])add_parameters.append(parameter)return multi_parameters, add_parametersmulti_parameters, add_parameters = get_calibration_parameters()cali_image = np.zeros_like(image, dtype=np.float32)for band in range(self.bands):gain = multi_parameters[band]offset = add_parameters[band]# 輻射定標像元 = DN * 增益 + 偏置,線性關系。cali_image[:, :, band] = image[:, :, band] * gain + offsetdel imagereturn cali_imageif __name__ == "__main__":reader = Landsat8Reader()image = reader.read()cali_image = reader.radiometric_calibration()file_path = r'D:\ProfessionalProfile\LandsatImage\1_RadiationCalibration0414\LC08.tif'reader.write(cali_image, file_path, reader.bands)(2)結果
兩幅圖基本沒什么變化,就是數值差了10倍。用ENVI軟件進行大氣校正的時候:點擊FLAASH Settings按鈕時,自動獲取輻射亮度單位轉換系數Scale Factor:0.1。估計這是差異的來源。
參考:https://blog.csdn.net/weixin_34472954/article/details/112359224
2.大氣校正
表觀反射率ρ=π*L*D2/(ESUN*cosθ)式中,ρ為大氣層頂( TOA) 表觀反射率( 無量綱),π為常量( 球面度sr),L?為大氣層頂進人衛星傳感器的光譜輻射亮度( W*1/(m2*sr*um)),D為日地之間距離( 天文單位),?ESUN為大氣層頂的平均太陽光譜輻照度(?W*1/(m2*sr*um))θ為太陽的天頂角。
見博文:https://blog.csdn.net/weixin_40501429/article/details/115856301
二、ENVI
1.輻射定標
?在 Radiometric Calibration 面板中,設置以下參數:
定標類型(Calibration Type) :輻射率數據 Radiance;Output Interleave:BIL;數據類型(Data Type) :Float;
?單擊 Apply FLAASH Settings 按鈕。
結果:
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?真彩色原圖&輻射定標圖
2.大氣校正
?
綠框框出來的找MTL設;。多光譜參數:設置?K-T 反演選擇默認模式:Defaults->Over-Land Retrieval standard(600:2100) ,其他默認。
結果
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 真彩色大氣校正圖
?
?
參考
1.博客
踏著路唱著歌博主:利用python進行Landsat8數據的輻射定標?https://blog.csdn.net/RSstudent/article/details/107521975
G1SLu博主:https://blog.csdn.net/ahlxt123/article/details/39568733?locationNum=9&fps=1
RS&Hydrology博主:https://blog.csdn.net/qq_32649321/article/details/111352217?utm_term=landsat%E6%95%B0%E6%8D%AE%E9%A2%84%E5%A4%84%E7%90%86&utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2~all~sobaiduweb~default-0-111352217&spm=3001.4430
Sam Qi博主:https://blog.csdn.net/weixin_34472954/article/details/112359224
2.論文
[1]初慶偉,張洪群,吳業煒,馮鐘葵,陳勃.Landsat-8衛星數據應用探討[J].遙感信息,2013,28(04):110-114.
[2]池宏康,周廣勝,許振柱,肖春旺,袁文平.表觀反射率及其在植被遙感中的應用[J].植物生態學報,2005(01):74-80.
總結
以上是生活随笔為你收集整理的python对landsat8数据进行辐射校正的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: [react] react怎么拿到组件对
- 下一篇: Python: 反方向迭代一个序列