Python进行.tif格式数据波段组合；多波段tif数据波段分离

阔别已久，我又来了；啊啊啊……
读了研究生，忙到天昏地暗（哈哈哈，这都是借口，主要是花了一定时间适应研究生生活）

这方面的内容我也主要是应用在遥感领域，我们知道在遥感应用中，大多数的图像数据格式为.tif格式的。进行波段组合使用ArcGIS里面的Composite Bands工具是可以实现的，波段分离也可以使用一些工具或巧妙的方法分开。但当你使用Python程序的话，就可以简化这个软件点击和等待的过程，只要修改2个路径就可以了，岂不乐载？

.tif格式数据波段组合

import os
 
import numpy as np
from osgeo import gdal
from skimage._shared.utils import safe_as_int
 
class IMAGE:
 
    # 读图像文件
    def read_img(self,filename):
 
        dataset = gdal.Open(filename)  # 打开文件
 
        im_width = dataset.RasterXSize  # 栅格矩阵的列数
        im_height = dataset.RasterYSize  # 栅格矩阵的行数
        # im_bands = dataset.RasterCount  # 波段数
        im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率
        im_proj = dataset.GetProjection()  # 地图投影信息，字符串表示
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
 
        del dataset   #关闭对象dataset，释放内存
        # return im_width, im_height, im_proj, im_geotrans, im_data,im_bands
        return  im_proj, im_geotrans, im_data, im_width,im_height
 
    # 遥感影像的存储
    # 写GeoTiff文件
    def write_img(self,filename, im_proj, im_geotrans, im_data):
        # 判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
 
        # 判读数组维数
        if len(im_data.shape) == 3:
            # 注意数据的存储波段顺序：im_bands, im_height, im_width
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape   #没看懂
 
        # 创建文件时 driver = gdal.GetDriverByName("GTiff")，数据类型必须要指定，因为要计算需要多大内存空间。
        driver = gdal.GetDriverByName("GTiff")
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
 
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影
 
        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
 
        del dataset
 
if __name__ == "__main__":
    # os.chdir(r'')   #切换路径到待处理图像所在文件夹
    run=IMAGE()
    # 第一步
    proj, geotrans, data1, row1, column1  = run.read_img(r'***.tif')  # 读数据   
    proj, geotrans, data2, row2, column2  = run.read_img(r'***.tif')  # 读数据
    proj, geotrans, data3, row3, column3  = run.read_img(r'***.tif')  # 读数据

    # 第二步:将上述读取的3个波段放到一个数组中
    data = np.array((data1,data2), dtype=data1.dtype)  # 按序将3个波段像元值放入
    # data = np.array((data1[0,::],data1[1,::],data2), dtype=data1.dtype) #这个是应对要进行波段组合的图像原本就是多波段的特点，因此我们可以用data[0,::],data[1,::],data[2,::]等表示，前面的0，1，2是图像波段号的索引
    # 第三步
    run.write_img(r'***.tif', proj, geotrans, data)  # 写为3波段数据,假彩色，nir,red,green

多波段.tif数据波段分离

import numpy as np 
import os
import random
from osgeo import gdal
# import cv2
import tifffile as tif
from skimage import data,exposure
from sklearn import preprocessing

#  读取图像像素矩阵
#  fileName 图像文件名
def readTif(fileName):
    dataset = gdal.Open(fileName)
    # width = dataset.RasterXSize
    # height = dataset.RasterYSize
    # GdalImg_data = dataset.ReadAsArray(0, 0, width, height)
    # return GdalImg_data
    return dataset

#保存tif文件函数
def writeTiff(im_data,im_geotrans,im_proj,path):
    if 'int8' in im_data.dtype.name:
        datatype=gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype=gdal.GDT_UInt16
    else:
        datatype=gdal.GDT_Float32
    if len(im_data.shape)==3:
        im_bands,im_height,im_width=im_data.shape
    elif len(im_data.shape)==2:
        im_data=np.array([im_data])
        im_bands,im_height,im_width=im_data.shape
    #创建文件
    driver=gdal.GetDriverByName("GTiff")
    dataset=driver.Create(path,int(im_width),int(im_height),int(im_bands),datatype)
    if(dataset!=None): 
        dataset.SetGeoTransform(im_geotrans)#写入仿射变换参数
        dataset.SetProjection(im_proj)#写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(im_data[i])
    del dataset

def bandsclip(path1,path2):
    dataset_img=readTif(path1)
    width=dataset_img.RasterXSize
    height=dataset_img.RasterYSize
    proj=dataset_img.GetProjection()
    geotrans=dataset_img.GetGeoTransform()
    img=dataset_img.ReadAsArray(0,0,width,height)

    img_out=[]
    #依次将各波段输出
    for i in range(img.shape[0]):
        img_out=np.array(img[i,::])
        #保存tiff格式文件数据
        writeTiff(img_out,geotrans,proj,path2+str(i)+'.tif') #输出波段的名称命名格式可以修改，结合传递的path2参数

os.chdir(r'G:/20211216/202109021')

path1=r'202109021.tif'  #要分离波段的原始图像数据名称
path2=r'202109021'      #分离的各波段结果图像部分名称
bandsclip(path1,path2)  #调用上面定义的波段分离函数
print('Bandsclip END!')