使用python对遥感波段进行合成

Question

Landsat8 遥感图像，原始影像为 16bits 数据，选择 将“landsat8 湖泊图像 ”中的 R 、G 、B 和全色波段影像，采用最大最小值拉伸 方式，转换为 8bits 影像（注意 nodata 为 65535，不能作为波段最大值），并将 R 、G 、B 三个波段图像合成为一幅 RGB 图像。

from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np
import os


def read_band(file_path, band_index):
    dataset = gdal.Open(file_path, gdal.GA_ReadOnly)
    band = dataset.GetRasterBand(band_index)
    array = band.ReadAsArray()
    return array


def stretch_8bit(band_array, nodata=65535):
    valid_pixels = band_array[band_array != nodata]
    min_val = valid_pixels.min()
    max_val = valid_pixels.max()
    stretched = ((band_array - min_val) / (max_val - min_val) * 255).astype(np.uint8)
    stretched[band_array == nodata] = 0
    return stretched


# 路径和波段索引可能需要调整
os.chdir(r'D:\2024test\landsat8')
r_band_path = "D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B4.TIF"
g_band_path = "D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B3.TIF"
b_band_path = "D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B2.TIF"

r_band = read_band(r_band_path, 1)  # 假设R波段在第一个波段
g_band = read_band(g_band_path, 2)  # 假设G波段在第二个波段
b_band = read_band(b_band_path, 3)  # 假设B波段在第三个波段

r_stretched = stretch_8bit(r_band)
g_stretched = stretch_8bit(g_band)
b_stretched = stretch_8bit(b_band)

rgb_image = np.dstack((r_stretched, g_stretched, b_stretched))#


# 保存RGB影像
driver = gdal.GetDriverByName('GTiff')
out_dataset = driver.Create('path_to_output_rgb_image.tif', rgb_image.shape[1], rgb_image.shape[0], 3, gdal.GDT_Byte)
for i in range(3):
    out_dataset.GetRasterBand(i + 1).WriteArray(rgb_image[:, :, i])
out_dataset.FlushCache()

出现了报错

ERROR 5: D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B3.TIF: GDALDataset::GetRasterBand(2) - Illegal band #

Traceback (most recent call last):
  File "D:\桌面\python\totaltext\课程\删格2.py", line 30, in <module>
    g_band = read_band(g_band_path, 2)  # 假设G波段在第二个波段
  File "D:\桌面\python\totaltext\课程\删格2.py", line 10, in read_band
    array = band.ReadAsArray()
AttributeError: 'NoneType' object has no attribute 'ReadAsArray'

怎么解决

Answer 1

【以下回答由 GPT 生成】

from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np
import os

def read_band(band_path, band_index):
    dataset = gdal.Open(band_path, gdal.GA_ReadOnly)
    band = dataset.GetRasterBand(band_index).ReadAsArray()
    return band

def stretch_8bit(band):
    min = 0
    max = 255
    nodata = 65535
    band = np.where(band == nodata, 0, band)  # 将nodata值替换为0
    stretched_band = (band - band.min()) * (max - min) / (band.max() - band.min()) + min
    return stretched_band.astype(np.uint8)

# 设置路径和波段索引
os.chdir(r'D:\2024\test\landsat8')
r_band_path = r"D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B4.TIF"
g_band_path = r"D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B3.TIF"
b_band_path = r"D:\桌面\遥感图象课程设计\课程\试验数据\landsat8湖泊图像\LC08_L1TP_121040_20210925_20211001_01_T1_B2.TIF"

r_band = read_band(r_band_path, 1)  
g_band = read_band(g_band_path, 1)  
b_band = read_band(b_band_path, 1)  

# 调用stretch_8bit函数
r_stretched = stretch_8bit(r_band)
g_stretched = stretch_8bit(g_band)
b_stretched = stretch_8bit(b_band)

# 将三个波段图像合成为RGB图像
rgb_image = np.dstack((r_stretched, g_stretched, b_stretched))

# 保存RGB影像
driver = gdal.GetDriverByName('GTiff')
path_to_output_rgb_image = 'path_to_output_rgb_image.tif'
out_dataset = driver.Create(path_to_output_rgb_image, rgb_image.shape[1], rgb_image.shape[0], 3, gdal.GDT_Byte)

for i in range(3):
    out_band = out_dataset.GetRasterBand(i + 1)
    out_band.WriteArray(rgb_image[:,:,i])

out_dataset.FlushCache()
out_dataset = None

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^