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

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

至此，大功告成。

标签: none

添加新评论