利用python进行Landsat8数据的辐射定标

最新推荐文章于 2024-01-07 18:05:42 发布

原创最新推荐文章于 2024-01-07 18:05:42 发布 · 5.3k 阅读

59 ·

CC 4.0 BY-SA版权

文章标签：

#python #opencv #RS #gdal #Landsat8

遥感图像处理专栏收录该内容

16 篇文章

订阅专栏

介绍了一种用于Landsat8遥感影像的辐射定标方法，该方法通过读取元文件中的定标参数，对可见光波段进行辐射定标处理，确保NaN数据的正确处理，并提供定标后数据的写入功能。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

上次用于处理Landsat8数据的类，增加了辐射定标的功能。目前是用于可见光波段的辐射定标。

这段辐射定标的代码，首先从元文件中读取定标参数，然后进行辐射定标。对NaN数据做了处理，使其一直保持为NaN。最后可以利用Landsat8类中的write函数，将辐射定标后的数据写出。

当然，这段代码的有效性，仅仅保证在我的数据上好用。

具体辐射定标的实现，参看Landsat8Reader类的radiometric_calibration函数。

import os
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 cv

class Landsat8Reader(object):
    def __init__(self):
        self.base_path =\
        "LC81220352018123LGN00/LC08_L1TP_122035_20180503_20180516_01_T1"
        self.bands = 7
        self.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).DataType
        image = 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.nan
        
        del ds    
        
        return image
    
    def write(self, image, file_name, bands, format = 'ENVI'):
        ds = gdal.Open(self.band_file_name[0])
        projection = ds.GetProjection()
        geotransform = ds.GetGeoTransform()  
        x_size = ds.RasterXSize
        y_size = ds.RasterYSize
        del ds
        band_count = bands
        dtype = gdal.GDT_Float32
        
        driver = 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_ds
    
    def 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')
            metadata = f.readlines()
            f.close()            
            multi_parameters = []
            add_parameters = []
            parameter_start_line = 0
            for lines in metadata:
                test_line = lines.split('=')
                if test_line[0] == '    RADIANCE_MULT_BAND_1 ':
                    break
                else:
                    parameter_start_line = parameter_start_line + 1

            for lines in range(parameter_start_line, parameter_start_line+11):
                parameter = float(metadata[lines].split('=')[1])
                multi_parameters.append(parameter)

            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_parameters
        
        multi_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]
            cali_image[:, :, band] = image[:, :, band]*gain + offset
            
        del image
        return cali_image

if __name__ == "__main__":
	reader = Landsat8Reader()
	image = reader.read()
	cali_image = reader.radiometric_calibration()
	file_path = './test/write_cali_image'
	reader.write(cali_image, file_path, reader.bands)

经过检验，与envi的定标结果完全相同。

额外用到的拉伸显示函数：

import numpy as np
import cv2
from matplotlib import pyplot as plt
from osgeo import gdal
from osgeo import gdal_array

def TwoPercentLinear(image, max_out=255, min_out=0):
    b, g, r = cv2.split(image)
    def gray_process(gray, maxout = max_out, minout = min_out):
        high_value = np.percentile(gray, 98)
        low_value = np.percentile(gray, 2)
        truncated_gray = np.clip(gray, a_min=low_value, a_max=high_value) 
        processed_gray = ((truncated_gray - low_value)/(high_value - low_value)) * (maxout - minout)
        return processed_gray
    r_p = gray_process(r)
    g_p = gray_process(g)
    b_p = gray_process(b)
    result = cv2.merge((b_p, g_p, r_p))
    return np.uint8(result)