Landsat影像的单独的tif文件合成为一个tif影像

今天需要在ENVI显示一景Landsat 9 影像，猛地发现竟然不能在envi直接打开，网上一搜大概是要下载插件，于是乎就打开俺以前写的基础代码，修修改改实现了Landsat 9影像的所有文件合成为一个波段。

把原始的L2级的Landsat 9 数据解压后，文件结构如下：xJgFaD.png

全部代码如下：

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2022/10/9 14:46 
# @File : combine_band.py
# l8影像的所有文件合成为一个波段
import os
from osgeo import gdal  # 非常重要是库, 先了解它的接口，后面玩得六之后再用C++

class Landsat8Combine:
    def __init__(self, targetPath, outPut):
        self.targetPath = targetPath
        self.outPut = outPut
        self.l8_to_tif()

    def type_of_image(self, path):
        # 判断影像是哪个卫星：Landsat 8 、 Sentinel 2、 Gaofen
        #  筛选

        image_fn = []

        type = 0
        for f_name in os.listdir(path):
            if f_name.endswith('.tiff'):
                image_fn.append(f_name)

            elif f_name.endswith('.tif'):
                image_fn.append(f_name)

            elif f_name.endswith('.TIF'):
                image_fn.append(f_name)

        if image_fn[0][:2] == image_fn[1][:2]:
            type = image_fn[0][:2]

        image_fn = []
        img = []
        if type == 'LC':
            for f_name in os.listdir(path):
                if f_name.endswith('.tiff'):

                    if "_B" in f_name:
                        image_fn.append(f_name)
                elif f_name.endswith('.tif'):

                    if "_B" in f_name:
                        image_fn.append(f_name)
                elif f_name.endswith('.TIF'):
                    if "_B" in f_name:
                        image_fn.append(f_name)
            # 排除B10波段
            for f in image_fn:
                if 'B10' not in f:
                    f = os.path.join(path, f)
                    img.append(f)
        return type, img

    def l8_to_tif(self):
        """
        把多个单波段的文件合成为一个新的tif文件
        """
        out_path = self.outPut
        type, img = self.type_of_image(self.targetPath)
        if type =='LC':
            # 设置输出波段
            IDataSet = gdal.Open(img[0])
            cols=IDataSet.RasterXSize
            rows = IDataSet.RasterYSize
            proj = IDataSet.GetProjection()
            Driver = IDataSet.GetDriver()
            geoTransform1 = IDataSet.GetGeoTransform()
            ListgeoTransform = list(geoTransform1)
            # ListgeoTransform[5] = -ListgeoTransform[5]

            newgeoTransform1 = tuple(ListgeoTransform)
            outFileName = os.path.splitext(os.path.basename(img[0]))[0][:-3]
            Outname = os.path.join(out_path, outFileName + ".tiff")

            # 获取数据类型
            self.im_data = IDataSet.ReadAsArray()
            if "int8" in self.im_data.dtype.name:
                datetype = gdal.GDT_Byte
            elif "int16" in self.im_data.dtype.name:
                datetype = gdal.GDT_UInt16
            else:
                datetype = gdal.GDT_Float32

            outDataset = Driver.Create(Outname, cols, rows, 7, datetype)
            outDataset.SetGeoTransform(newgeoTransform1)
            outDataset.SetProjection(proj)

            for m in range(1, 8):
                print('读取第%d波段' % m)
                dataset = gdal.Open(img[m - 1])
                band = dataset.GetRasterBand(1)
                data = band.ReadAsArray(0, 0, cols, rows)
                outband = outDataset.GetRasterBand(m)
                outband.SetNoDataValue(-9999)  # 无效值设置为-9999
                outband.WriteArray(data, 0, 0)
        else:
            print("Is not Landsat 8 image")
        return Outname
if __name__ == '__main__':
    file = '' # 解压后的路径
    out = ''  # 输出路径
    Landsat8Combine(file, out)

打个广告喔，新手阶段的同学可以关注我的个人公众号remote sensing。

xJ5YAP.png