8. 4波段遥感影像获取前3波段和计算添加ndvi

gdal读取了遥感数据后就变成了numpy数组了，之后就是numpy的操作了，numpy操作了之后再加坐标保存就行了，实际上很简单读取和写入看下面的连接就够了
(4条消息) 5. gdal实现对遥感影像的读写，信息统计和波段选择_xiaotiig的博客-CSDN博客
https://blog.csdn.net/xiaotiig/article/details/118721427

1 提取4波段的前3波段

# 将4波段的遥感影像提取出前3波段

import os
import gdal
import numpy as np

#  读取tif数据集
def readTif(fileName):
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName + "文件无法打开")
    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")
    # 这个写入函数比较特殊，跟数据的shape不一样，不一样才对，数据对应行列，这里必须是列行，先写列，再写行
    # https://vimsky.com/examples/detail/python-method-gdal.GetDriverByName.html
    dataset = driver.Create(path, int(im_width), int(im_height), int(3), datatype)
    if (dataset != None):
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影
    # 这里写入前3个波段，全部输出im_bands
    for i in range(3):
        dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset


def bands4to3(TifPath, SavePath):
    i = 0
    filenames = os.listdir(TifPath)
    print(len(filenames))
    for filename in filenames:
        i +=1
        if i%1000==0:
            print("number：%d"%i)
        # 输出文件名
        file_path = os.path.join(TifPath, filename)
        # 读取栅格文件
        dataset_img = readTif(file_path)
        proj = dataset_img.GetProjection()
        geotrans = dataset_img.GetGeoTransform()
        img = dataset_img.ReadAsArray()  # 获取数据
        writeTiff(img, geotrans, proj, os.path.join(SavePath, filename))

if __name__ == '__main__':
    bands4to3(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\2arcgislLuotian\3imageDataset\1lt497908\images",
              r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\2arcgislLuotian\3imageDataset\1lt497908\images3"
              )

2 进行波段的组合计算

在4波段的基础上加入ndvi波段，如果之后要保存为tif文件，看
(4条消息) 9.将数组存储为带坐标的tif文件_xiaotiig的博客-CSDN博客
https://blog.csdn.net/xiaotiig/article/details/122308138

def imgread(fileName, addNDVI=True):
    dataset = gdal.Open(fileName)
    width = dataset.RasterXSize
    height = dataset.RasterYSize
    data = dataset.ReadAsArray(0, 0, width, height)   # data, (4, 36786, 37239) ,波段，行，列
    # 如果是image的话,才进行处理，因为label是单通道，
    if(len(data.shape) == 3):
        # 添加归一化植被指数NDVI特征
        if(addNDVI):
            # 这景影像的波段顺序为，R,G,B，NIR
            nir, r = data[3], data[0]
            ndvi = (nir - r) / (nir + r + 0.00001) * 1.0    # 这里已经计算好了ndvi了
            # 和其他波段保持统一,归到0-255,在深度学习中后面的totensor会/255统一归一化
            # -1到1，先变成0到2，再乘
            ndvi = (ndvi+1) * 127.5
            # 去除小于0和大于255的值,这步不需要
            # ndvi = np.clip(ndvi, 0, 255)
            data_add_ndvi = np.zeros((5, 256, 256), np.uint8)
            data_add_ndvi[0:4] = data
            # 先进行四舍五入，再转换为8位深度
            data_add_ndvi[4] = np.uint8(np.around(ndvi))
            data = data_add_ndvi
        # (C,H,W)->(H,W,C) 这里进行了波段转换，为了深度学习处理，如果单纯的要保存，不用转换用(C,H,W)的数组
        data = data.swapaxes(1, 0).swapaxes(1, 2)
    return data