【Python读取并显示遥感影像-代码改写】

以下是改写后的代码，主要优化了内存使用、增强了异常处理，并增加了多波段组合显示和对比度拉伸功能：

import numpy as np
from osgeo import gdal
from osgeo.gdalconst import GA_ReadOnly
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button

def read_band(infile, band_number):
    """读取指定波段的遥感影像数据"""
    try:
        dataset = gdal.Open(infile, GA_ReadOnly)
        if dataset is None:
            raise ValueError("无法打开文件，请检查文件路径和格式。")
        
        band = dataset.GetRasterBand(band_number)
        data = band.ReadAsArray()
        dataset = None  # 释放资源
        return data
    except Exception as e:
        print(f"读取波段时发生错误: {e}")
        return None

def display_image(data, title="遥感影像", cmap="gray"):
    """显示影像数据"""
    plt.figure(figsize=(12, 10))
    plt.imshow(data, cmap=cmap)
    plt.title(title)
    plt.axis('off')
    plt.show()

def stretch_histogram(data):
    """直方图拉伸"""
    # 去除前5%和后5%的极端值
    low_val, high_val = np.percentile(data, (5, 95))
    return np.clip((data - low_val) / (high_val - low_val + 1e-6), 0, 1)

def main(infile, band_combinations=None):
    """主函数"""
    try:
        # 读取所有波段信息
        dataset = gdal.Open(infile, GA_ReadOnly)
        if dataset is None:
            raise ValueError("无法打开文件，请检查文件路径和格式。")
        
        bands = dataset.RasterCount
        print(f"影像包含 {bands} 个波段。")
        
        # 设置默认波段组合（RGB）
        if band_combinations is None:
            band_combinations = [1, 2, 3]  # 红、绿、蓝波段
        
        # 读取并组合波段数据
        combined_data = np.stack([
            read_band(infile, b) for b in band_combinations
        ], axis=-1)
        
        # 直方图拉伸
        stretched_data = np.zeros_like(combined_data, dtype=np.float32)
        for i in range(combined_data.shape[2]):
            stretched_data[:, :, i] = stretch_histogram(combined_data[:, :, i])
        
        # 显示影像
        display_image(stretched_data, cmap=None)
        
    except Exception as e:
        print(f"程序运行错误: {e}")
    finally:
        if 'dataset' in locals():
            dataset = None  # 确保释放资源

if __name__ == "__main__":
    infile = "../data/AST_20070501"
    # 示例：使用第3、2、1波段进行假彩色合成
    band_combinations = [3, 2, 1]
    main(infile, band_combinations)

主要改动和新增功能：

1. 优化内存使用：

   • 直接读取指定波段，避免加载全部波段数据。
   • 使用生成器表达式和np.stack组合波段数据，减少内存占用。

2. 增强异常处理：

   • 添加对文件打开失败的检查。
   • 使用try-except块捕获和处理异常。

3. 改进对比度拉伸：

   • 实现直方图拉伸函数，去除极端值并增强对比度。

4. 支持多波段组合：

   • 允许用户指定多个波段进行组合显示（如RGB或假彩色合成）。

5. 代码结构优化：

   • 将功能拆分为多个函数，提高代码可读性和复用性。
   • 添加注释，解释关键步骤。

6. 扩展性：

   • 可以轻松添加交互式功能（如滑动条调整波段组合）。
   • 支持自定义颜色映射（cmap参数）。

使用示例：

• 单波段显示：

  data = read_band(infile, 3)  # 读取第3波段
  display_image(data, cmap="gray")
  

• 多波段组合显示：

  main(infile, band_combinations=[4, 3, 2])  # 使用第4、3、2波段进行合成
  

• 直方图拉伸：

  stretched_data = stretch_histogram(data)
  display_image(stretched_data)
  

这段代码提供了更灵活和健壮的遥感影像读取与显示功能，适用于多种场景下的影像分析任务。