利用GDAL(python)读取Landsat8数据

本人编写了一个利用GDAL读取Landsat8数据的代码，现在已经拥有了读取、真彩色、假彩色显示等功能。准备日后再进行扩展。

注意：路径是相对路径，使用需自行修改。
这个代码是将压缩包解压出来的数据的可见光、近红外和短波红外部分一齐读进一个三维数组里面。

下面是代码部分：

首先是主程序：

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 = []
    
    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 =\
                         gdal_array.GDALTypeCodeToNumericTypeCode(image_dt))
        for band in range(self.bands):
            ds = gdal.Open(self.band_file_name[band])
            band_image = ds.GetRasterBand(1)
            image[:,:, band] = band_image.ReadAsArray()
        
        return image
    
    def show_true_color(self, image):
        index = np.array([3, 2, 1])
        true_color_image = image[:, : , index].astype(np.float64)
        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.float64)
        strech_image = TwoPercentLinear(true_color_image)
        plt.imshow(strech_image)

if __name__ == "__main__":
    reader = Landsat8Reader()
    image = reader.read()
    reader.show_true_color(image)

还有一个函数TwoPercentLinear，命名规则不大对，是我早期编写的。这个代码是用于拉伸显示的。

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)