Python计算NDVI公式?波段数组咋操作?

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

很多同学在做遥感植被指数时会直接搜索“Python计算NDVI公式?波段数组咋操作?”，真正卡住的往往不是 NDVI 公式本身，而是红光波段、近红外波段读成数组以后，如何处理数据类型、NoData、除零和结果保存。

引言：Python计算NDVI公式到底怎么落到数组上

NDVI，全称归一化植被指数，常用于快速判断植被覆盖和长势。它的核心公式很简单：

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

其中 NIR 是近红外波段，Red 是红光波段。用 Python 计算 NDVI 时，关键就是把这两个波段读取为 NumPy 数组，然后对数组逐像元套用这个公式。

本文以 GeoTIFF 影像为例，使用 Rasterio + NumPy 演示完整流程：读取波段、数组计算、处理无效值、保存 NDVI 结果，并解释为什么很多人算出来会全是 0、全是空值，或者结果范围不对。

背景：为什么Python计算NDVI经常不是公式问题

在 GIS 和遥感处理中，NDVI 公式本身并不难，真正容易出错的是影像波段和数组细节。

• 不同卫星影像的红光波段和近红外波段编号不同。
• 原始像元值可能是整数反射率，需要转换为浮点数。
• 影像中可能存在 NoData 区域，直接计算会污染结果。
• NIR + Red 等于 0 时会出现除零问题。
• 如果使用整数数组计算，结果可能被截断或异常。
• 输出 GeoTIFF 时，数据类型、坐标系和仿射变换需要继承原始影像。

因此，Python计算NDVI公式的正确姿势不是只写一行公式，而是建立一套完整的波段数组操作流程。

原理：NDVI公式、波段数组和结果范围

NDVI公式含义

NDVI 利用绿色植被在近红外波段反射强、在红光波段吸收强的特征来表达植被信息。公式如下：

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

理论上，NDVI 的结果范围通常在 -1 到 1 之间：

• 接近 1：植被较旺盛。
• 接近 0：裸地、建筑物、低植被覆盖等。
• 小于 0：水体、阴影、云等区域较常见。

波段数组咋操作

在 Python 中，Rasterio 读取单个波段后，返回的是二维 NumPy 数组。假设红光波段数组为 red，近红外波段数组为 nir，那么 NDVI 就是两个同尺寸数组之间的逐像元运算。

ndvi = (nir - red) / (nir + red)

这行代码看起来简单，但必须保证 red 和 nir 是浮点数组，并且已经处理了无效值和除零情况。

步骤：用Rasterio和NumPy完整计算NDVI

步骤1：确认红光和近红外波段编号

不同数据源的波段顺序不一样，不能盲目套用固定编号。常见情况如下：

如果你不确定波段编号，可以先查看影像元数据，或者在 QGIS 中打开图层属性，检查波段名称和显示组合。

步骤2：安装需要的Python库

推荐使用 Rasterio 读取和写入 GeoTIFF，使用 NumPy 进行数组计算。

pip install rasterio numpy

如果你使用的是 Conda 环境，也可以使用：

conda install -c conda-forge rasterio numpy

步骤3：读取红光和近红外波段数组

下面示例假设输入影像是一个多波段 GeoTIFF，并且红光波段是第 4 波段，近红外波段是第 5 波段。这个编号适合部分 Landsat 8/9 数据，但你需要根据自己的数据调整。

import rasterio
import numpy as np

input_tif = "input_multiband.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")
    profile = src.profile
    nodata = src.nodata

这里的 astype("float32") 很重要。NDVI 是小数结果，如果继续使用整数数组，容易造成计算精度问题。

步骤4：处理NoData和除零问题

计算 NDVI 前，建议先构建一个有效像元掩膜。常见无效条件包括：原始 NoData、红光和近红外同时为 0、分母为 0。

denominator = nir + red

valid_mask = denominator != 0

if nodata is not None:
    valid_mask = valid_mask & (red != nodata) & (nir != nodata)

ndvi = np.full(red.shape, -9999, dtype="float32")
ndvi[valid_mask] = (nir[valid_mask] - red[valid_mask]) / denominator[valid_mask]

这段代码的思路是：先把输出 NDVI 数组全部填为 -9999，然后只对有效像元计算 NDVI。这样可以避免除零警告，也能防止无效区域参与统计。

步骤5：保存NDVI为GeoTIFF

输出 NDVI 栅格时，需要继承原始影像的坐标系、范围、分辨率和仿射变换，同时把波段数改成 1，把数据类型改成 float32。

output_tif = "ndvi_result.tif"

profile.update(
    dtype=rasterio.float32,
    count=1,
    nodata=-9999,
    compress="lzw"
)

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

保存后可以在 QGIS 或 ArcGIS Pro 中打开 ndvi_result.tif，设置单波段伪彩色渲染，检查结果是否符合预期。

完整代码：Python计算NDVI公式和波段数组操作

import rasterio
import numpy as np

input_tif = "input_multiband.tif"
output_tif = "ndvi_result.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")
    profile = src.profile
    nodata = src.nodata

denominator = nir + red
valid_mask = denominator != 0

if nodata is not None:
    valid_mask = valid_mask & (red != nodata) & (nir != nodata)

ndvi = np.full(red.shape, -9999, dtype="float32")
ndvi[valid_mask] = (nir[valid_mask] - red[valid_mask]) / denominator[valid_mask]

profile.update(
    dtype=rasterio.float32,
    count=1,
    nodata=-9999,
    compress="lzw"
)

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

print("NDVI计算完成：", output_tif)

常见坑：NDVI结果不对通常查这几项

1. 红光和近红外波段选反了

如果把 Red 和 NIR 选反，NDVI 值会整体偏负，植被区域反而变成低值。遇到这种情况，第一步应检查波段编号，而不是急着改公式。

2. 波段不是同一分辨率

Sentinel-2 的不同波段可能有 10 米、20 米、60 米等不同分辨率。NDVI 常用 Band 4 和 Band 8，它们通常同为 10 米；但如果你换用其他红边或短波红外波段，就可能需要先重采样。

3. 没有转成浮点数组

Python计算NDVI公式时，建议明确使用 float32 或 float64。这样可以保证数组除法得到小数结果，也方便后续保存连续型栅格。

4. 没有处理NoData

如果原始影像边缘、云区或裁剪外区域存在 NoData，直接参与计算会产生异常值。正确做法是使用掩膜，只对有效像元计算，并把结果 NoData 写入输出文件元数据。

5. 原始值是缩放后的整数反射率

有些遥感产品把反射率乘以比例因子后保存为整数。例如像元值可能需要乘以 0.0001 才是真实反射率。对于 NDVI 这种比值指数，如果 Red 和 NIR 使用同一比例因子，比例因子会在公式中抵消；但如果存在偏移量或不同缩放规则，就必须按产品说明先还原。

6. 输出结果看起来一片黑

这不一定是计算错了。NDVI 是 -1 到 1 左右的小数，如果软件按默认灰度拉伸显示，可能不直观。建议在 QGIS 中使用单波段伪彩色，并设置合理的最小值和最大值，例如 -0.2 到 0.8。

方法比较：Rasterio、GDAL、QGIS谁更适合计算NDVI

如果你正在学习 Python GIS，推荐从 Rasterio + NumPy 开始。它能让你真正理解“波段数组咋操作”，以后扩展到 EVI、NDBI、NDWI 等指数也很自然。

检查清单：运行前后逐项确认

• 确认输入影像是多波段 GeoTIFF，或者能分别读取 Red 和 NIR 两个单波段文件。
• 确认红光波段和近红外波段编号没有选错。
• 确认两个波段尺寸、坐标系、分辨率一致。
• 读取数组后使用 astype("float32") 转为浮点类型。
• 计算前处理 NoData 和 NIR + Red = 0 的像元。
• 输出 NDVI 时设置 count=1、dtype=float32、nodata=-9999。
• 在 QGIS 或 ArcGIS Pro 中检查 NDVI 统计范围是否大致在 -1 到 1。
• 使用伪彩色渲染查看植被区域是否表现为高值。

FAQ：Python计算NDVI公式常见问题

Python计算NDVI公式一定要用Rasterio吗？

不一定。也可以使用 GDAL、rioxarray、ArcPy 或 QGIS 处理工具。但 Rasterio 对 GeoTIFF 支持好，和 NumPy 数组结合自然，非常适合学习和编写批处理脚本。

NDVI的红光和近红外波段怎么确定？

应以数据产品说明为准。Landsat 8/9 常用 Band 4 作为 Red，Band 5 作为 NIR；Sentinel-2 常用 Band 4 作为 Red，Band 8 作为 NIR。无人机或自定义多光谱影像必须查看传感器和导出说明。

波段数组操作时为什么要转float32？

NDVI 是小数指数，使用浮点数组可以避免整数计算带来的精度问题。输出为 float32 通常已经足够，也能控制文件体积。

NDVI结果有大于1或小于-1的值正常吗？

少量异常值可能来自 NoData、云、阴影、传感器异常或缩放规则处理不当。建议先检查无效值掩膜、波段编号和原始数据范围。如果大面积超出 -1 到 1，通常说明输入波段或预处理存在问题。

两个单波段文件能不能计算NDVI？

可以，但需要确保两个文件范围、分辨率、行列数和坐标系一致。读取方式可以分别打开 Red 和 NIR 文件，再用同样的数组公式计算。如果不一致，需要先重投影、重采样或对齐栅格。

Python计算NDVI后如何在QGIS中显示？

把输出的 ndvi_result.tif 拖入 QGIS，图层样式选择“单波段伪彩色”，设置合适的颜色带，并手动设置最小值和最大值，例如 -0.2 和 0.8。这样比默认灰度显示更容易判断植被分布。

结论：把公式、数组和栅格元数据一起处理才算完整

“Python计算NDVI公式?波段数组咋操作?”这个问题的核心答案是：先确认 Red 和 NIR 波段，再用 Rasterio 读取为浮点 NumPy 数组，按 (NIR - Red) / (NIR + Red) 逐像元计算，同时处理 NoData 和除零，最后按原始影像的空间参考保存 GeoTIFF。

对于 GIS 初学者来说，NDVI 是理解遥感栅格数组运算的很好入口。掌握这套流程后，你可以用同样思路继续计算 NDWI、NDBI、EVI 等指数，也可以把脚本扩展为批量处理多个影像的自动化工具。