4、高光谱数据读取：ENVI格式解析、HDF5文件读取、数据可视化

好，咱们进入实战环节的第一个硬骨头——数据读取。

说实话，我刚开始接触高光谱数据时，被各种文件格式搞得一头雾水。什么ENVI、HDF5、GeoTIFF……每个格式都有自己的脾气。你想想看，一个高光谱图像文件，动辄几百兆甚至几个G，里面不光有图像数据，还带着波长信息、头文件、元数据。读不对，后面全白搭。

这一章，我就带你搞定两种最常见的格式：ENVI格式和HDF5格式。顺便把数据可视化也一并拿下。

4.1 ENVI格式解析——老牌遥感格式的硬核玩法

ENVI格式在遥感圈子里用了二十多年了。它的结构其实很简单：一个头文件（.hdr）加一个数据文件。但坑就坑在，数据文件的后缀名五花八门——.dat、.img、.raw、.bsq……甚至没有后缀。

我当年第一次处理AVIRIS数据时，拿到一个叫"f970619t01p02_r02"的文件，没有后缀，差点以为是损坏的。后来一查头文件才知道，这是标准的ENVI格式。

4.1.1 头文件里藏着什么？

头文件是纯文本，用记事本就能打开。里面记录了关键信息：

ENVI
description = {AVIRIS scene}
samples = 614
lines = 512
bands = 224
header offset = 0
file type = ENVI Standard
data type = 2           ; 2=16位有符号整型
interleave = bsq        ; 波段顺序：BSQ/BIL/BIP
byte order = 0          ; 0=小端序，1=大端序
wavelength = {366.0, 376.0, ...}

这里有几个关键参数，我重点说一下：

• samples/lines/bands：图像的宽、高、波段数。224个波段，每个波段就是一张图。
• interleave：波段存储顺序。BSQ是每个波段单独存，BIL是按行交叉，BIP是按像素交叉。这个决定了你读取时怎么reshape数据。
• data type：数据类型。1=字节，2=16位整型，4=32位浮点。读错了类型，图像会花掉。
• header offset：头偏移量。有些文件在数据前面会加一些头信息，这个值告诉你要跳过多少字节。

4.1.2 用Python读取ENVI格式

我推荐两种方式：一是用spectral库，封装得比较好；二是手动解析，适合定制化需求。

方式一：用spectral库（推荐）

import spectral
from spectral import open_image

# 读取ENVI文件（自动解析头文件）
img = open_image('scene.hdr')

# 查看基本信息
print(img.shape)        # (lines, samples, bands)
print(img.bands)        # 波段数
print(img.metadata)     # 头文件元数据

# 加载全部数据到内存
data = img.load()
print(data.shape)       # (512, 614, 224)

方式二：手动解析（适合定制需求）

import numpy as np
import struct

def read_envi_raw(hdr_file, data_file):
    # 解析头文件
    metadata = {}
    with open(hdr_file, 'r') as f:
        for line in f:
            if '=' in line:
                key, value = line.split('=', 1)
                metadata[key.strip()] = value.strip()
    
    samples = int(metadata['samples'])
    lines = int(metadata['lines'])
    bands = int(metadata['bands'])
    offset = int(metadata.get('header offset', 0))
    
    # 根据data type确定字节数
    dtype_map = {1: np.uint8, 2: np.int16, 4: np.float32}
    dtype = dtype_map[int(metadata['data type'])]
    
    # 读取二进制数据
    with open(data_file, 'rb') as f:
        f.seek(offset)
        raw_data = np.fromfile(f, dtype=dtype)
    
    # 根据interleave reshape
    interleave = metadata.get('interleave', 'bsq').lower()
    if interleave == 'bsq':
        data = raw_data.reshape((bands, lines, samples))
        data = np.transpose(data, (1, 2, 0))  # 转为 (lines, samples, bands)
    elif interleave == 'bil':
        data = raw_data.reshape((lines, bands, samples))
        data = np.transpose(data, (0, 2, 1))
    elif interleave == 'bip':
        data = raw_data.reshape((lines, samples, bands))
    
    return data, metadata

4.2 HDF5文件读取——大规模数据的首选

HDF5是Hierarchical Data Format的缩写。说白了，它就像一个文件系统，里面可以放数据集（dataset）和组（group）。NASA的很多卫星数据，比如MODIS、VIIRS，都用这个格式。

为什么用HDF5？因为高光谱数据太大了。一个224波段的AVIRIS场景，如果存成ENVI格式，光数据文件就500MB以上。HDF5支持压缩、分块读取，还能存元数据，特别适合大规模数据。

4.2.1 用h5py读取HDF5文件

import h5py
import numpy as np

# 打开HDF5文件
with h5py.File('hyperspectral_data.h5', 'r') as f:
    # 查看文件结构
    def print_structure(name, obj):
        print(name, type(obj))
    f.visititems(print_structure)
    
    # 读取数据集
    # 假设数据集路径为 '/data/hsi_cube'
    hsi_data = f['/data/hsi_cube'][:]
    wavelengths = f['/data/wavelengths'][:]
    
    print(f"数据形状: {hsi_data.shape}")
    print(f"波长范围: {wavelengths[0]:.1f} - {wavelengths[-1]:.1f} nm")

这里有个关键点：HDF5文件里的数据集路径要搞清楚。不同数据源的组织方式不一样。我处理过的一个国产高光谱数据，数据藏在/Measurement/SpectralCube下面，找了半天。

4.2.2 读取属性元数据

HDF5的每个数据集和组都可以挂属性（attributes）。这些属性里往往藏着关键信息，比如波长、增益、偏移量。

with h5py.File('hyperspectral_data.h5', 'r') as f:
    dset = f['/data/hsi_cube']
    
    # 读取属性
    for key in dset.attrs:
        print(f"{key}: {dset.attrs[key]}")
    
    # 常见属性
    if 'wavelength' in dset.attrs:
        wavelengths = dset.attrs['wavelength']
    elif 'wavelengths' in dset.attrs:
        wavelengths = dset.attrs['wavelengths']

4.3 数据可视化——让高光谱数据"看得见"

数据读进来了，下一步就是看看它长什么样。高光谱数据有几十上百个波段，不能像普通RGB图像那样直接显示。我一般从三个角度入手：单波段、伪彩色、光谱曲线。

4.3.1 单波段灰度显示

先挑一个波段看看。比如近红外波段（第50波段左右）通常植被很亮，水体很暗。

import matplotlib.pyplot as plt

# 假设 data 是 (lines, samples, bands) 的numpy数组
band_index = 50  # 第50波段
band_data = data[:, :, band_index]

plt.figure(figsize=(8, 6))
plt.imshow(band_data, cmap='gray')
plt.colorbar(label='DN值')
plt.title(f'波段 {band_index+1} 灰度图')
plt.show()

4.3.2 伪彩色合成

选三个波段，分别赋给R、G、B通道。比如用近红外、红、绿波段做标准假彩色合成，植被会显示为红色。

# 标准假彩色合成：R=近红外，G=红，B=绿
red_band = data[:, :, 29]    # 红波段（约650nm）
green_band = data[:, :, 19]  # 绿波段（约550nm）
nir_band = data[:, :, 49]    # 近红外波段（约800nm）

# 归一化到0-1范围
def normalize(band):
    band_min = band.min()
    band_max = band.max()
    return (band - band_min) / (band_max - band_min)

rgb_image = np.stack([
    normalize(nir_band),   # R通道
    normalize(red_band),   # G通道
    normalize(green_band)  # B通道
], axis=2)

plt.figure(figsize=(8, 6))
plt.imshow(rgb_image)
plt.title('标准假彩色合成 (NIR-R-G)')
plt.show()

4.3.3 光谱曲线绘制

选几个像素点，画出它们的光谱曲线。这是高光谱分析最基础的操作。

# 选取几个典型像素
pixels = {
    '植被': (200, 300),
    '水体': (100, 150),
    '土壤': (350, 400)
}

plt.figure(figsize=(10, 6))
for label, (row, col) in pixels.items():
    spectrum = data[row, col, :]
    plt.plot(wavelengths, spectrum, label=label, linewidth=2)

plt.xlabel('波长 (nm)')
plt.ylabel('反射率')
plt.title('典型地物光谱曲线')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

4.4 实战小贴士

最后，分享几个我在项目中踩过的坑：

• 内存爆炸：一个224波段、1024x1024像素的数据，float32类型，占内存约900MB。如果电脑内存不够，用np.memmap做内存映射，或者分块读取。
• 波长单位：有的数据用纳米（nm），有的用微米（μm）。读出来先检查一下，别把单位搞混了。
• 数据缩放因子：很多高光谱数据存的是DN值，需要乘以缩放因子才能得到反射率。这个信息通常在头文件或HDF5属性里。

嗯，数据读取这块就讲到这里。你把这些代码跑一遍，基本就能搞定大部分高光谱数据了。下一节，咱们聊聊数据预处理——坏波段剔除、大气校正这些实战操作。

---

公众号:蓝海资料掘金营，微信deep3321