Rasterio读取遥感影像?波段计算怎么搞?

wowwwai GIS研习社 · 工具流程与项目排障

Rasterio读取遥感影像?波段计算怎么搞? 这个问题很适合用一个完整流程来理解：先打开遥感影像，确认坐标系、尺寸、波段数量和 NoData 值，再按数组方式读取波段，最后完成 NDVI、波段差值、比例指数或自定义栅格计算，并把结果正确写回 GeoTIFF。

引言：Rasterio读取遥感影像与波段计算的基本思路

在 Python GIS 工作流里，Rasterio读取遥感影像常用于处理 GeoTIFF、遥感分类结果、DEM、单波段指数图和多波段卫星影像。它的优势是直接面向栅格数据，既能读取空间参考信息，也能把像元值转成 NumPy 数组进行计算。

很多初学者卡住的地方不是代码语法，而是这几个问题：

• 不知道 Rasterio 的波段编号为什么从 1 开始。
• 读取出来的数组为什么没有地理坐标。
• 做 NDVI 等波段计算时，为什么结果全是异常值。
• 写出新影像后，为什么在 QGIS 或 ArcGIS Pro 里位置不对。
• NoData、数据类型、除零和坐标变换没有处理好。

下面以 GeoTIFF 遥感影像为例，演示一个可复用的 Rasterio 波段读取和计算模板。

背景：为什么遥感影像不能只当普通图片读取

遥感影像虽然看起来像图片，但它和普通 PNG、JPG 有明显区别。GeoTIFF 这类 GIS 栅格数据不仅保存像元值，还保存坐标系、仿射变换、分辨率、范围、NoData、波段数量和数据类型。

如果只用普通图像库读取影像，可能会丢失这些空间信息。结果就是：计算值看起来没问题，但输出影像无法在地图中正确叠加。

Rasterio 的作用就是在读取像元数组的同时保留 GIS 栅格元数据。典型遥感影像波段计算包括：

• NDVI：利用近红外和红光波段计算植被指数。
• NDWI：利用绿光和近红外或短波红外计算水体指数。
• NBR：利用近红外和短波红外计算火烧迹地指数。
• 波段差值：用于变化检测或简单增强。
• 波段比例：用于突出某类地物光谱差异。

原理：Rasterio读取遥感影像时到底读到了什么

Rasterio 打开一个栅格文件后，主要可以获取两类信息：一类是元数据，一类是像元数组。

元数据通常包括：

• crs：坐标参考系统，例如 EPSG:4326 或 EPSG:32650。
• transform：仿射变换参数，用于把行列号转换为空间坐标。
• width 和 height：影像列数和行数。
• count：波段数量。
• dtype：像元数据类型，例如 uint16、int16、float32。
• nodata：无效值标记。
• profile：写出新栅格时最常用的完整参数字典。

像元数组则是通过 src.read(band_index) 读取的 NumPy 数组。需要注意，Rasterio 的波段编号从 1 开始，不是从 0 开始。例如第 1 波段写作 src.read(1)。

理解这一点很重要：Rasterio 负责地理空间信息，NumPy 负责数值计算。波段计算本质上就是在保持栅格空间元数据一致的前提下，对数组进行数学运算。

步骤：用 Rasterio 读取遥感影像并查看基本信息

先安装依赖。建议在独立 Python 环境中安装 Rasterio，避免和已有 GDAL 环境冲突。

pip install rasterio numpy

下面代码用于打开一个 GeoTIFF 并检查基础信息：

import rasterio

input_tif = "input_multiband.tif"

with rasterio.open(input_tif) as src:
    print("文件路径:", src.name)
    print("宽度:", src.width)
    print("高度:", src.height)
    print("波段数:", src.count)
    print("坐标系:", src.crs)
    print("仿射变换:", src.transform)
    print("数据类型:", src.dtypes)
    print("NoData:", src.nodata)
    print("范围:", src.bounds)
    print("Profile:", src.profile)

如果这里发现 src.crs 是 None，说明影像没有写入坐标系。后续即使波段计算成功，输出结果也可能无法正确叠加到底图上。

步骤：Rasterio读取单个波段和多个波段

读取单个波段时，使用 src.read(1)。返回结果是二维数组，形状通常为 (height, width)。

import rasterio

with rasterio.open("input_multiband.tif") as src:
    band1 = src.read(1)
    print(band1.shape)
    print(band1.dtype)
    print(band1.min(), band1.max())

读取多个波段时，可以分别读取，也可以一次读取全部波段。

import rasterio

with rasterio.open("input_multiband.tif") as src:
    all_bands = src.read()
    print(all_bands.shape)

src.read() 读取全部波段时，数组形状通常是 (band_count, height, width)。这和很多图像库常见的 (height, width, channel) 不一样，写代码时要特别注意。

步骤：用 Rasterio 做 NDVI 波段计算

NDVI 是最常见的遥感波段计算之一，公式为：

NDVI = (NIR - Red) / (NIR + Red)

不同卫星影像的红光和近红外波段编号不同。以常见多光谱数据为例，你必须先确认影像的波段说明，而不是盲目套用代码。

下面是一个完整的 Rasterio NDVI 计算示例。假设红光是第 4 波段，近红外是第 5 波段。

import numpy as np
import rasterio

input_tif = "input_multiband.tif"
output_tif = "output_ndvi.tif"

red_band_index = 4
nir_band_index = 5

with rasterio.open(input_tif) as src:
    red = src.read(red_band_index).astype("float32")
    nir = src.read(nir_band_index).astype("float32")

    nodata = src.nodata
    profile = src.profile.copy()

    denominator = nir + red

    ndvi = np.where(
        denominator == 0,
        np.nan,
        (nir - red) / denominator
    )

    if nodata is not None:
        mask = (red == nodata) | (nir == nodata)
        ndvi[mask] = np.nan

    profile.update(
        count=1,
        dtype="float32",
        nodata=np.nan,
        compress="lzw"
    )

    with rasterio.open(output_tif, "w", **profile) as dst:
        dst.write(ndvi.astype("float32"), 1)

这段代码做了几个关键处理：

• 把红光和近红外转为 float32，避免整数除法和数据类型溢出。
• 对分母为 0 的位置写入 np.nan，避免无穷值。
• 把原始 NoData 区域排除在计算之外。
• 更新输出 profile，确保输出是单波段 float32 GeoTIFF。
• 保留原始影像的坐标系、范围、分辨率和仿射变换。

步骤：Rasterio波段计算后如何写出 GeoTIFF

写出结果时，最稳妥的方法是复制原始影像的 profile，然后只修改必要字段。这样可以最大程度保留原始空间参考。

profile = src.profile.copy()
profile.update(
    count=1,
    dtype="float32",
    nodata=np.nan,
    compress="lzw"
)

常见需要修改的字段包括：

• count：输出波段数。NDVI 通常是 1。
• dtype：输出数据类型。指数计算通常建议使用 float32。
• nodata：输出无效值。可以使用 np.nan 或指定数值。
• compress：压缩方式，例如 lzw，可减小文件体积。

如果输出影像在 QGIS 中打开后显示为纯黑或纯白，不一定是计算错了。指数值通常在 -1 到 1 之间，需要在符号系统中调整拉伸方式或设置合适的色带。

常见坑：Rasterio读取遥感影像波段计算为什么出错

1. 波段编号写错

Rasterio 的波段从 1 开始。如果写 src.read(0) 会报错。如果把红光和近红外波段弄反，NDVI 结果可能整体偏负。

2. 没有转成 float32

遥感原始波段常见类型是 uint16 或 int16。做比例计算时建议先转为 float32，否则容易出现精度问题或异常结果。

3. 没处理 NoData

NoData 像元参与计算会污染结果。例如 NoData 值为 0、-9999 或 65535 时，如果直接计算指数，可能得到大量异常值。

4. 分母为 0

NDVI、NDWI 等指数都有除法。如果分母为 0，需要用 np.where 或掩膜处理，否则结果中会出现 inf 或 nan。

5. Sentinel-2 不同波段分辨率不同

Sentinel-2 的部分波段是 10 米、20 米或 60 米分辨率。如果直接把不同尺寸的波段数组相减，会出现数组形状不一致。需要先重采样到同一分辨率。

6. 输出 profile 没更新

如果输入是多波段影像，输出 NDVI 却忘记设置 count=1，写出时容易出错。输出数据类型也要和计算结果匹配。

方法比较：Rasterio、GDAL、QGIS字段计算器该怎么选

如果只是临时计算一幅 NDVI，QGIS 栅格计算器很方便。如果要处理几十景或几百景影像，Rasterio读取遥感影像再结合 NumPy 做波段计算更适合自动化。

检查清单：运行 Rasterio 波段计算前后要确认什么

• 确认输入文件是否为 GeoTIFF 或 Rasterio 支持的栅格格式。
• 确认 src.crs 不为空，空间参考正确。
• 确认 src.count 满足计算所需波段数量。
• 确认红光、近红外、绿光、短波红外等波段编号没有写错。
• 确认参与计算的波段数组尺寸一致。
• 确认原始 NoData 值，并在计算中排除。
• 确认除法公式中分母为 0 的像元已处理。
• 确认输出 profile 的 count、dtype、nodata 设置正确。
• 确认输出 GeoTIFF 能在 QGIS 或 ArcGIS Pro 中正确叠加。
• 确认可视化拉伸方式正确，不要把显示异常误判为计算错误。

FAQ：Rasterio读取遥感影像波段计算常见问题

Rasterio读取遥感影像时，为什么 read(1) 不是第一行而是第一波段？

因为 src.read(1) 的参数表示波段编号，不是行号。Rasterio 的波段编号从 1 开始，读取后得到的二维数组才包含行和列。

Rasterio波段计算可以直接处理 JPG 或 PNG 吗？

技术上可以读取部分普通图像格式，但普通 JPG 或 PNG 通常没有完整 GIS 空间参考。做遥感分析时，建议使用 GeoTIFF 或其他带空间信息的栅格格式。

为什么 NDVI 结果不在 -1 到 1 之间？

常见原因包括波段选错、NoData 没处理、分母为 0、原始数据需要先做缩放或反射率转换。不同数据产品的像元值含义不同，计算前应查看数据说明。

输出的 NDVI 在 QGIS 里一片黑，是不是算错了？

不一定。NDVI 是 float32 指数值，范围通常较窄。需要在 QGIS 中调整渲染方式，例如使用单波段伪彩色，设置最小最大值为 -1 到 1，或使用累计计数裁剪拉伸。

多个波段尺寸不一样，Rasterio能直接计算吗？

不能直接计算。NumPy 数组形状不一致时无法直接相加相减。需要先用重采样方法把不同波段统一到同一分辨率、同一范围和同一行列数。

Rasterio读取大影像时内存不够怎么办？

不要一次性读取整幅影像。可以使用 Rasterio 的窗口读取方式按块处理，也可以先裁剪研究区，或使用 Cloud Optimized GeoTIFF、Dask、xarray 等更适合大数据的方案。

结论：把读取、计算和写出连成一个可靠流程

Rasterio读取遥感影像做波段计算，核心不是简单写一行公式，而是保证波段编号、数据类型、NoData、空间参考和输出 profile 都正确。只要把这些环节检查清楚，NDVI、NDWI、波段差值和比例指数都可以用同一套模板稳定完成。

实际项目中建议先用一景小范围影像测试流程，在 QGIS 或 ArcGIS Pro 中验证输出位置和数值范围，再扩展到批量处理。这样既能减少错误，也能让 Python GIS 遥感处理流程更容易复用。