编程 | Landsat9影像转真彩色RGBA缩略图

Landsat 9

从USGS下载Landsat9影像，Landsat9影像和国产影像不一样的地方在于：

1.Landsat9影像已做好了几何校正，无rpc文件。

2.Landsat9影像以单波段的形式，将每个波段保存为单独的文件。

Landsat9 的真彩色RGBA

众所周知，Landsat9的第2波段是蓝色，Blue,简写B

Landsat9的第3波段是绿色，Green,简写G

Landsat9的第4波段是红色，Red,简写R

接下来，要写一个函数，实现这个功能：

输入Landsat9解压后的文件夹，得到Landsat9的第2、3、4波段的文件名

这个函数的作用是，有利于自动化，每景影像的名字是唯一的，如果每次都要自己手动指定Landsat9的第2、3、4波段的文件名是费时费力的。

注意分析Landsat9的tif名字的格式，提炼出共同的命名规则。

image-20231013172543805

LC09_L2SR_018113_20230209_20230310_02_T2_SR_B1.TIF 的特点是

B1代表band 1 ，第一波段。

该文件的后缀是大写字母 TIF。

根据这两点可以自动获取第2、3、4波段的名字。

细看以下代码：

class landsat9_to_webp:
    def __init__(self, path):
       self.path = path
    def get_tif_file(self, file_dir, type = '.TIF'):
        # 获取TIF文件
        """
        搜索 后缀名为type的文件  不包括子目录的文件
        """
        L = []
        if type(type) == str:
                if len([type]) == 1:
                    filelist = os.listdir(file_dir)
                    for file in filelist:
                        if os.path.splitext(file)[1] == type:
                            L.append(os.path.join(file_dir, file))
        if type(type) != str:
                if len(type) > 1:
                    for i in range(len(type)):
                        filelist = os.listdir(file_dir)
                        for file in filelist:
                            if os.path.splitext(file)[1] == type[i]:
                                L.append(os.path.join(file_dir, file))
        return L

我们先写一个类，这个类叫做landsat9_to_webp，我们常常用self代指这个类。landsat9_to_webp这个类的输入是lansat9文件夹路径，所以这个类的其中一个属性是输入路径。

这个类 有一个类方法，叫做 get_tif_file，这个也是landsat9_to_webp是属性。get_tif_file方法实现查找 指定文件夹下 后缀为TIF的文件，若文件存在则返回符合条件的文件名字的列表。

下面是 示例

path = r'X:\LC09_L2SP_121038_20220616_20230411_02_T1'
filelist = landsat9_to_webp(path).get_tif_file()

我们把filelist循环打印出来，瞧瞧它到底长啥样。

import os
for file in filelist:
    print(os.path.basename(file))

控制台打印一下信息：

LC09_L2SP_121038_20220616_20230411_02_T1_QA_PIXEL.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_QA_RADSAT.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B1.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B2.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B3.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B4.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B5.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B6.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B7.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_QA_AEROSOL.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_ATRAN.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_B10.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_CDIST.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_DRAD.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_EMIS.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_EMSD.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_QA.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_TRAD.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_ST_URAD.TIF

显然，filelist是一个长度为19的列表，下图是在pycharm的截图

image-20231013180036859

这时候，我们已经成功一半了。

下一步，是对这个列表进行数据清洗，目的是稳定地找到我们所需要的波段。

有基础的朋友会注意到，filelist的第4个元素、第5个元素、第6个元素就是我们所需要的B2、B3、B4波段。为什么还有做数据清理呢？ 理由有以下：

不能肯定、不能确定每次获得的filelist的第4个元素、第5个元素、第6个元素就是我们所需要的B2、B3、B4波段。

比如，我用sentinel-2举例，sentinel-2 L1C数据和sentinel-2 L2A数据就存在文件结构差异，所以不能保证每次都是filelist[3]、filelist[4]、filelist[5] 就是我们所要的文件名字。

这里我们再写几行代码，确保下一步读取是文件，就是Landsat9的第2、3、4波段。

笨方法就是循环filelist列表，加上一个if语句，如果当前元素存在B2.TIF，则该元素就是 第2波段的名字。

我们再写一个 类方法，

再次，整合代码，如下：

class landsat9_to_webp:
    def __init__(self, path):
       self.path = path
    def get_tif_file(self, filetype = '.TIF'):
        # 获取TIF文件
        """
        搜索 后缀名为type的文件  不包括子目录的文件
        """
        L = []
        if type(filetype) == str:
                if len([filetype]) == 1:
                    filelist = os.listdir(self.path)
                    for file in filelist:
                        if os.path.splitext(file)[1] == filetype:
                            L.append(os.path.join(self.path, file))
        if type(filetype) != str:
                if len(filetype) > 1:
                    for i in range(len(filetype)):
                        filelist = os.listdir(self.path)
                        for file in filelist:
                            if os.path.splitext(file)[1] == type[i]:
                                L.append(os.path.join(self.path, file))
        return L
    
    
    def finded_2_3_4_band(self):
        filelist = self.get_tif_file()
        for file in filelist:
            if 'B2.TIF' in file:
                B2file = file
            elif  'B3.TIF' in file:
                B3file = file
            elif  'B4.TIF' in file:
                B4file = file
        return [B2file, B3file, B4file]

ok，调用这个类，这个类的方法有两个finded_2_3_4_band，get_tif_file

我们解析一下目前这个类的运行逻辑。

首先实例化这个类。

其次调用finded_2_3_4_band方法。

然后finded_2_3_4_band方法 ，会调用get_tif_file这个方法。

最后返回一个列表，这个列表包含了B2file, B3file, B4file三个元素。

path = r'X:\LC09_L2SP_121038_20220616_20230411_02_T1'
filelist = landsat9_to_webp(path).finded_2_3_4_band()

我们把filelist循环打印出来，瞧瞧它到底长啥样。

import os
for file in filelist:
    print(os.path.basename(file))

控制台打印一下信息：

LC09_L2SP_121038_20220616_20230411_02_T1_SR_B2.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B3.TIF
LC09_L2SP_121038_20220616_20230411_02_T1_SR_B4.TIF

完成一个功能了，通过文件夹获取Landsat9影像的第2、3、4波段的文件名

GDAL 读取tif

得到对应的文件名后，我们通过gdal把tif读取到内存里。

利用numpy创建一个为【X,Y,4】的矩阵。

将对应的RGB波段数据放到矩阵的第1维到第3维。

第4维放置 透明度。

最后，返回矩阵。

    ds1 = gdal.Open(self.bluefile)
    blue = ds1.ReadAsArray()

    ds2 = gdal.Open(self.greenfile)
    green = ds2.ReadAsArray()

    ds3 = gdal.Open(self.redfile)
    red = ds3.ReadAsArray()

    x, y = blue.shape
    new_arr = np.zeros(shape=[x, y, 4]) # 构造一个新的零矩阵
    new_arr[:, :, 0] = linear_stretch(blue) # 数据拉伸
    new_arr[:, :, 1] = linear_stretch(green)
    new_arr[:, :, 2] = linear_stretch(red)

    new_arr[:, :, 3][blue != 0] = 255  # 透明度.

我们在生成new_arr的时候，使用了

new_arr = np.zeros(shape=[x, y, 4]) # 构造一个新的零矩阵

所以new_arr是一个长为x,宽为y，维度（深度）为4的矩阵，并且全部数值均为0。

透明度解析

python的索引是从0开始计数，

所以new_arr[:, :, 3]代表new_arr的第4维。

new_arr[:, :, 3][blue != 0] = 255  # 透明度.

以上代码的意思是，标记所有blue矩阵中不等于0像素的位置，暂时赋予给new_arr的第4维，将这些位置上的数值改写为255。

简单来说，就是将蓝色波段中存在不为零的像素的透明度设置为 不透明，将蓝色波段中存在为零的像素的透明度设置为透明。

数据拉伸

为什么要做数据拉伸，我曾经写过一篇公众号进行解释。

这里简单说一下，因为Landsat9的数据类型是16bit，而计算机屏幕显示的范围是8bit，16bit转为8bit这个过程叫做数据拉伸。

线性拉伸，也是等比例拉伸，具体代码如下：

def linear_stretch(data, num=1):
    '''

    @param data: 待拉伸的矩阵
    @param num: 拉伸系数，一般为1-5
    @return: 拉伸后的矩阵
    '''
    x, y = np.shape(data)
    data_8bit = data
    data_8bit[data_8bit == -9999] = 0

    # 把数据中的nan转为某个具体数值，例如
    # data_8bit[np.isnan(data_8bit)] = 0
    d2 = np.percentile(data_8bit, num)
    u98 = np.percentile(data_8bit, 100 - num)

    maxout = 255
    minout = 0
    data_8bit_new = minout + ((data_8bit - d2) / (u98 - d2)) * (maxout - minout)
    data_8bit_new[data_8bit_new < minout] = minout
    data_8bit_new[data_8bit_new > maxout] = maxout
    data_8bit_new = data_8bit_new.astype(np.int32)
    return data_8bit_new

整合代码

项目地址是在：

https://github.com/ytkz11/landsat9_to_rgba

如果对你有帮助，点赞是我最大的创作动力，

image-20231010182504647