??本文介紹基于Python中的gdal模塊，批量基于大量多波段遙感影像文件，計算其每1景圖像各自的NDVI數(shù)值，并將多景結(jié)果依次保存為柵格文件的方法。

??如下圖所示，現(xiàn)在有大量.tif格式的遙感影像文件，其中均含有紅光波段與近紅外波段（此外也可以含有其他光譜波段，有沒有都不影響）；我們希望，批量計算其每1景遙感影像的NDVI。

??在之前的文章中，我們多次介紹過在不同軟件或平臺中計算NDVI的方法，大家可以參考文章ArcGIS自動計算單波段、多波段柵格圖像的NDVI，或者文章Google Earth Engine谷歌地球引擎柵格代數(shù)與波段計算求取NDVI。而在本文中，我們就介紹一下基于Python中的gdal模塊，實現(xiàn)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!")

??代碼整體也非常簡單。首先，我們定義輸入文件與輸入結(jié)果文件的路徑，前者就是待計算NDVI的遙感影像文件路徑，后者則是NDVI結(jié)果的遙感影像文件路徑。

??接下來，遍歷original_folder文件夾中的文件。其中，os.listdir()用于獲取文件夾中的文件列表，其后的endswith('.tif')用于篩選出以.tif擴(kuò)展名結(jié)尾的文件。

??隨后，對于每個以.tif結(jié)尾的文件，首先使用gdal.Open()打開文件——其中的os.path.join()用于構(gòu)建完整的文件路徑；接下來獲取影像數(shù)據(jù)集的寬度和高度，并使用gdal.GetDriverByName()獲取GTiff驅(qū)動程序，用于創(chuàng)建輸出影像文件；同時，使用driver.Create()創(chuàng)建一個與原始影像具有相同大小的輸出影像文件。

??緊接著，從數(shù)據(jù)集中獲取紅光和近紅外波段的數(shù)據(jù)。dataset.GetRasterBand()用以獲取指定的柵格波段，而band.ReadAsArray()則將波段數(shù)據(jù)讀取為數(shù)組；同時，我這里還用了astype()轉(zhuǎn)換數(shù)組的格式，避免原本遙感影像的數(shù)據(jù)格式帶來的問題——例如，假如原本遙感影像是無符號整型的數(shù)據(jù)格式，那么這里不加astype()計算NDVI就會有問題。

??其次，即可計算NDVI。使用獲取的紅光和近紅外波段數(shù)據(jù)計算NDVI，并將NDVI數(shù)據(jù)保存在data_ndvi數(shù)組中。

??最后，將NDVI數(shù)據(jù)寫入輸出影像文件。output_dataset.GetRasterBand()獲取輸出影像文件的波段，band.WriteArray()將數(shù)據(jù)寫入波段，band.FlushCache()刷新波段緩存。

??此外，記得通過output_dataset.SetGeoTransform()和output_dataset.SetProjection()設(shè)置輸出影像文件的地理變換和投影信息。

??同時，需要清理和關(guān)閉數(shù)據(jù)集，將數(shù)據(jù)集和輸出數(shù)據(jù)集設(shè)置為None以釋放資源。還可以打印文件名和finished!，表示當(dāng)前文件處理完成。

??執(zhí)行上述代碼，我們即可在結(jié)果文件夾中看到計算得到的NDVI數(shù)據(jù)；如下圖所示。

??至此，大功告成。