基於Python GDAL為長時間序列遙感影像繪製時相變化曲線圖

疯狂学习GIS發表於2024-02-28
Python

  本文介紹基於Python中gdal模組，對大量多時相柵格影像，批次繪製像元時間序列折線圖的方法。

  首先，明確一下本文需要實現的需求：現有三個資料夾，其中第一個資料夾存放了某一研究區域原始的多時相柵格遙感影像資料（每一景遙感影像對應一個時相，資料夾中有多景遙感影像），每一景遙感影像都是.tif格式；第二個資料夾與第三個資料夾則分別存放了前述第一個資料夾中原始遙感影像基於2種不同濾波方法處理後的遙感影像（同樣是每一景遙感影像對應一個時相，資料夾中有多景遙感影像），每一景遙感影像同樣也都是.tif格式。我們希望分別針對這三個資料夾中的多張遙感影像資料，隨機繪製部分像元對應的時間序列曲線圖（每一個像元對應一張曲線圖，一張曲線圖中有三條曲線）；每一張曲線圖的最終結果都是如下所示的類似的樣式，X軸表示時間節點，Y軸就是具體的畫素值。

  知道了需求，我們便開始程式碼的書寫。具體程式碼如下：

# -*- coding: utf-8 -*-
"""
Created on Wed Dec 14 00:48:48 2022

@author: fkxxgis
"""

import os
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal

original_file_path = r"E:\AllYear\Original"
hants_file_path = r"E:\AllYear\Reconstruction"
sg_file_path = r"E:\AllYear\SG"
pic_file_path = r"E:\AllYear\Pic"
pic_num = 50
np.random.seed(6)

original_file_list = os.listdir(original_file_path)
tem_raster = gdal.Open(os.path.join(original_file_path, original_file_list[0]))
col_num = tem_raster.RasterXSize
row_num = tem_raster.RasterYSize
col_point_array = np.random.randint(0, col_num, pic_num)
row_point_array = np.random.randint(0, row_num, pic_num)
del tem_raster

hants_file_list = os.listdir(hants_file_path)
start_day = hants_file_list[0][12:15]
end_day = hants_file_list[-1][12:15]
day_list = [x for x in range(int(start_day), int(end_day) + 20, 10)]

for i in range(pic_num):
    original_pixel_list, hants_pixel_list, sg_pixel_list = [[] for x in range(3)]
    
    for tif in original_file_list:
        original_raster = gdal.Open(os.path.join(original_file_path, tif))
        original_array = original_raster.ReadAsArray()
        original_pixel_list.append(original_array[row_point_array[i],col_point_array[i]])
        
    for tif in hants_file_list:
        hants_raster = gdal.Open(os.path.join(hants_file_path, tif))
        hants_array = hants_raster.ReadAsArray()
        hants_pixel_list.append(hants_array[1, row_point_array[i],col_point_array[i]])
        
    sg_file_list = os.listdir(sg_file_path)
    for tif in sg_file_list:
        sg_raster = gdal.Open(os.path.join(sg_file_path, tif))
        sg_array = sg_raster.ReadAsArray()
        sg_pixel_list.append(sg_array[1, row_point_array[i],col_point_array[i]])
    
    pic_file_name = str(col_point_array[i]) + "_" + str(row_point_array[i]) + ".png"
    plt.figure(dpi = 300)
    plt.plot(original_pixel_list,color = "red", label = "Original")
    plt.plot(hants_pixel_list,color = "green", label = "HANTS")
    plt.plot(sg_pixel_list,color = "blue", label = "SG")
    plt.legend()
    plt.xticks(range(len(day_list)), day_list, fontsize = 11)
    plt.xticks(rotation = 45)
    plt.title(str(col_point_array[i]) + "_" + str(row_point_array[i]), fontweight = "bold")
    plt.savefig(os.path.join(pic_file_path, pic_file_name))
    plt.show()
    plt.clf()
    
    del original_raster
    del hants_raster
    del sg_raster

  其中，E:\AllYear\Original為原始多時相遙感影像資料存放路徑，也就是前述的第一個資料夾的路徑；而E:\AllYear\R與E:\AllYear\S則是前述第二個資料夾和第三個資料夾對應的路徑；E:\AllYear\Pic則是批次繪圖後，圖片儲存的路徑。這裡請注意，在執行程式碼前我們需要在資源管理器中，將上述三個路徑下的各檔案以“名稱”排序的方式進行排序（每一景遙感影像都是按照成像時間命名的）。此外，pic_num則是需要加以繪圖的像元個數，也就表明後期我們所生成的曲線圖的張數為50。

  程式碼的整體思路也非常簡單。首先，我們藉助os.listdir()函式獲取original_file_path路徑下的所有柵格遙感影像檔案，在基於gdal.Open()函式將這一檔案下的第一景遙感影像開啟後，獲取其行數與列數；隨後，透過np.random.randint()函式生成兩個隨機數陣列，分別對應著後期我們繪圖的像元的行號與列號。

  在程式碼的下一部分（就是hants_file_list開頭的這一部分），我們是透過擷取資料夾中影像的名稱，來確定後期我們生成的時間序列曲線圖中X軸的標籤（也就是每一個x對應的時間節點是什麼）——其中，這裡的[12:15]就表示對於我的柵格影像而言，其檔名的第13到15個字元表示了遙感影像的成像時間；大家在使用程式碼時依據自己的實際情況加以修改即可。在這裡，我們得到的day_list，就是後期曲線圖中X軸各個標籤的內容。

  隨後，程式碼中最外層的for迴圈部分，即為批次繪圖工作的開始。我們前面選擇好了50個隨機位置的像元，此時就可以遍歷這些像元，對每一個像元在不同時相中的數值加以讀取——透過.ReadAsArray()函式將柵格影像各波段的資訊讀取為Array格式，並透過對應的行號與列號加以畫素值的獲取；隨後，將獲取得到的像元在不同時相的數值透過.append()函式依次放入前面新生成的列表中。

  在接下來，即可開始繪圖的工作。其中，pic_file_name表示每一張曲線圖的檔名稱，這是透過當前像元對應的行號與列號來命名的；plt.figure(dpi = 300)表示設定繪圖的DPI為300。隨後，再對每一張曲線圖的圖名、圖例與座標軸標籤等加以配置，並透過plt.savefig()函式將生成的圖片儲存在指定路徑下。

  最終，我們得到的多張曲線圖結果如下圖所示，其檔名透過列號與行號分別表示了當前這張圖是基於哪一個像元繪製得到的；其中，每一張圖的具體樣式就是本文開頭所展示的那一張圖片的樣子。

  至此，大功告成。

相關文章