第四章：遥感影像读取与可视化

各位同学，今天我们来聊聊遥感影像的读取和可视化。说实话，这是整个农业金融遥感应用中最基础、也最容易被忽视的一环。我见过太多人拿着GeoTIFF文件，却不知道怎么把它变成能用的数据。

我个人习惯用Rasterio来处理遥感影像。为什么选它？因为它对GeoTIFF的支持特别好，而且跟NumPy、Matplotlib这些Python生态库配合得天衣无缝。你想想看，一个工具能搞定读取、处理、可视化，多省心。

4.1 用Rasterio读取GeoTIFF

先说说GeoTIFF。它本质上就是TIFF图片，但多了地理空间信息。比如经纬度、投影坐标系、像元大小等等。这些信息对农业金融分析至关重要——你得知道这块地到底在哪儿，面积多大。

读取代码其实很简单：

import rasterio

# 打开GeoTIFF文件
with rasterio.open('sentinel2_2024_08_15.tif') as src:
    # 读取所有波段
    data = src.read()
    # 获取元数据
    meta = src.meta
    # 获取地理变换信息
    transform = src.transform
    
print(f"影像尺寸: {src.width} x {src.height}")
print(f"波段数: {src.count}")
print(f"坐标系: {src.crs}")

嗯，这里要注意一点。我在项目中遇到过有人直接用src.read()不指定波段，结果内存爆了。为什么？因为高分辨率遥感影像动辄几十个波段，每个波段几千万像素，全读进来内存肯定扛不住。

4.2 波段合成

波段合成说白了就是把多个单波段影像叠在一起，变成多波段影像。这在农业金融里很常见——比如你要把不同时相的影像合成一个时间序列，或者把不同传感器的数据融合。

我给你们看个实际案例。去年做玉米种植面积核查时，我需要把红、绿、蓝、近红外四个波段合成一个四波段影像：

import numpy as np
import rasterio

# 假设我们已经读取了四个单波段
red = src.read(4)    # 红波段
green = src.read(3)  # 绿波段  
blue = src.read(2)   # 蓝波段
nir = src.read(8)    # 近红外波段

# 波段合成 - 堆叠成三维数组
composite = np.stack([red, green, blue, nir], axis=0)

# 保存为新的GeoTIFF
with rasterio.open(
    'composite_4band.tif', 'w',
    driver='GTiff',
    height=composite.shape[1],
    width=composite.shape[2],
    count=4,
    dtype=composite.dtype,
    crs=src.crs,
    transform=src.transform
) as dst:
    dst.write(composite)

4.3 使用Matplotlib展示影像

影像读进来了，怎么展示？Matplotlib是最直接的选择。但遥感影像跟普通照片不一样，它的像素值范围可能很大（比如16位整型，0-65535），直接显示会一片漆黑。

我一般会做拉伸处理：

import matplotlib.pyplot as plt
import numpy as np

def stretch_histogram(band, lower_percent=2, upper_percent=98):
    """直方图拉伸，增强显示效果"""
    band = band.astype(np.float32)
    lower = np.percentile(band, lower_percent)
    upper = np.percentile(band, upper_percent)
    band = np.clip(band, lower, upper)
    band = (band - lower) / (upper - lower)
    return band

# 读取RGB波段
with rasterio.open('sentinel2_2024_08_15.tif') as src:
    r = stretch_histogram(src.read(4))
    g = stretch_histogram(src.read(3))
    b = stretch_histogram(src.read(2))

# 合成RGB影像
rgb = np.stack([r, g, b], axis=-1)

# 显示
plt.figure(figsize=(12, 8))
plt.imshow(rgb)
plt.title('真彩色合成 - 2024年8月15日')
plt.axis('off')
plt.show()

你看，这样显示出来的影像就清晰多了。为什么用2%和98%的百分位？说白了就是去掉两头的极端值，让中间的主体部分显示得更清楚。这个技巧我在做农田长势监测时经常用。

4.4 真彩色与假彩色合成

真彩色合成就是用红、绿、蓝波段模拟人眼看到的颜色。但遥感影像的价值远不止于此。假彩色合成能揭示肉眼看不到的信息。

在农业金融里，最常用的是近红外-红-绿假彩色合成。为什么？因为植被在近红外波段反射率特别高，显示出来是红色。红色越亮，说明植被越茂盛。这对判断作物长势、估算产量太有用了。

# 假彩色合成：近红外(R) - 红(G) - 绿(B)
with rasterio.open('sentinel2_2024_08_15.tif') as src:
    nir = stretch_histogram(src.read(8))  # 近红外 -> 红色通道
    red = stretch_histogram(src.read(4))  # 红 -> 绿色通道
    green = stretch_histogram(src.read(3)) # 绿 -> 蓝色通道

# 合成假彩色影像
false_color = np.stack([nir, red, green], axis=-1)

plt.figure(figsize=(12, 8))
plt.imshow(false_color)
plt.title('假彩色合成 (NIR-R-G) - 2024年8月15日')
plt.axis('off')
plt.show()

这里有个表格，总结了不同合成方式的用途：

嗯，最后说一句。可视化不只是为了好看，它是你理解数据的第一步。我每次拿到新的遥感数据，第一件事就是做各种合成看看。有时候肉眼发现的问题，比跑一堆算法来得更快。