3. 遥感图像数据读取与元数据解析

拿到一景遥感数据，真正决定后续预处理能否做对的，往往不是图像本身，而是它“附带的那堆说明”。同一幅 Landsat 影像，GeoTIFF 里可能只存了 DN 值和少量标签，而 HDF 里则把每个波段、定标系数、观测几何、质量波段分门别类放在不同数据集里。如果读错了波段顺序、漏掉了定标增益、忽略了填充值，后面的辐射定标和大气校正就会整体偏移。因此本节的核心任务是：认识常见格式的组织方式，学会从头文件和元数据中提取真正参与计算的字段，并把这些字段正确映射到波段数据上。

一、四种常见格式的组织差异

GeoTIFF 是最通用的栅格格式，本质是 TIFF 加上地理标签。它的波段通常按“波段交错”方式存储，常见三种布局：BSQ（按波段分块）、BIL（按行交错）、BIP（按像元交错）。读取时若把 BIP 当成 BSQ 解释，图像会出现规律性条纹。GeoTIFF 的地理信息放在 GeoKey 标签中，包括投影坐标系、仿射变换六参数、NoData 值等，用 GDAL 可直接读取。

HDF（HDF4/HDF5）是分层数据格式，MODIS、VIIRS 等产品大量使用。它像一棵目录树，根下挂多个数据集（Dataset）和属性（Attribute）。例如 MOD09GA 中，sur_refl_b01 到 sur_refl_b07 是七个反射率波段，state_1km 是质量状态位，SensorZenith 是观测天顶角。HDF 的波段不是“第几个”，而是“叫什么名字”，必须按数据集名称索引。

ENVI 格式由一对文件组成：数据文件（如 .img 或 .dat）加头文件 .hdr。头文件是纯文本，字段清晰，例如 samples、lines、bands、data type、interleave、byte order、map info、wavelength。ENVI 的优点是元数据可读可改，缺点是头文件与数据文件必须同名同目录，缺一不可。

PDS（行星数据系统）格式常见于月球、火星探测数据，采用“标签 + 数据”结构，标签为 ASCII 键值对，描述数据对象的行列数、采样方式、缩放因子等。PDS4 进一步用 XML 描述，读取时需先解析标签再定位数据块。

二、头文件与元数据字段的读取

无论哪种格式，读取流程都可以归纳为“先读元数据，再读数据”。以 ENVI 为例，用 Python 读取头文件的关键字段：

import numpy as np

# 解析 ENVI .hdr 头文件
hdr = {}
with open('LC08.hdr', 'r') as f:
    for line in f:
        if '=' in line:
            k, v = line.split('=', 1)
            hdr[k.strip().lower()] = v.strip()

samples = int(hdr['samples'])      # 列数
lines   = int(hdr['lines'])        # 行数
bands   = int(hdr['bands'])        # 波段数
dtype   = hdr['data type']         # 数据类型编码
interleave = hdr['interleave']     # bsq / bil / bip

# 按 BSQ 读取：形状为 (bands, lines, samples)
data = np.fromfile('LC08.img', dtype=np.float32)
data = data.reshape(bands, lines, samples)
print(data.shape, data[0].mean())  # 查看第一波段均值

这里 data type 是编码值，ENVI 中 1 表示 8 位无符号、2 表示 16 位有符号、4 表示 32 位浮点，必须查表转换，不能直接当字节数用。byte order 为 0 表示小端、1 表示大端，读错会导致数值完全错乱。

GeoTIFF 用 GDAL 读取更省事：

from osgeo import gdal
ds = gdal.Open('LC08_B4.tif')
gt = ds.GetGeoTransform()   # (x0, dx, 0, y0, 0, dy)
proj = ds.GetProjection()   # WKT 投影字符串
band = ds.GetRasterBand(1)
nodata = band.GetNoDataValue()
arr = band.ReadAsArray()

其中 GetGeoTransform() 返回的六参数定义了像元到地理坐标的仿射关系：x = x0 + col*dx + row*0，y = y0 + col*0 + row*dy。注意 dy 通常为负值，表示行号增大时纬度减小，这是栅格坐标与地理坐标方向相反的约定，全文保持一致。

三、波段组织与辅助数据的使用

多光谱数据常按波长分组。以 Landsat 8 OLI 为例，波段 1 到 7 为可见光近红外，波段 8 为全色，波段 9 为卷云，波段 10、11 为热红外。读取时要确认波段索引与物理波段的对应关系，不能默认“第 1 个就是蓝光”。HDF 产品则通过数据集名称区分，如 MODIS 的 sur_refl_b03 对应蓝光。

辅助数据同样关键。质量波段（QA）用位标志记录云、云影、雪、水体等状态，使用时需按位解析：

qa = ds_qa.ReadAsArray()
cloud_bit = (qa >> 10) & 1   # 示例：第 10 位表示云
clear = cloud_bit == 0
print('晴空像元比例:', clear.mean())

观测几何数据（太阳天顶角、方位角、传感器天顶角）用于大气校正和 BRDF 归一化。定标系数（增益 gain、偏移 offset）用于把 DN 转为辐亮度：L = gain * DN + offset。这些字段通常存放在元数据属性中，读取后应单独保存为变量，避免与图像数组混淆。

四、读取结果的检查与判断

数据读进来后，不要急着做校正，先做三项检查。第一，看形状：data.shape 是否等于 (bands, lines, samples)，若维度顺序不对，说明 interleave 判断有误。第二，看数值范围：16 位无符号数据应在 0 到 65535 之间，若出现负数或极大值，多半是 byte order 或 data type 读错。第三，看 NoData：统计等于 NoData 值的像元比例，若比例异常高，可能是填充值设置错误。

以 Landsat 8 为例，读取后可用 data[3].mean() 查看红光波段均值，正常景应在数千量级；若接近 0 或超过 60000，需回头检查定标和数据类型。对于 HDF 产品，还要确认数据集名称是否与文档一致，例如 MOD09GA 的 sur_refl_b01 不能写成 sur_refl_b1。

掌握格式差异、字段含义和波段映射关系后，读取就不再是“打开文件”这么简单，而是为后续每一步预处理建立正确的数据基础。元数据解析得越准确，辐射定标、几何校正和大气校正的误差来源就越少。