批次計算遙感影像NDVI：Python程式碼

疯狂学习GIS發表於2024-11-09
Python

  本文介紹基於Python中的gdal模組，批次基於大量多波段遙感影像檔案，計算其每1景影像各自的NDVI數值，並將多景結果依次儲存為柵格檔案的方法。

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

image

  在之前的文章中，我們多次介紹過在不同軟體或平臺中計算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資料；如下圖所示。

  至此，大功告成。

相關文章