Python GDAL绘制栅格像元数值变化曲线图

  本文介绍基于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()函数将生成的图片保存在指定路径下。

  最终，我们得到的多张曲线图结果如下图所示，其文件名通过列号与行号分别表示了当前这张图是基于哪一个像元绘制得到的；其中，每一张图的具体样式就是本文开头所展示的那一张图片的样子。

  至此，大功告成。