3、直方图可视化:使用 Matplotlib 绘制单波段与多波段直方图

直方图这东西,说白了就是一张“像素值分布地图”。你拿到一张遥感影像,第一件事往往不是直接做增强,而是先看看它的直方图长什么样。我个人的习惯是:先看图,再动手。直方图能告诉你影像的亮度范围、对比度高低、有没有过曝或欠曝,甚至能帮你判断地物类型的分布情况。

这一章我们就来实战——用 Python 的 Matplotlib 库,把单波段和多波段的直方图画出来。嗯,代码我会给全,但更重要的是理解每张图背后的含义。

3.1 单波段直方图:从一维数组到可视化

单波段直方图是最基础的形式。你有一个波段的数据(比如近红外波段),把它拉成一维数组,然后统计每个灰度值出现的次数,再画成柱状图——就这么简单。

但这里有个坑:遥感影像通常是 16 位或 32 位浮点型,不像普通照片是 0-255 的 8 位整数。如果你直接拿原始数据去画直方图,横坐标会变得非常稀疏,看起来像一堆刺。我刚开始做的时候就被这个坑过,画出来的图根本没法看。

我的经验:画直方图前,先把数据缩放到 0-255 范围,或者用 Matplotlib 的 bins 参数控制分组数量。对于 16 位影像,bins=256 是个不错的起点。

来看一段核心代码:

import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal

# 读取单波段数据
ds = gdal.Open('landsat_b4.tif')
band = ds.GetRasterBand(1).ReadAsArray()

# 展平为一维数组
data_flat = band.flatten()

# 过滤掉无效值(比如 NoData 值)
data_valid = data_flat[data_flat > 0]

# 绘制直方图
plt.figure(figsize=(10, 6))
plt.hist(data_valid, bins=256, color='steelblue', edgecolor='none', alpha=0.8)
plt.title('单波段直方图 - Landsat 8 Band 4 (近红外)')
plt.xlabel('像素值')
plt.ylabel('频数')
plt.grid(axis='y', alpha=0.3)
plt.show()

这段代码里,我特别加了一行过滤无效值的操作。为什么?因为遥感影像的边缘通常有黑色填充(值为 0),或者有 NoData 标记值。如果不过滤,这些无效像素会严重干扰直方图的形态,让你误以为影像整体偏暗。

注意:不同传感器的 NoData 值不一样。Landsat 常用 0 或 -9999,Sentinel-2 可能是 NaN。画图前一定要确认你的数据里有没有特殊值。

3.2 多波段直方图:叠加与对比

单波段直方图只能看一个波段的信息。但遥感影像通常是多波段的——比如 Landsat 8 有 11 个波段。把多个波段的直方图画在一起,能直观地对比不同波段的光谱响应范围。

我记得有一次做植被覆盖度分析,需要对比红波段和近红外波段的分布差异。把两个直方图叠在一起画,一眼就能看出近红外波段的像素值整体偏高——这就是植被在近红外波段强反射的典型特征。

多波段直方图有两种常见画法:

画法 适用场景 优点 缺点
叠加直方图 波段数量少(2-4 个) 对比直观,一眼看出差异 波段多了会重叠看不清
子图排列 波段数量多(5 个以上) 每个波段独立显示,清晰 对比不够直接

叠加直方图的代码实现:

# 读取三个波段
bands_data = []
band_names = ['蓝波段', '绿波段', '红波段']
colors = ['blue', 'green', 'red']

for i in [1, 2, 3]:  # Landsat 8 前三个波段
    band = ds.GetRasterBand(i).ReadAsArray()
    data_valid = band[band > 0].flatten()
    bands_data.append(data_valid)

# 叠加绘制
plt.figure(figsize=(12, 6))
for data, name, color in zip(bands_data, band_names, colors):
    plt.hist(data, bins=256, alpha=0.5, label=name, color=color, edgecolor='none')

plt.title('多波段直方图叠加 - Landsat 8 RGB')
plt.xlabel('像素值')
plt.ylabel('频数')
plt.legend()
plt.grid(axis='y', alpha=0.3)
plt.show()

这里 alpha=0.5 是关键——透明度让重叠区域能透出下面的颜色。如果不设透明度,后面的波段会把前面的完全盖住,那就白画了。

3.3 直方图解读:从图形到信息

画图不是目的,读懂图才是。我总结了几个直方图形态对应的影像特征:

  • 直方图集中在左侧:影像整体偏暗,可能是阴影区域多,或者传感器增益设置偏低。
  • 直方图集中在右侧:影像整体偏亮,可能是雪地、沙漠等高反射地物占比大。
  • 直方图呈双峰或多峰:影像中存在两种或多种明显不同的地物类型。比如一个峰对应水体(低值),另一个峰对应植被(中高值)。
  • 直方图很窄:影像对比度低,地物之间的光谱差异小,需要做拉伸增强。
  • 直方图有断崖或空洞:可能是数据压缩或量化造成的,也可能是传感器故障。
实战案例:有一次我处理某地区的多光谱影像,发现红波段的直方图在值 100 附近有个明显的“缺口”。排查后发现是传感器在采集时出现了条带噪声。直方图帮我快速定位了问题,省去了逐像素检查的麻烦。

3.4 直方图可视化流程框架

下面这张图总结了从原始影像到直方图输出的完整流程,以及不同场景下的选择路径:

直方图可视化流程框架 原始遥感影像 数据预处理(去无效值、缩放) 波段数量判断 单波段 多波段 单波段直方图绘制 多波段直方图绘制 叠加直方图 子图排列 直方图输出与解读 关键步骤 去无效值 → 缩放 → 选择画法 → 输出

3.5 常见问题与避坑指南

画直方图看起来简单,但实际操作中会遇到各种问题。我把自己踩过的坑整理了一下:

  • 直方图全是尖刺,看不出分布:bins 数量设置太多。对于 16 位影像,bins=256 通常够用;对于 8 位影像,bins=64 或 128 更合适。
  • 直方图两端有异常高峰:大概率是 NoData 值或边界填充值没过滤。用 np.percentile() 先看看数据的有效范围。
  • 多波段叠加时颜色混在一起看不清:试试用不同的线型(histtype='step')代替填充柱状图,或者改用子图排列。
  • 直方图横坐标范围太大:用 plt.xlim() 手动设置显示范围,聚焦在有效数据区间。
我曾经踩过的坑:有一次画 Sentinel-2 的直方图,发现所有波段的直方图都集中在 0-100 之间,我以为数据有问题。后来才发现 Sentinel-2 的原始数据是除以了 10000 的反射率值,范围本来就是 0-1 左右。所以画图前一定要搞清楚数据的物理含义和量纲。