用Python和GDAL处理高分二号卫星遥感数据：从TIF读取到归一化的保姆级教程

用Python和GDAL处理高分二号卫星遥感数据：从TIF读取到归一化的保姆级教程

遥感数据处理是地理信息科学和深度学习交叉领域的关键环节。对于刚接触卫星影像分析的开发者而言，如何正确解析多波段TIF文件并将其转化为适合神经网络训练的输入格式，往往成为项目推进的第一个技术门槛。本文将构建一套完整的处理流程，特别针对高分二号卫星数据的特点，从底层数据解析到实用预处理技巧，逐步拆解每个技术细节。

1. 卫星遥感数据的特殊性认知

与普通JPG/PNG图像不同，卫星遥感数据是以科学数据格式存储的多维信息集合。以高分二号多光谱数据为例，其TIFF文件通常包含四个波段：

这些波段的数值范围差异显著：

• 可见光波段（1-3）通常取值在0-2000DN值
• 近红外波段（4）可能达到3000-4000DN值

import gdal import numpy as np # 查看波段数值范围示例 dataset = gdal.Open('GF2.tif') for i in range(1, 5): band = dataset.GetRasterBand(i) arr = band.ReadAsArray() print(f'波段{i} - 最小值:{arr.min()} 最大值:{arr.max()}')

注意：直接使用原始DN值进行可视化或模型训练会导致结果失真，必须经过专业预处理

2. GDAL高效读取多波段数据

GDAL作为地理空间数据的"瑞士军刀"，其Python绑定提供了直接操作卫星数据的能力。以下是优化后的读取流程：

def read_tif_advanced(img_path, bands_order=[3,2,1], nodata=None): """ 增强版TIF读取函数 :param img_path: 文件路径 :param bands_order: 波段顺序，默认[3,2,1]对应RGB :param nodata: 自定义无效值处理 :return: 三维数组 + 元数据字典 """ dataset = gdal.Open(img_path) if not dataset: raise ValueError("文件无法打开，请检查路径") # 获取基础信息 cols = dataset.RasterXSize rows = dataset.RasterYSize geotrans = dataset.GetGeoTransform() projection = dataset.GetProjection() # 多波段读取优化 data = np.empty((rows, cols, len(bands_order)), dtype=np.float32) for i, band_idx in enumerate(bands_order): band = dataset.GetRasterBand(band_idx) if nodata is not None: band.SetNoDataValue(nodata) arr = band.ReadAsArray() data[:,:,i] = arr meta = { 'geotransform': geotrans, 'projection': projection, 'driver': dataset.GetDriver().ShortName } return data, meta

关键改进点：

1. 增加异常处理机制
2. 保留地理参考信息
3. 支持无效值处理
4. 使用float32节省内存

3. 面向深度学习的归一化策略

传统图像处理常用的/255归一化方法完全不适用于遥感数据。我们推荐两种科学归一化方法：

3.1 百分比截断拉伸法

def percentile_normalization(data, lower=2, upper=98): """ 基于百分位的自适应归一化 :param data: 输入数组 :param lower: 下限百分比 :param upper: 上限百分比 :return: 归一化后的数组 """ per_low = np.percentile(data, lower) per_high = np.percentile(data, upper) norm_data = (data - per_low) / (per_high - per_low) norm_data = np.clip(norm_data, 0, 1) return norm_data # 多波段分别处理示例 rgb_data = read_tif_advanced('GF2.tif')[0] for i in range(3): rgb_data[:,:,i] = percentile_normalization(rgb_data[:,:,i])

3.2 标准差归一化法

def sigma_normalization(data, sigma=2): """ 基于标准差的正态分布归一化 :param data: 输入数组 :param sigma: 标准差倍数 :return: 归一化后的数组 """ mean = np.mean(data) std = np.std(data) norm_data = (data - mean) / (sigma * std) norm_data = (norm_data + 1) / 2 # 映射到0-1范围 norm_data = np.clip(norm_data, 0, 1) return norm_data

两种方法对比：

4. 多波段数据可视化技巧

将多波段数据转换为JPG需要特殊的波段组合策略。以下是常见可视化方案：

def visualize_band_combination(data, combination='natural_color'): """ 专业级波段组合可视化 :param data: 原始多波段数据 :param combination: 波段组合模式 :return: 可视化数组(0-255) """ band_combinations = { 'natural_color': [3,2,1], # 真彩色 'vegetation': [4,3,2], # 植被增强 'urban': [4,2,1], # 城市特征 'water': [2,3,4] # 水体特征 } if combination not in band_combinations: raise ValueError(f"不支持的组合模式，可选: {list(band_combinations.keys())}") selected = band_combinations[combination] vis_data = data[:,:,selected] # 对各波段单独拉伸 for i in range(3): vis_data[:,:,i] = percentile_normalization(vis_data[:,:,i]) return (vis_data * 255).astype(np.uint8)

典型应用场景：

• 自然资源监测：使用[4,3,2]组合突出植被
• 水体识别：采用[2,3,4]组合增强水体特征
• 城市扩张分析：[4,2,1]组合适合建筑识别

5. 生产环境优化实践

在实际项目中，我们需要考虑大规模数据处理的效率和稳定性：

5.1 内存映射处理大文件

def process_large_tif(input_path, output_path, chunk_size=1024): """ 分块处理大尺寸TIF文件 :param input_path: 输入文件路径 :param output_path: 输出文件路径 :param chunk_size: 分块大小 """ dataset = gdal.Open(input_path) cols = dataset.RasterXSize rows = dataset.RasterYSize driver = gdal.GetDriverByName('GTiff') out_dataset = driver.Create(output_path, cols, rows, 3, gdal.GDT_Float32) for i in range(0, rows, chunk_size): # 计算当前块的实际高度 current_chunk = min(chunk_size, rows - i) # 读取数据块 data = np.empty((current_chunk, cols, 3)) for b in range(3): band = dataset.GetRasterBand(b+1) data[:,:,b] = band.ReadAsArray(0, i, cols, current_chunk) # 处理数据块 processed = percentile_normalization(data) # 写入输出文件 for b in range(3): out_band = out_dataset.GetRasterBand(b+1) out_band.WriteArray(processed[:,:,b], 0, i) out_dataset.FlushCache()

5.2 并行处理加速

from multiprocessing import Pool def parallel_convert(file_list): """ 多进程并行转换 :param file_list: 文件路径列表 """ with Pool(processes=4) as pool: pool.map(process_single_file, file_list) def process_single_file(file_path): """单个文件处理函数""" try: data = read_tif_advanced(file_path) norm_data = percentile_normalization(data) output_path = file_path.replace('.tif', '_norm.tif') save_as_tif(norm_data, output_path) except Exception as e: print(f"处理{file_path}时出错: {str(e)}")

6. 与深度学习框架的集成

将处理好的数据输入PyTorch/TensorFlow模型时，推荐使用自定义Dataset类：

import torch from torch.utils.data import Dataset class SatelliteDataset(Dataset): def __init__(self, file_list, transform=None): self.file_list = file_list self.transform = transform def __len__(self): return len(self.file_list) def __getitem__(self, idx): img_path = self.file_list[idx] # 读取并预处理 data, _ = read_tif_advanced(img_path) norm_data = percentile_normalization(data) # 转为Tensor tensor_data = torch.from_numpy(norm_data).float() tensor_data = tensor_data.permute(2, 0, 1) # HWC -> CHW if self.transform: tensor_data = self.transform(tensor_data) return tensor_data

提示：在数据加载时使用torch.utils.data.DataLoader可实现自动批处理和异步加载

实际项目中，建议将原始数据、预处理代码和模型训练整合为完整pipeline：

# 完整训练示例 train_files = [...] # 训练集文件列表 val_files = [...] # 验证集文件列表 train_dataset = SatelliteDataset(train_files, transform=augmentation_transform) val_dataset = SatelliteDataset(val_files) train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True) val_loader = DataLoader(val_dataset, batch_size=32) model = create_model() optimizer = torch.optim.Adam(model.parameters()) for epoch in range(100): train_one_epoch(model, train_loader, optimizer) validate(model, val_loader)