Python ArcPy批次計算多時相遙感影像的各項元平均值

瘋狂學習GIS發表於2023-04-18

  本文介紹基於PythonArcPy模組,對大量長時間序列柵格遙感影像檔案的每一個像元進行多時序平均值的求取。

  在遙感應用中,我們經常需要對某一景遙感影像中的全部像元的畫素值進行平均值求取——這一操作很好實現,基於ArcMap軟體或者簡單的Python程式碼就可以實現;但有時候,我們會需要結合同一地區、不同時相多景遙感影像,求取每一個像元全部時相中畫素值的平均值——這一需求的實現較之前者就有些麻煩,本文對此加以介紹。

  首先,我們來明確一下本文的具體需求。現有一個儲存有大量.tif格式遙感影像的資料夾,其中每一個遙感影像的檔名中都包含有該影像的成像時間,如下圖所示。且其中除了.tif格式的遙感影像檔案外,還具有其它格式的檔案。

  我們希望,對於同一年成像的遙感影像進行逐像元平均值的求取。例如,上圖中具有2001年第185天成像、第193天成像、第201天成像……等等遙感影像8幅,每一幅都是這一年不同時間在同一空間位置的成像;同時,還有2005年不同時間成像的遙感影像9幅。我們希望,首先將2001年成像的8幅遙感影像加以逐像元平均值的求取,即求取每一個像元在這8景影像中畫素值的平均;隨後再對2005年成像的9幅遙感影像加以逐像元平均值的求取,以此類推。

  明確了需求後,我們就可以開始具體的操作。首先,本文所需用到的程式碼如下。

# -*- coding: utf-8 -*-
"""
Created on Sat Apr 16 10:48:37 2022

@author: fkxxgis
"""

import arcpy
from arcpy.sa import *

tif_file_path="E:/LST/Data/MODIS/05_Resample/"
average_file_path="E:/LST/Data/MODIS/06_Average/"
arcpy.env.workspace=tif_file_path

tif_file_name=arcpy.ListRasters("*","tif")
tif_file_year=tif_file_name[0][0:4]
one_year_tif_list=[]
sum_pic=0

for tif_file in tif_file_name:
    if tif_file[0:4]==tif_file_year:
        one_year_tif_list.append(tif_file)
        tif_file_temp=tif_file
        if tif_file==tif_file_name[len(tif_file_name)-1]:
            pic_num=len(one_year_tif_list)
            for tif_file_new in one_year_tif_list:
                sum_pic=sum_pic+Raster(tif_file_new)
            (sum_pic/pic_num).save(average_file_path+tif_file_year+"_Ave.tif")
    else:
        pic_num=len(one_year_tif_list)
        for tif_file_new in one_year_tif_list:
            sum_pic=sum_pic+Raster(tif_file_new)
        (sum_pic/pic_num).save(average_file_path+tif_file_year+"_Ave.tif")
        one_year_tif_list=[]
        sum_pic=0
        one_year_tif_list.append(tif_file)
        tif_file_year=tif_file[0:4]

  其中,tif_file_path是原有計算平均值前遙感影像的儲存路徑,average_file_path是我們新生成的求取平均值後遙感影像的儲存路徑,也就是結果儲存路徑。

  在這裡,和我們前期的部落格Python ArcPy批次拼接長時間序列柵格影像類似,需要首先在資源管理器中,將tif_file_path路徑下的各檔案以“名稱”排序的方式進行排序;隨後,利用arcpy.ListRasters()函式,獲取路徑下原有的全部.tif格式的影像檔案,並擷取第一個檔案的部分檔名,從而獲取其成像時間的具體年份。

  接下來,遍歷tif_file_path路徑下全部.tif格式影像檔案。其中,我們透過一個簡單的判斷語句if tif_file[0:4]==tif_file_year:,來確定某一年的遙感影像是否已經讀取完畢——如果已經讀取完畢,例如假如2001年成像的8幅遙感影像都已經遍歷過了,那麼就對這8景遙感影像加以逐像元的平均值求取,並開始對下一個年份(即2005年)成像的遙感影像繼續加以計算;如果還沒有讀取完畢,例如假如2001年成像的8幅遙感影像目前僅遍歷到了第5幅,那麼就不求平均值,繼續往下遍歷,直到遍歷完2001年成像的8幅遙感影像。

  這裡相信大家也看到了為什麼我們要在前期先將資料夾中的檔案按照“名稱”排序——是為了保證同一年成像的所有遙感影像都排列在一起,遍歷時只要遇到一個新的年份,程式就知道上一個年份的所有影像都已經遍歷完畢了,就可以將上一個年份的所有柵格影像加以平均值求取。

  在這裡,逐像元的平均值求取其實也非常簡單——我們對每一個像元分別執行以下操作:首先將該像元在當前年份裡所有遙感影像的畫素值相加,隨後除以這一年份的遙感影像的數量,得到的就是該像元在這一年中畫素值的平均值

  最後,透過if tif_file==tif_file_name[len(tif_file_name)-1]:這個判斷,來確認是否目前已經遍歷到資料夾中的最後一個影像檔案。如果是的話,就需要將當前成像年份的所有影像進行平均值的求取,並宣告程式碼完成執行。

  在 IDLE (Python GUI) 中執行程式碼。程式碼執行完畢後,我們看一下結果資料夾。可以看到,其中的影像已經是按照成像時間,分別完成平均值求取後的結果了。

  在最後,還需要說明一點——用以上程式碼來求取長時間序列遙感影像的像元平均值,對於任意一個像元,只要該像元在任意一個時相的影像中是無效值(即為NoData),那麼該像元在最終求出的平均值結果圖中,畫素值也將會是無效值NoData。針對這一問題的解決,我們將在下一篇部落格中介紹。

相關文章