批量计算遥感影像NDVI:Python代码
本文介绍基于Python中的gdal
模块,批量基于大量多波段遥感影像文件,计算其每1
景图像各自的NDVI数值,并将多景结果依次保存为栅格文件的方法。
如下图所示,现在有大量.tif
格式的遥感影像文件,其中均含有红光波段与近红外波段(此外也可以含有其他光谱波段,有没有都不影响);我们希望,批量计算其每1
景遥感影像的NDVI。
在之前的文章中,我们多次介绍过在不同软件或平台中计算NDVI的方法;而在本文中,我们就介绍一下基于Python中的gdal
模块,实现NDVI批量计算的方法。
这里所需的代码如下。
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 18 12:37:22 2024
@author: fkxxgis
"""
import os
from osgeo import gdal
original_folder = r"E:\04_Reconstruction\99_MODIS\new_data\2021_48STA_Result\Original"
output_folder = r"E:\04_Reconstruction\99_MODIS\new_data\2021_48STA_Result\NDVI_Original"
for filename in os.listdir(original_folder):
if filename.endswith('.tif'):
dataset = gdal.Open(os.path.join(original_folder, filename), gdal.GA_ReadOnly)
width = dataset.RasterXSize
height = dataset.RasterYSize
driver = gdal.GetDriverByName('GTiff')
output_dataset = driver.Create(os.path.join(output_folder, "NDVI_" + filename), width, height, 1, gdal.GDT_Float32)
band_red = dataset.GetRasterBand(3)
data_red = band_red.ReadAsArray()
data_red = data_red.astype(float)
band_nir = dataset.GetRasterBand(4)
data_nir = band_nir.ReadAsArray()
data_nir = data_nir.astype(float)
data_ndvi = (data_nir - data_red) / (data_nir + data_red)
output_band = output_dataset.GetRasterBand(1)
output_band.WriteArray(data_ndvi)
output_band.FlushCache()
output_dataset.SetGeoTransform(dataset.GetGeoTransform())
output_dataset.SetProjection(dataset.GetProjection())
dataset = None
output_dataset = None
print(filename, "finished!")
代码整体也非常简单。首先,我们定义输入文件与输入结果文件的路径,前者就是待计算NDVI的遥感影像文件路径,后者则是NDVI结果的遥感影像文件路径。
接下来,遍历original_folder
文件夹中的文件。其中,os.listdir()
用于获取文件夹中的文件列表,其后的endswith('.tif')
用于筛选出以.tif
扩展名结尾的文件。
随后,对于每个以.tif
结尾的文件,首先使用gdal.Open()
打开文件——其中的os.path.join()
用于构建完整的文件路径;接下来获取影像数据集的宽度和高度,并使用gdal.GetDriverByName()
获取GTiff驱动程序,用于创建输出影像文件;同时,使用driver.Create()
创建一个与原始影像具有相同大小的输出影像文件。
紧接着,从数据集中获取红光和近红外波段的数据。dataset.GetRasterBand()
用以获取指定的栅格波段,而band.ReadAsArray()
则将波段数据读取为数组;同时,我这里还用了astype()
转换数组的格式,避免原本遥感影像的数据格式带来的问题——例如,假如原本遥感影像是无符号整型的数据格式,那么这里不加astype()
计算NDVI就会有问题。
其次,即可计算NDVI。使用获取的红光和近红外波段数据计算NDVI,并将NDVI数据保存在data_ndvi
数组中。
最后,将NDVI数据写入输出影像文件。output_dataset.GetRasterBand()
获取输出影像文件的波段,band.WriteArray()
将数据写入波段,band.FlushCache()
刷新波段缓存。
此外,记得通过output_dataset.SetGeoTransform()
和output_dataset.SetProjection()
设置输出影像文件的地理变换和投影信息。
同时,需要清理和关闭数据集,将数据集和输出数据集设置为None
以释放资源。还可以打印文件名和finished!
,表示当前文件处理完成。
执行上述代码,我们即可在结果文件夹中看到计算得到的NDVI数据;如下图所示。
至此,大功告成。