3. 多光谱数据格式与读取：常见数据格式（GeoTIFF、HDF、NetCDF），使用Python（GDAL、Rasterio）读取与可视化

各位同学，今天我们来聊聊多光谱数据的“门牌号”——数据格式与读取。说实话，搞遥感的人，一半时间花在分析上，另一半时间就花在跟数据格式较劲上。我自己刚入行那会儿，拿到一个HDF文件，愣是打不开，急得满头大汗。后来才明白，每种格式都有自己的脾气，摸透了就好办了。

3.1 三大主流格式：GeoTIFF、HDF、NetCDF

多光谱数据不像普通照片，它不光有像素值，还带着地理位置、波段信息、甚至卫星姿态。所以，存储格式必须“能装”。目前主流的就是这三种：

格式	特点	典型应用	我踩过的坑
GeoTIFF	基于TIFF，嵌入地理坐标信息	Landsat、高分系列、无人机影像	坐标参考系统（CRS）不一致，叠加时错位
HDF	分层数据格式，支持多维数组	MODIS、VIIRS、部分国产卫星	子数据集（Subdataset）藏得深，容易漏掉
NetCDF	网络通用数据格式，自描述性强	气象卫星、海洋遥感、气候模型	维度顺序搞反，读取出来全是乱码

你想想看，GeoTIFF就像一本带地图的相册，翻开就能看到照片在哪拍的。HDF更像一个多层抽屉的柜子，每层放着不同的数据。NetCDF呢，有点像带说明书的数据包，每个变量都标注得清清楚楚。

3.2 用GDAL读取数据——老牌工具，稳如老狗

GDAL（Geospatial Data Abstraction Library）是遥感界的“瑞士军刀”。我个人的习惯是，不管什么格式，先用GDAL探探路。它支持超过200种栅格格式，几乎能搞定所有遥感数据。

下面这段代码，是我在项目中反复使用的“三板斧”：

from osgeo import gdal
import numpy as np

# 第一步：打开文件
ds = gdal.Open('landsat8_multispectral.tif')
if ds is None:
    print("文件打不开，检查路径或格式")
    exit()

# 第二步：获取基本信息
cols = ds.RasterXSize      # 列数（宽度）
rows = ds.RasterYSize      # 行数（高度）
bands = ds.RasterCount     # 波段数
print(f"影像尺寸：{cols} x {rows}，波段数：{bands}")

# 第三步：读取某个波段（比如第3波段）
band3 = ds.GetRasterBand(3)
data = band3.ReadAsArray()
print(f"波段3的统计：最小值={data.min()}，最大值={data.max()}")

# 第四步：获取地理变换参数
geotransform = ds.GetGeoTransform()
print(f"左上角坐标：({geotransform[0]}, {geotransform[3]})")
print(f"像元分辨率：{geotransform[1]} x {geotransform[5]}")

# 记得关闭数据集
ds = None

3.3 用Rasterio读取——更Pythonic，更优雅

Rasterio是GDAL的Python封装，但用起来舒服多了。它遵循Python的上下文管理器（with语句），代码更简洁，也更不容易出错。我个人现在更倾向于用Rasterio做快速原型开发。

import rasterio
import matplotlib.pyplot as plt

# 用with语句自动管理资源
with rasterio.open('sentinel2_multispectral.tif') as src:
    # 读取元数据
    print(f"CRS: {src.crs}")
    print(f"边界: {src.bounds}")
    print(f"波段索引: {src.indexes}")
    
    # 读取所有波段（返回三维数组：波段×行×列）
    data = src.read()
    print(f"数据形状: {data.shape}")
    
    # 读取单个波段
    red_band = src.read(4)  # Sentinel-2的第4波段是红波段
    nir_band = src.read(8)  # 第8波段是近红外
    
    # 计算NDVI（归一化植被指数）
    ndvi = (nir_band.astype(float) - red_band.astype(float)) / \
           (nir_band + red_band + 1e-10)  # 加小值防除零

# 可视化NDVI
plt.figure(figsize=(10, 8))
plt.imshow(ndvi, cmap='RdYlGn', vmin=-0.2, vmax=0.8)
plt.colorbar(label='NDVI')
plt.title('植被指数 - 基于Sentinel-2数据')
plt.show()

3.4 三种格式的读取对比

为了让你更直观地理解，我整理了一个对比表。嗯，这里要注意，不同格式的读取方式差异很大，但核心逻辑是一样的：

格式	GDAL打开方式	Rasterio打开方式	常见坑点
GeoTIFF	gdal.Open('file.tif')	rasterio.open('file.tif')	压缩方式（LZW、DEFLATE）可能导致读取慢
HDF	gdal.Open('HDF4_EOS:EOS_GRID:"file.hdf":Grid:band')	需先转GeoTIFF或使用xarray	子数据集路径复杂，容易写错
NetCDF	gdal.Open('NETCDF:"file.nc":variable')	rasterio.open('file.nc', variable='sst')	维度顺序（time, lat, lon）与预期不符

说白了，GeoTIFF是最省心的，HDF和NetCDF需要多花点心思。我建议初学者先从GeoTIFF入手，等熟悉了再挑战后两种。

3.5 可视化技巧——让数据“说话”

数据读出来了，怎么展示？多光谱数据不像普通RGB照片，它有多个波段，需要合理组合才能看到有用的信息。

我常用的可视化方法有三种：

1. 真彩色合成：用红、绿、蓝波段（比如Landsat的4、3、2波段）合成，看起来像真实照片。
2. 假彩色合成：用近红外、红、绿波段（比如5、4、3），植被会显示为红色，特别适合分析植被覆盖。
3. 单波段灰度图：直接显示某个波段的灰度值，适合分析特定地物的反射特征。

import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt

with rasterio.open('landsat8_multispectral.tif') as src:
    # 读取三个波段用于假彩色合成
    red = src.read(4)    # 红波段
    green = src.read(3)  # 绿波段
    blue = src.read(2)   # 蓝波段
    nir = src.read(5)    # 近红外波段
    
    # 创建子图
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 真彩色
    rgb = np.stack([red, green, blue], axis=0)
    # 需要做拉伸增强
    rgb_stretched = np.clip((rgb / rgb.max()) * 255, 0, 255).astype(np.uint8)
    axes[0, 0].imshow(np.transpose(rgb_stretched, (1, 2, 0)))
    axes[0, 0].set_title('真彩色合成')
    
    # 假彩色（近红外、红、绿）
    false_color = np.stack([nir, red, green], axis=0)
    false_stretched = np.clip((false_color / false_color.max()) * 255, 0, 255).astype(np.uint8)
    axes[0, 1].imshow(np.transpose(false_stretched, (1, 2, 0)))
    axes[0, 1].set_title('假彩色合成（NIR-R-G）')
    
    # 单波段显示
    axes[1, 0].imshow(red, cmap='gray')
    axes[1, 0].set_title('红波段（灰度）')
    
    axes[1, 1].imshow(nir, cmap='gray')
    axes[1, 1].set_title('近红外波段（灰度）')
    
    plt.tight_layout()
    plt.show()

3.6 知识体系图

下面这张图，是我梳理的本章知识结构。你可以把它当作一个“导航图”，遇到问题就知道该查哪部分。

3.7 避坑指南——我踩过的那些坑

最后，分享几个我亲身经历过的坑，希望能帮你少走弯路：

• 坐标参考系统（CRS）不一致：有一次我把Landsat数据和Sentinel-2数据叠加，结果错位了十几公里。后来发现一个是WGS84，一个是UTM投影。解决办法是用 rasterio.warp.reproject() 统一CRS。
• 内存爆炸：读取高分辨率多光谱数据时，直接 src.read() 可能会把内存撑爆。我建议用 src.read(window=window) 分块读取，或者用 src.read(1, out_shape=(rows//10, cols//10)) 降采样。
• 数据范围问题：多光谱数据通常是16位整型（0-65535），但有些数据是浮点型（比如反射率0-1）。可视化时如果不做归一化，显示出来就是一片白或一片黑。

---

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