百度360必应搜狗淘宝本站头条
当前位置:网站首页 > 技术资源 > 正文

Python GDAL绘制栅格图像的时间序列折线图

off999 2024-10-31 14:04 19 浏览 0 评论

??本文介绍基于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

相关推荐

推荐一款Python的GUI可视化工具(python 可视化工具)

在Python基础语法学习完成后,进一步开发应用界面时,就需要涉及到GUI了,GUI全称是图形用户界面(GraphicalUserInterface,又称图形用户接口),采用图形方式显示的计算机操...

教你用Python绘制谷歌浏览器的3种图标

前两天在浏览matplotlib官方网站时,笔者无意中看到一个挺有意思的图片,就是用matplotlib制作的火狐浏览器的logo,也就是下面这个东东(网页地址是https://matplotlib....

小白学Python笔记:第二章 Python安装

Windows操作系统的python安装:Python提供Windows、Linux/UNIX、macOS及其他操作系统的安装包版本,结合自己的使用情况,此处仅记录windows操作系统的python...

Python程序开发之简单小程序实例(9)利用Canvas绘制图形和文字

Python程序开发之简单小程序实例(9)利用Canvas绘制图形和文字一、项目功能利用Tkinter组件中的Canvas绘制图形和文字。二、项目分析要在窗体中绘制图形和文字,需先导入Tkinter组...

一文吃透Python虚拟环境(python虚拟环境安装和配置)

摘要在Python开发中,虚拟环境是一种重要的工具,用于隔离不同项目的依赖关系和环境配置。本文将基于windows平台介绍四种常用的Python虚拟环境创建工具:venv、virtualenv、pip...

小白也可以玩的Python爬虫库,收藏一下

最近,微软开源了一个项目叫「playwright-python」,作为一个兴起项目,出现后受到了大家热烈的欢迎,那它到底是什么样的存在呢?今天为你介绍一下这个传说中的小白神器。Playwright是...

python环境安装+配置教程(python安装后怎么配置环境变量)

安装python双击以下软件:弹出一下窗口需选择一些特定的选项默认选项不需要更改,点击next勾选以上选项,点击install进度条安装完毕即可。到以下界面,证明安装成功。接下来安装库文件返回电脑桌面...

colorama,一个超好用的 Python 库!

大家好,今天为大家分享一个超好用的Python库-colorama。Github地址:https://github.com/tartley/coloramaPythoncolorama库是一...

python制作仪表盘图(python绘制仪表盘)

今天教大家用pyecharts画仪表盘仪表盘(Gauge)是一种拟物化的图表,刻度表示度量,指针表示维度,指针角度表示数值。仪表盘图表就像汽车的速度表一样,有一个圆形的表盘及相应的刻度,有一个指针...

总结90条写Python程序的建议(python写作)

  1.首先  建议1、理解Pythonic概念—-详见Python中的《Python之禅》  建议2、编写Pythonic代码  (1)避免不规范代码,比如只用大小写区分变量、使用容易...

[oeasy]python0137_相加运算_python之禅_import_this_显式转化

变量类型相加运算回忆上次内容上次讲了是从键盘输入变量input函数可以有提示字符串需要有具体的变量接收输入的字符串输入单个变量没有问题但是输入两个变量之后一相加就非常离谱添加图片注释,不超过1...

Python入门学习记录之一:变量(python中变量的规则)

写这个,主要是对自己学习python知识的一个总结,也是加深自己的印象。变量(英文:variable),也叫标识符。在python中,变量的命名规则有以下三点:>变量名只能包含字母、数字和下划线...

掌握Python的"魔法":特殊方法与属性完全指南

在Python的世界里,以双下划线开头和结尾的"魔法成员"(如__init__、__str__)是面向对象编程的核心。它们赋予开发者定制类行为的超能力,让自定义对象像内置类型一样优雅工...

11个Python技巧 不Pythonic 实用大于纯粹

虽然Python有一套强大的设计哲学(体现在“Python之禅”中),但总有一些情况需要我们“打破规则”来解决特定问题。这触及了Python哲学中一个非常核心的理念:“实用主义胜于纯粹主义”...

Python 从入门到精通 第三课 诗意的Python之禅

导言:Python之禅,英文名是TheZenOfPython。最早由TimPeters在Python邮件列表中发表,它包含了影响Python编程语言设计的20条软件编写原则。它作为复活节彩蛋...

取消回复欢迎 发表评论: