Python批量计算遥感图像NDVI（归一化植被指数）

归一化植被指数：遥感影像中，近红外波段的反射值与红光波段的反射值之差比上两者之和。

计算公式：NDVI = (NIR-R)/(NIR+R)

注：NIR为近红外波段的反射值 

       R为红光波段的反射值。

由于学习需要，从网上寻找方法试图批量计算该指数，防止以后忘记。附代码如下：

（主要使用gdal库）

import os
from PIL import Image
import numpy as np
from osgeo import gdal
import glob
import cv2

list_tif = glob.glob('H:/gdal/test-data/cut-test/*.tif')
out_path = 'H:/gdal/test-data/ndvi-test/'
 
for tif in list_tif:
    in_ds = gdal.Open(tif)
    # 获取文件所在路径以及不带后缀的文件名
    (filepath, fullname) = os.path.split(tif)
    (prename, suffix) = os.path.splitext(fullname)
    if in_ds is None:
        print('Could not open the file ' + tif)
    else:
        # 将MODIS原始数据类型转化为反射率
        red = in_ds.GetRasterBand(1).ReadAsArray() * 0.0001
        nir = in_ds.GetRasterBand(2).ReadAsArray() * 0.0001
        ndvi = (nir - red) / (nir + red)
        # 将NAN转化为0值
        nan_index = np.isnan(ndvi)
        ndvi[nan_index] = 0
        ndvi = ndvi.astype(np.float32)
        # 将计算好的NDVI保存为GeoTiff文件
        gtiff_driver = gdal.GetDriverByName('GTiff')
        # 批量处理需要注意文件名是变量，这里截取对应原始文件的不带后缀的文件名
        out_ds = gtiff_driver.Create(out_path + prename + '.tif',
                         ndvi.shape[1], ndvi.shape[0], 1, gdal.GDT_Float32)
        # 将NDVI数据坐标投影设置为原始坐标投影
        out_ds.SetProjection(in_ds.GetProjection())
        out_ds.SetGeoTransform(in_ds.GetGeoTransform())
        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(ndvi)
        out_band.FlushCache()

list_tif：需要计算的tif图像所在文件夹

out_path：输出NDVI图像所在文件夹

 

 

参考文章：

1.https://baike.baidu.com/item/%E5%BD%92%E4%B8%80%E5%8C%96%E6%A4%8D%E8%A2%AB%E6%8C%87%E6%95%B0/625296?fr=aladdin

2.https://blog.csdn.net/tanlangqie/article/details/79486144