1 文章背景在遥感影像处理中归一化植被指数 NDVI 是一个非常常见的植被指数。它可以根据近红外波段和红光波段之间的反射特征反映地表植被覆盖和植被生长状况。本文主要整理 NDVI 的基本原理并使用 Python、GDAL 和 NumPy 读取多波段遥感影像计算 NDVI 结果最后将结果保存为新的 GeoTIFF 文件。2 NDVI 基本含义NDVI 的全称是 Normalized Difference Vegetation Index中文通常称为归一化植被指数。它主要利用了植被在不同波段上的反射差异① 健康植被在近红外波段通常具有较强反射。② 健康植被在红光波段通常具有较强吸收。③ 当近红外反射较强、红光反射较弱时通常表示植被覆盖较好或长势较好。因此可以通过近红外波段和红光波段构造 NDVI用于初步判断植被分布和植被生长情况。3 NDVI 计算公式NDVI 的计算公式为NDVINIR−RedNIRRed NDVI\frac{NIR-Red}{NIRRed}NDVINIRRedNIR−Red​其中NIRNIRNIR表示近红外波段RedRedRed表示红光波段。从公式可以看出NDVI 本质上是一个归一化差值指数。理论上NDVI 的取值范围通常在[−1,1][-1,1][−1,1]之间。一般可以粗略理解为① NDVI 越接近1植被覆盖或植被生长状况通常越好。② NDVI 接近0可能表示裸地、建筑物或植被覆盖较少区域。③ NDVI 小于0通常可能对应水体、云、雪或其他非植被地物。需要注意NDVI 的具体解释会受到影像类型、传感器、季节、区域地物和预处理方式影响不能只依靠一个固定阈值机械判断。4 波段选择说明计算 NDVI 的关键是找到影像中的近红外波段和红光波段。不同遥感影像的波段顺序可能不同。例如某些影像中第 1 个波段可能是近红外第 2 个波段可能是红光但在其他数据中红光和近红外可能位于不同波段。因此计算前需要先确认影像的波段含义。例如假设当前影像读取为数组后①tif_array[0]表示近红外波段。②tif_array[1]表示红光波段。则可以写成nir_bandtif_array[0]red_bandtif_array[1]如果影像波段顺序不同需要根据实际数据调整索引。5 Python 实现思路整体计算流程如下① 使用 GDAL 打开遥感影像。② 读取影像宽度、高度、投影和仿射变换信息。③ 读取近红外波段和红光波段。④ 将波段数据转换为浮点类型。⑤ 根据公式计算 NDVI。⑥ 处理分母为0或无效值的像元。⑦ 创建新的 GeoTIFF 文件。⑧ 写入 NDVI 数组并保留原始影像的空间参考信息。6 完整代码下面代码使用 GDAL 和 NumPy 计算 NDVI并将结果保存为新的 GeoTIFF 文件。fromosgeoimportgdalimportnumpyasnp input_tif../jiaozuo_jiefangqu_crop.tifoutput_tif../jiaozuo_jiefangqu_crop_ndvi.tifdatasetgdal.Open(input_tif)ifdatasetisNone:raiseFileNotFoundError(f无法打开影像文件{input_tif})tif_widthdataset.RasterXSize tif_heightdataset.RasterYSize tif_arraydataset.ReadAsArray(0,0,tif_width,tif_height)# 需要根据实际影像波段顺序确认近红外和红光波段索引nir_bandtif_array[0].astype(np.float64)red_bandtif_array[1].astype(np.float64)denominatornir_bandred_band nodata_value9999ndvi_arraynp.full((tif_height,tif_width),nodata_value,dtypenp.float64)valid_maskdenominator!0ndvi_array[valid_mask]((nir_band[valid_mask]-red_band[valid_mask])/denominator[valid_mask])ndvi_array[(ndvi_array-1)|(ndvi_array1)]nodata_value drivergdal.GetDriverByName(GTiff)new_tiffdriver.Create(output_tif,xsizetif_width,ysizetif_height,bands1,eTypegdal.GDT_Float64)new_tiff.SetGeoTransform(dataset.GetGeoTransform())new_tiff.SetProjection(dataset.GetProjection())new_tiff.GetRasterBand(1).SetNoDataValue(nodata_value)new_tiff.GetRasterBand(1).WriteArray(ndvi_array)new_tiff.FlushCache()new_tiffNonedatasetNone7 代码说明7.1 打开影像datasetgdal.Open(input_tif)gdal.Open()用于打开栅格影像。如果路径错误或文件无法读取返回值可能为None因此代码中增加了判断ifdatasetisNone:raiseFileNotFoundError(f无法打开影像文件{input_tif})这样可以避免后续读取影像属性时报错不清晰。7.2 读取影像数组tif_arraydataset.ReadAsArray(0,0,tif_width,tif_height)对于多波段影像ReadAsArray()读取后的数组通常为三维数组结构可以理解为波段数 × 影像高度 × 影像宽度因此可以通过tif_array[0]、tif_array[1]这种方式读取不同波段。7.3 选择近红外和红光波段nir_bandtif_array[0].astype(np.float64)red_bandtif_array[1].astype(np.float64)这里假设第 1 个波段为近红外波段第 2 个波段为红光波段。需要注意实际使用时必须结合影像数据说明或波段信息确认波段顺序。如果波段顺序选错计算出来的 NDVI 就没有实际意义。7.4 处理分母为 0 的情况NDVI 公式的分母是NIRRedNIRRedNIRRed。如果某些像元中NIRRed0NIRRed0NIRRed0直接计算会导致除零问题。因此可以先构造有效像元掩膜valid_maskdenominator!0然后只对有效像元计算 NDVIndvi_array[valid_mask]((nir_band[valid_mask]-red_band[valid_mask])/denominator[valid_mask])无效像元保留为 NoData 值。7.5 保留空间参考信息输出 NDVI 影像时需要保留原始影像的空间位置和投影信息。new_tiff.SetGeoTransform(dataset.GetGeoTransform())new_tiff.SetProjection(dataset.GetProjection())其中①GetGeoTransform()用于获取仿射变换参数。②GetProjection()用于获取影像投影信息。这样输出的 NDVI 结果才能在 GIS 软件中正确叠加到原始影像位置。8 小结NDVI 是遥感影像处理中最常用的植被指数之一。它通过近红外波段和红光波段之间的差异来反映植被信息公式简单计算方便。在实际使用 Python 计算 NDVI 时需要重点注意三个问题确认近红外和红光波段顺序处理分母为0的无效像元并在输出结果时保留原始影像的空间参考信息。