归一化植被指数:遥感影像中,近红外波段的反射值与红光波段的反射值之差比上两者之和。
计算公式:NDVI = (NIR-R)/(NIR+R)
注:NIR为近红外波段的反射值
R为红光波段的反射值。
由于学习需要,从网上寻找方法试图批量计算该指数,防止以后忘记。附代码如下:
(主要使用gdal库)
import osfrom PIL import Imageimport numpy as npfrom osgeo import gdalimport globimport cv2list_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.0001nir = in_ds.GetRasterBand(2).ReadAsArray() * 0.0001ndvi = (nir - red) / (nir + red)# 将NAN转化为0值nan_index = np.isnan(ndvi)ndvi[nan_index] = 0ndvi = 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./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./tanlangqie/article/details/79486144