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

NDVI一片黑？别慌，可能是你波段顺序搞反了

上周一个研究生私信我：‘老师，我用Python算出来的NDVI全黑，拉伸也没用，是不是代码写错了？’——结果一查，红光和近红外波段位置颠倒了。这问题太常见了！今天我就手把手带你从原理到实战，彻底搞懂NDVI计算中那些‘数组操作’的坑。

NDVI本质：一场植被与光的‘恋爱指数’

NDVI（归一化植被指数）公式很简单：(NIR - Red) / (NIR + Red)。但为什么是‘近红外减红光’？因为健康植被会强烈反射近红外光（像在说‘我爱你’），同时吸收红光（像在说‘别靠近我’）。这个比值越高，说明植被越茂盛——就像恋爱中一方热情似火、一方保持距离，反而感情稳定（笑）。

Dr. Gis经验谈：我在做黄土高原生态修复评估时发现，农田休耕期NDVI骤降不是算法问题，而是作物真被收割了！别总怀疑代码，先看数据合理性。

波段数组操作：把卫星影像当‘乐高积木’拆解

遥感影像本质是三维数组：行×列×波段。假设你用Landsat 8数据：

• 波段4 = 红光（Red）
• 波段5 = 近红外（NIR）

关键来了！Python读取时默认按波段顺序存储，如果你用band[3]取红光（第4波段），band[4]取近红外（第5波段）——恭喜，你成功制造了一片黑图（负值溢出）。正确姿势：

import numpy as np
from osgeo import gdal

# 读取影像
dataset = gdal.Open('LC08_L1TP.tif')
red_band = dataset.GetRasterBand(4).ReadAsArray()  # 波段4=红光
nir_band = dataset.GetRasterBand(5).ReadAsArray()  # 波段5=近红外

# 计算NDVI（重点：避免除零！）
ndvi = (nir_band.astype(float) - red_band.astype(float)) / 
       (nir_band.astype(float) + red_band.astype(float) + 1e-8)

这里加1e-8是防止分母为零——就像给恋爱双方留个台阶，避免‘除零崩溃’。

避坑指南：三个必须检查的致命细节

进阶技巧：用GeoPandas批量处理多景影像

当你有100景影像要处理时，手动改路径会疯掉。试试这个自动化脚本：

import glob
from osgeo import gdal

for img_path in glob.glob('*.tif'):
    ds = gdal.Open(img_path)
    red = ds.GetRasterBand(4).ReadAsArray().astype(float)
    nir = ds.GetRasterBand(5).ReadAsArray().astype(float)
    ndvi = (nir - red) / (nir + red + 1e-8)
    
    # 输出结果（略去驱动创建步骤）
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(f'NDVI_{img_path}', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32)
    out_ds.GetRasterBand(1).WriteArray(ndvi)
    out_ds.SetProjection(ds.GetProjection())
    out_ds.SetGeoTransform(ds.GetGeoTransform())
    out_ds = None  # 关闭文件

总结：NDVI计算就三步走

1. 确认波段：打印元数据核对红光/NIR位置
2. 转浮点数：避免整数运算截断
3. 防除零：分母加极小值保平安

现在轮到你了！你在计算NDVI时踩过什么坑？是在波段顺序、数据类型还是投影上翻车？评论区留下你的血泪史，点赞最高的送《遥感指数计算避坑手册》PDF！