?首先展示单次针对波段分开的TIFF影像的NDVI植被指数的计算
代码如下：（代码中已经将NDVI指数异常值进行了剔除，取值范围最终在[-1, 1]区间内）
from osgeo import gdal, gdalconst
import os
import numpy as np
def NDVI_calculte(Red_band_file, NIR_band_file, NDVI_output_file):
"""
红光与近红外波段分别存储为tif文件的NDVI指数计算
:param Red_band_file: 红光波段文件
:param NIR_band_file: 近红外波段文件
:param NDVI_output_file: 输出NDVI影像
:return:
"""
band_Red = gdal.Open(Red_band_file, gdalconst.GA_ReadOnly)
band_NIR = gdal.Open(NIR_band_file, gdalconst.GA_ReadOnly)
cols_Red = band_Red.RasterXSize # 列数
rows_Red = band_Red.RasterYSize # 行数
cols_NIR = band_NIR.RasterXSize # 列数
rows_NIR = band_NIR.RasterYSize # 行数
data_type_Red = band_Red.GetRasterBand(1).DataType
data_type_NIR = band_NIR.GetRasterBand(1).DataType
geotrans = list(band_Red.GetGeoTransform())
if (cols_Red != cols_NIR or rows_Red != rows_NIR or data_type_Red != data_type_NIR):
print("输入红光波段与近红外波段影像参数不对应")
return
if os.path.exists(NDVI_output_file) and os.path.isfile(NDVI_output_file): # 如果已存在同名影像
os.remove(NDVI_output_file) # 则删除之
NDVI = band_Red.GetDriver().Create(NDVI_output_file, xsize=cols_Red, ysize=rows_Red, bands=1,
eType=data_type_Red) # 创建空的与待处理的tif文件
NDVI.SetProjection(band_Red.GetProjection()) # 设置投影坐标
NDVI.SetGeoTransform(geotrans) # 设置地理变换参数
data_Red = band_Red.ReadAsArray(buf_xsize=cols_Red, buf_ysize=rows_Red)
data_NIR = band_NIR.ReadAsArray(buf_xsize=cols_NIR, buf_ysize=rows_NIR)
data_NDVI = data_NIR
count = 0
list_x = []
list_y = []
for i in range(0, rows_Red):
for j in range(0, cols_Red):
if (data_NIR[i][j]+data_Red[i][j] != 0 and data_Red[i][j]!=None and data_NIR[i][j]!=None):
data_NDVI[i][j] = (data_NIR[i][j]-data_Red[i][j])/(data_NIR[i][j]+data_Red[i][j])
if (data_NDVI[i][j]<-1 or data_NDVI[i][j]>1): #剔除异常值，使最终NDVI范围在[-1, 1]内
data_NDVI[i][j] = None
count += 1 #记录异常值的个数，
list_x.append(i+1)
list_y.append(j+1) #记录异常值像素的X,Y坐标
else:
data_NDVI[i][j] = None
NDVI.GetRasterBand(1).WriteArray(data_NDVI) # 写入数据到新影像中
NDVI.GetRasterBand(1).FlushCache()
NDVI.GetRasterBand(1).ComputeBandStats(False) # 计算统计信息
print("已剔除异常值，异常值个数为"+str(count)+"个")
print("正在写入完成")
del band_Red
del band_NIR
if __name__ == "__main__":
band_Red = "D:/MOD09A1_06_09/Result_Link/mean_Bo_1000/Bo1_mean_2019_06_1000.tif" #红光波段文件
band_NIR = "D:/MOD09A1_06_09/Result_Link/mean_Bo_1000/Bo2_mean_2019_06_1000.tif" #近红外波段文件
NDVI_output_file = "D:/MOD09A1_06_09/Result_Link/NDVI_6789/NDVI_2019_06.tif " #输出NDVI指数文件
NDVI_calculte(band_Red, band_NIR, NDVI_output_file)展示计算后ArcGIS内的影像截图：
注：该图为Sin投影数据框属性切换为WGS_84后的影像
查看其NDVI直方图如下：
可以看出计算后效果很好。
针对波段合成后的TIFF影像计算NDVI处理影像
只需要修改主函数内代码即可
对于MODIS——MOD09A1影像（band1为红光波段，band2为近红外波段）
?
if __name__ == "__main__":
band_All = "D:/MOD09A1_06_09/Result_Link/mean_Bo_1000/Bo1_mean_2019_06_1000.tif" #输入合成多波段后文件
NDVI_output_file = "D:/MOD09A1_06_09/Result_Link/NDVI_6789/NDVI_2019_06.tif " #输出NDVI指数文件
NDVI_calculte(band_All.GetRasterBand(1), band_All.GetRasterBand(2), NDVI_output_file)对于Landsat 4/5 TM波段影像（band3为红光波段，band4为近红外波段）?
if __name__ == "__main__":
band_All = "D:/MOD09A1_06_09/Result_Link/mean_Bo_1000/Bo1_mean_2019_06_1000.tif" #输入合成多波段后文件
NDVI_output_file = "D:/MOD09A1_06_09/Result_Link/NDVI_6789/NDVI_2019_06.tif " #输出NDVI指数文件
NDVI_calculte(band_All.GetRasterBand(3), band_All.GetRasterBand(4), NDVI_output_file)针对波段合成后的批量TIFF影像计算NDVI处理
?（以MODIS——MOD09A1影像为例）
if __name__ == "__main__":
file_dir = "D:/MOD09A1_06_09/Result_Link/mean_Bo_1000" #输入合成多波段后文件路径
list_dir = os.listdir(file_dir) #以列表形式返回存储tif文件的文件夹内所有文件
number_tif = len(list_dir) #待处理影像数目
for i in range(0, number_tif):
band_All = list_dir[i]
NDVI_output_file = "D:/MOD09A1_06_09/Result_Link/NDVI_6789/NDVI_2019_0"+str(i+1)+".tif " #输出NDVI指数文件
NDVI_calculte(band_All.GetRasterBand(1), band_All.GetRasterBand(2), NDVI_output_file)以上就是关于tif影像的分开波段和多波段植被指数（NDVI）计算方法！