??本文介紹基于Python中的gdal模塊，對大量長時間序列的柵格遙感影像文件，繪制其每一個波段中、若干隨機指定的像元的時間序列曲線圖的方法。

??在之前的文章Python GDAL繪制柵格像元數(shù)值變化曲線圖（http://www.itdecent.cn/p/344e2ccba0e1）中，我們就已經(jīng)介紹過基于gdal模塊，對大量多時相柵格圖像，批量繪制像元時間序列折線圖的方法。不過當時文章中的需求，每1個時相都對應著3個不同的遙感影像文件，而每1個遙感影像文件則都僅僅只有1個波段；而在本文中，我們每1景遙感影像都對應著2個波段，我們最終繪制的多條曲線圖，也都來自于這每1景遙感影像的不同波段。

??我們再來明確一下本文的需求?，F(xiàn)在有一個文件夾，其中放置了大量的遙感影像文件，如下圖所示。其中，所有遙感影像都是同一地區(qū)、不同成像時間的圖像，其各自的空間參考信息、像元行數(shù)與列數(shù)等都是一致的，文件名中有表示成像日期的具體字段；且每1景遙感影像都具有2個波段。現(xiàn)在我們希望，在遙感影像覆蓋的區(qū)域內，隨機選取若干的像元，基于這些像元，我們繪制其隨時間變化的曲線圖。因為我們的每個遙感影像都有2個波段，且都希望繪制出曲線圖，因此最終的曲線圖一共就有2條曲線。

??明確了需求，我們就可以開始代碼的撰寫。本文用到的代碼如下。

# -*- coding: utf-8 -*-
"""
Created on Tue Jul 25 23:04:41 2023

@author: fkxxgis
"""

import os
import random
import matplotlib.pyplot as plt
from osgeo import gdal

def load_image(image_path):
    dataset = gdal.Open(image_path)
    band1 = dataset.GetRasterBand(1).ReadAsArray()
    band2 = dataset.GetRasterBand(2).ReadAsArray()
    del dataset
    return band1, band2

def plot_time_series(image_folder, pic_folder, num_pixels):
    image_files = [file for file in os.listdir(image_folder) if file.endswith(".tif")]
    band1_merge, band2_merge = [], []
    i = 0
    
    for image_file in image_files:
        image_path = os.path.join(image_folder, image_file)
        band1, band2 = load_image(image_path)
        band1_merge.append(band1)
        band2_merge.append(band2)
        i += 1

    x_size, y_size = band1.shape
    pixel_indices = random.sample(range(x_size * y_size), num_pixels)

    for pixel_index in pixel_indices:
        x, y = divmod(pixel_index, y_size)
        band_list_1, band_list_2 = [], []
        for i in range(len(band1_merge)):
            band_data_1 = band1_merge[i]
            band_list_1.append(band_data_1[x, y])
            band_data_2 = band2_merge[i]
            band_list_2.append(band_data_2[x, y])

        plt.figure()
        plt.plot(range(len(band1_merge)), band_list_1, label="Band 1")
        plt.plot(range(len(band1_merge)), band_list_2, label="Band 2")
        plt.xlabel("Time")
        plt.ylabel("NDVI")
        plt.ylim(0, 1000)
        plt.title(f"Time Series for Pixel at ({x}, {y})")
        plt.legend()
        plt.savefig(os.path.join(pic_folder, str(x) + "_" + str(y)))
        plt.show()

image_folder_path = "E:/02_Project/Result/test"
pic_folder_path = "E:/02_Project/TIFF/Plot"
num_pixels = 50
plot_time_series(image_folder_path, pic_folder_path, num_pixels)

??上述代碼的具體含義如下。

??首先，我們導入了需要使用的庫；其中，os用于處理文件路徑和目錄操作，random用于隨機選擇像素，matplotlib.pyplot則用于繪制圖像。

??隨后，我們定義函數(shù)load_image(image_path)；這個函數(shù)接收一個影像文件路徑image_path作為輸入?yún)?shù)。隨后，在函數(shù)內使用gdal庫打開該影像文件，然后提取其第一個和第二個波段的數(shù)據(jù)，并分別存儲在band1和band2中。最后，函數(shù)返回這兩個波段的數(shù)據(jù)。

??接下來，我們定義函數(shù)plot_time_series(image_folder, pic_folder, num_pixels)；這個函數(shù)接收三個輸入?yún)?shù)，分別為image_folder、pic_folder和num_pixels。其中，image_folder為包含多個.tif格式的影像文件的文件夾路徑，pic_folder是保存生成的時間序列圖像的文件夾路徑，而num_pixels則指定了隨機選擇的像素數(shù)量，用于繪制時間序列圖——這個參數(shù)設置為幾，我們最后就會得到幾張結果圖像。

??在這個函數(shù)的內部，我們通過os.listdir函數(shù)獲取image_folder中所有以.tif結尾的影像文件，并將這些文件名存儲在image_files列表中。然后，我們創(chuàng)建兩個空列表band1_merge和band2_merge，用于存儲所有影像文件的2個波段數(shù)據(jù)。接下來，我們遍歷所有影像文件，逐個加載每個影像文件的全部波段數(shù)據(jù)，并將它們添加到對應的列表中。其次，使用random.sample函數(shù)從像素索引的范圍中隨機選擇num_pixels個像素的索引，并保存在pixel_indices列表中。接下來，我們遍歷并恢復pixel_indices中的每個像素索引，計算該像素在每個影像中的每個波段的時間序列數(shù)據(jù)，并存儲在band_list_1、band_list_2列表中。

??隨后，我們即可繪制兩個時間序列圖，分別表示2個波段在不同影像日期上的數(shù)值。最后，我們將圖像保存到指定的文件夾pic_folder中，命名規(guī)則為x_y，其中x與y分別代表像素的橫、縱坐標。

??執(zhí)行上述代碼，我們即可在指定的文件夾路徑下看到我們生成的多張曲線圖；如下圖所示。

??其中，每1張圖像都表示了我們2個波段在這段時間內的數(shù)值走勢；如下圖所示。

??至此，大功告成。