Python中GDAL批量绘制多时相栅格遥感影像的像元时间序列曲线图

简介: Python中GDAL批量绘制多时相栅格遥感影像的像元时间序列曲线图

  本文介绍基于Pythongdal模块,对大量多时相栅格图像,批量绘制像元时间序列折线图的方法。

  首先,明确一下本文需要实现的需求:现有三个文件夹,其中第一个文件夹存放了某一研究区域原始的多时相栅格遥感影像数据(每一景遥感影像对应一个时相,文件夹中有多景遥感影像),每一景遥感影像都是.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\RE:\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]就表示对于我的栅格图像而言,其文件名的第1315个字符表示了遥感影像的成像时间;大家在使用代码时依据自己的实际情况加以修改即可。在这里,我们得到的day_list,就是后期曲线图中X轴各个标签的内容。

  随后,代码中最外层的for循环部分,即为批量绘图工作的开始。我们前面选择好了50个随机位置的像元,此时就可以遍历这些像元,对每一个像元在不同时相中的数值加以读取——通过.ReadAsArray()函数将栅格图像各波段的信息读取为Array格式,并通过对应的行号列号加以像素值的获取;随后,将获取得到的像元在不同时相的数值通过.append()函数依次放入前面新生成的列表中。

  在接下来,即可开始绘图的工作。其中,pic_file_name表示每一张曲线图的文件名称,这是通过当前像元对应的行号列号来命名的;plt.figure(dpi = 300)表示设置绘图的DPI300。随后,再对每一张曲线图的图名、图例与坐标轴标签等加以配置,并通过plt.savefig()函数将生成的图片保存在指定路径下。

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

  至此,大功告成。

欢迎关注:疯狂学习GIS

相关文章
|
2月前
|
机器学习/深度学习 数据采集 算法
时间序列结构变化分析:Python实现时间序列变化点检测
在时间序列分析和预测中,准确检测结构变化至关重要。新出现的分布模式往往会导致历史数据失去代表性,进而影响基于这些数据训练的模型的有效性。
149 1
|
3月前
|
机器学习/深度学习 算法 数据挖掘
6种有效的时间序列数据特征工程技术(使用Python)
在本文中,我们将探讨使用日期时间列提取有用信息的各种特征工程技术。
131 0
|
22天前
|
机器学习/深度学习 数据采集 数据挖掘
11种经典时间序列预测方法:理论、Python实现与应用
本文将总结11种经典的时间序列预测方法,并提供它们在Python中的实现示例。
58 2
11种经典时间序列预测方法:理论、Python实现与应用
|
14天前
|
数据可视化 算法 JavaScript
基于图论的时间序列数据平稳性与连通性分析:利用图形、数学和 Python 揭示时间序列数据中的隐藏模式
本文探讨了如何利用图论分析时间序列数据的平稳性和连通性。通过将时间序列数据转换为图结构,计算片段间的相似性,并构建连通图,可以揭示数据中的隐藏模式。文章介绍了平稳性的概念,提出了基于图的平稳性度量,并展示了图分区在可视化平稳性中的应用。此外,还模拟了不同平稳性和非平稳性程度的信号,分析了图度量的变化,为时间序列数据分析提供了新视角。
32 0
基于图论的时间序列数据平稳性与连通性分析:利用图形、数学和 Python 揭示时间序列数据中的隐藏模式
|
2月前
|
机器学习/深度学习 索引 Python
python之序列
python之序列
143 59
|
23天前
|
存储 编译器 索引
Python 序列类型(2)
【10月更文挑战第8天】
Python 序列类型(2)
|
24天前
|
存储 C++ 索引
Python 序列类型(1)
【10月更文挑战第8天】
|
2月前
|
机器学习/深度学习 数据采集 算法
数据稀缺条件下的时间序列微分:符号回归(Symbolic Regression)方法介绍与Python示例
有多种方法可以处理时间序列数据中的噪声。本文将介绍一种在我们的研究项目中表现良好的方法,特别适用于时间序列概况中数据点较少的情况。
48 1
数据稀缺条件下的时间序列微分:符号回归(Symbolic Regression)方法介绍与Python示例
|
2月前
|
机器学习/深度学习 测试技术 数据处理
KAN专家混合模型在高性能时间序列预测中的应用:RMoK模型架构探析与Python代码实验
Kolmogorov-Arnold网络(KAN)作为一种多层感知器(MLP)的替代方案,为深度学习领域带来新可能。尽管初期测试显示KAN在时间序列预测中的表现不佳,近期提出的可逆KAN混合模型(RMoK)显著提升了其性能。RMoK结合了Wav-KAN、JacobiKAN和TaylorKAN等多种专家层,通过门控网络动态选择最适合的专家层,从而灵活应对各种时间序列模式。实验结果显示,RMoK在多个数据集上表现出色,尤其是在长期预测任务中。未来研究将进一步探索RMoK在不同领域的应用潜力及其与其他先进技术的结合。
86 4
|
9天前
|
设计模式 开发者 Python
Python编程中的设计模式:工厂方法模式###
本文深入浅出地探讨了Python编程中的一种重要设计模式——工厂方法模式。通过具体案例和代码示例,我们将了解工厂方法模式的定义、应用场景、实现步骤以及其优势与潜在缺点。无论你是Python新手还是有经验的开发者,都能从本文中获得关于如何在实际项目中有效应用工厂方法模式的启发。 ###