4. NumPy与遥感数据矩阵运算：数组基础、矩阵切片、统计计算、掩膜处理

说实话，刚接触遥感数据处理那会儿，我总觉得NumPy就是个数学工具库，跟遥感没啥关系。直到有一次处理Landsat影像，老板丢给我一个3000×3000的多波段矩阵，让我做植被指数计算。我愣是用Excel手动算了三天……后来同事告诉我，用NumPy三行代码就搞定了。嗯，从那以后，我再也不敢小看这个库了。

今天咱们就聊聊，怎么用NumPy这把“瑞士军刀”来切遥感数据这块“硬骨头”。

4.1 数组基础：遥感数据的“骨架”

遥感影像本质上就是个多维数组。你想啊，一张RGB影像，宽1000像素、高800像素，那就是一个800×1000×3的矩阵。每个像素点存着红、绿、蓝三个波段的反射率值。

我习惯把NumPy数组理解成“带标签的Excel表格”——只不过它更快、更省内存、而且能直接做数学运算。

import numpy as np

# 模拟一个100×100×4的多光谱影像（蓝、绿、红、近红外）
image = np.random.randint(0, 256, (100, 100, 4), dtype=np.uint8)
print(f"影像形状: {image.shape}")   # (100, 100, 4)
print(f"数据类型: {image.dtype}")    # uint8
print(f"总像素数: {image.size}")     # 40000

4.2 矩阵切片：像切蛋糕一样切影像

切片是NumPy最爽的功能之一。说白了，就是按行、列、波段去“抠”出你想要的子区域。

举个例子，我想从上面那个100×100的影像里，取出第20到40行、第30到50列、只保留红和近红外两个波段：

# 切片：行20:40，列30:50，波段2和3（红、近红外）
sub_image = image[20:40, 30:50, [2, 3]]
print(f"切片后形状: {sub_image.shape}")  # (20, 20, 2)

你想想看，如果不用切片，你得写多少层for循环？而且切片返回的是原数组的视图，不是拷贝——这意味着修改切片会直接影响原数据。我在项目中就踩过这个坑：

# 正确做法：需要独立副本时
sub_image_copy = image[20:40, 30:50, [2, 3]].copy()

4.3 统计计算：从像素到知识

遥感数据最常用的统计量就那几个：均值、标准差、最大值、最小值、百分位数。别小看这些，它们能告诉你很多信息。

比如，计算NDVI（归一化植被指数）之前，我习惯先看看近红外和红波段的统计分布：

nir = image[:, :, 3].astype(np.float32)  # 近红外
red = image[:, :, 2].astype(np.float32)  # 红波段

print(f"近红外: 均值={nir.mean():.2f}, 标准差={nir.std():.2f}")
print(f"红波段: 均值={red.mean():.2f}, 标准差={red.std():.2f}")

# 计算NDVI
ndvi = (nir - red) / (nir + red + 1e-10)  # 加小量防除零
print(f"NDVI范围: [{ndvi.min():.3f}, {ndvi.max():.3f}]")

另外，np.percentile 在遥感中特别有用。比如做影像拉伸显示时，我通常取2%和98%的百分位作为拉伸范围：

p2, p98 = np.percentile(nir, [2, 98])
print(f"2%分位值: {p2:.2f}, 98%分位值: {p98:.2f}")

4.4 掩膜处理：给数据“戴面具”

掩膜，说白了就是用一个布尔矩阵来筛选数据。在遥感里，最常见的场景就是去云、去水体、去无效值。

我记得有一次做土地利用分类，影像上有大片云层。如果不做掩膜，分类器会把云当成“白色地物”去学习，结果可想而知。

# 假设我们有一个云检测结果：1表示云，0表示晴空
cloud_mask = np.random.choice([0, 1], size=(100, 100), p=[0.8, 0.2])

# 创建一个布尔掩膜：True表示有效像素
valid_mask = cloud_mask == 0

# 应用掩膜：只保留有效像素
nir_valid = nir[valid_mask]
red_valid = red[valid_mask]

print(f"原始像素数: {nir.size}")
print(f"有效像素数: {nir_valid.size}")

掩膜还能做更复杂的操作。比如，我只想计算水体区域的NDVI均值：

# 假设有一个水体掩膜（1表示水体，0表示非水体）
water_mask = np.random.choice([0, 1], size=(100, 100), p=[0.7, 0.3])

# 计算水体区域的NDVI均值
ndvi_water = ndvi[water_mask == 1]
print(f"水体区域NDVI均值: {ndvi_water.mean():.3f}")

4.5 本章知识体系

下面这张图，是我自己总结的NumPy遥感数据处理流程。你看一眼，基本就明白今天讲了啥：

说白了，NumPy处理遥感数据的核心就四个字：向量化操作。别写for循环，别一个一个像素去算。用矩阵运算，一次搞定。

我个人觉得，掌握好这四块内容，你就能处理90%的遥感数据预处理工作了。剩下的10%，无非是更复杂的掩膜逻辑或者更精细的统计策略——但底层原理，还是今天讲的这些。

---

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