??本文介紹基于Python中gdal模塊，對(duì)大量多時(shí)相柵格圖像，批量繪制像元時(shí)間序列折線圖的方法。

??首先，明確一下本文需要實(shí)現(xiàn)的需求：現(xiàn)有三個(gè)文件夾，其中第一個(gè)文件夾存放了某一研究區(qū)域原始的多時(shí)相柵格遙感影像數(shù)據(jù)（每一景遙感影像對(duì)應(yīng)一個(gè)時(shí)相，文件夾中有多景遙感影像），每一景遙感影像都是.tif格式；第二個(gè)文件夾與第三個(gè)文件夾則分別存放了前述第一個(gè)文件夾中原始遙感影像基于2種不同濾波方法處理后的遙感影像（同樣是每一景遙感影像對(duì)應(yīng)一個(gè)時(shí)相，文件夾中有多景遙感影像），每一景遙感影像同樣也都是.tif格式。我們希望分別針對(duì)這三個(gè)文件夾中的多張遙感影像數(shù)據(jù)，隨機(jī)繪制部分像元對(duì)應(yīng)的時(shí)間序列曲線圖（每一個(gè)像元對(duì)應(yīng)一張曲線圖，一張曲線圖中有三條曲線）；每一張曲線圖的最終結(jié)果都是如下所示的類似的樣式，X軸表示時(shí)間節(jié)點(diǎn)，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為原始多時(shí)相遙感影像數(shù)據(jù)存放路徑，也就是前述的第一個(gè)文件夾的路徑；而E:\AllYear\R與E:\AllYear\S則是前述第二個(gè)文件夾和第三個(gè)文件夾對(duì)應(yīng)的路徑；E:\AllYear\Pic則是批量繪圖后，圖片保存的路徑。這里請(qǐng)注意，在運(yùn)行代碼前我們需要在資源管理器中，將上述三個(gè)路徑下的各文件以“名稱”排序的方式進(jìn)行排序（每一景遙感影像都是按照成像時(shí)間命名的）。此外，pic_num則是需要加以繪圖的像元個(gè)數(shù)，也就表明后期我們所生成的曲線圖的張數(shù)為50。

??代碼的整體思路也非常簡單。首先，我們借助os.listdir()函數(shù)獲取original_file_path路徑下的所有柵格遙感影像文件，在基于gdal.Open()函數(shù)將這一文件下的第一景遙感影像打開后，獲取其行數(shù)與列數(shù)；隨后，通過np.random.randint()函數(shù)生成兩個(gè)隨機(jī)數(shù)數(shù)組，分別對(duì)應(yīng)著后期我們繪圖的像元的行號(hào)與列號(hào)。

??在代碼的下一部分（就是hants_file_list開頭的這一部分），我們是通過截取文件夾中圖像的名稱，來確定后期我們生成的時(shí)間序列曲線圖中X軸的標(biāo)簽（也就是每一個(gè)x對(duì)應(yīng)的時(shí)間節(jié)點(diǎn)是什么）——其中，這里的[12:15]就表示對(duì)于我的柵格圖像而言，其文件名的第13到15個(gè)字符表示了遙感影像的成像時(shí)間；大家在使用代碼時(shí)依據(jù)自己的實(shí)際情況加以修改即可。在這里，我們得到的day_list，就是后期曲線圖中X軸各個(gè)標(biāo)簽的內(nèi)容。

??隨后，代碼中最外層的for循環(huán)部分，即為批量繪圖工作的開始。我們前面選擇好了50個(gè)隨機(jī)位置的像元，此時(shí)就可以遍歷這些像元，對(duì)每一個(gè)像元在不同時(shí)相中的數(shù)值加以讀取——通過.ReadAsArray()函數(shù)將柵格圖像各波段的信息讀取為Array格式，并通過對(duì)應(yīng)的行號(hào)與列號(hào)加以像素值的獲??；隨后，將獲取得到的像元在不同時(shí)相的數(shù)值通過.append()函數(shù)依次放入前面新生成的列表中。

??在接下來，即可開始繪圖的工作。其中，pic_file_name表示每一張曲線圖的文件名稱，這是通過當(dāng)前像元對(duì)應(yīng)的行號(hào)與列號(hào)來命名的；plt.figure(dpi = 300)表示設(shè)置繪圖的DPI為300。隨后，再對(duì)每一張曲線圖的圖名、圖例與坐標(biāo)軸標(biāo)簽等加以配置，并通過plt.savefig()函數(shù)將生成的圖片保存在指定路徑下。

??最終，我們得到的多張曲線圖結(jié)果如下圖所示，其文件名通過列號(hào)與行號(hào)分別表示了當(dāng)前這張圖是基于哪一個(gè)像元繪制得到的；其中，每一張圖的具體樣式就是本文開頭所展示的那一張圖片的樣子。

??至此，大功告成。