Python的GDAL库绘制多波段、长时序遥感影像时间曲线图

申建兵

  本文介绍基于Python中的gdal模块，对大量长时间序列的栅格遥感影像文件，绘制其每一个波段中、若干随机指定的像元的时间序列曲线图的方法。
  在之前的文章中，我们就已经介绍过基于gdal模块，对大量多时相栅格图像，批量绘制像元时间序列折线图的方法。不过当时文章中的需求，每1个时相都对应着3个不同的遥感影像文件，而每1个遥感影像文件则都仅仅只有1个波段；而在本文中，我们每1景遥感影像都对应着2个波段，我们最终绘制的多条曲线图，也都来自于这每1景遥感影像的不同波段。
  我们再来明确一下本文的需求。现在有一个文件夹，其中放置了大量的遥感影像文件，如下图所示。其中，所有遥感影像都是同一地区、不同成像时间的图像，其各自的空间参考信息、像元行数与列数等都是一致的，文件名中有表示成像日期的具体字段；且每1景遥感影像都具有2个波段。现在我们希望，在遥感影像覆盖的区域内，随机选取若干的像元，基于这些像元，我们绘制其随时间变化的曲线图。因为我们的每个遥感影像都有2个波段，且都希望绘制出曲线图，因此最终的曲线图一共就有2条曲线。

  明确了需求，我们就可以开始代码的撰写。本文用到的代码如下。

1. # -*- coding: utf-8 -*-
2. """
3. Created on Tue Jul 25 23:04:41 2023
5. @author: fkxxgis
6. """
8. import os
9. import random
10. import matplotlib.pyplot as plt
11. from osgeo import gdal
13. def load_image(image_path):
14.     dataset = gdal.Open(image_path)
15.     band1 = dataset.GetRasterBand(1).ReadAsArray()
16.     band2 = dataset.GetRasterBand(2).ReadAsArray()
17.     del dataset
18.     return band1, band2
20. def plot_time_series(image_folder, pic_folder, num_pixels):
21.     image_files = [file for file in os.listdir(image_folder) if file.endswith(".tif")]
22.     band1_merge, band2_merge = [], []
23.     i = 0
24.    
25.     for image_file in image_files:
26.         image_path = os.path.join(image_folder, image_file)
27.         band1, band2 = load_image(image_path)
28.         band1_merge.append(band1)
29.         band2_merge.append(band2)
30.         i += 1
32.     x_size, y_size = band1.shape
33.     pixel_indices = random.sample(range(x_size * y_size), num_pixels)
35.     for pixel_index in pixel_indices:
36.         x, y = divmod(pixel_index, y_size)
37.         band_list_1, band_list_2 = [], []
38.         for i in range(len(band1_merge)):
39.             band_data_1 = band1_merge[i]
40.             band_list_1.append(band_data_1[x, y])
41.             band_data_2 = band2_merge[i]
42.             band_list_2.append(band_data_2[x, y])
44.         plt.figure()
45.         plt.plot(range(len(band1_merge)), band_list_1, label="Band 1")
46.         plt.plot(range(len(band1_merge)), band_list_2, label="Band 2")
47.         plt.xlabel("Time")
48.         plt.ylabel("NDVI")
49.         plt.ylim(0, 1000)
50.         plt.title(f"Time Series for Pixel at ({x}, {y})")
51.         plt.legend()
52.         plt.savefig(os.path.join(pic_folder, str(x) + "_" + str(y)))
53.         plt.show()
55. image_folder_path = "E:/02_Project/Result/test"
56. pic_folder_path = "E:/02_Project/TIFF/Plot"
57. num_pixels = 50
58. plot_time_series(image_folder_path, pic_folder_path, num_pixels)

复制代码

  上述代码的具体含义如下。
  首先，我们导入了需要使用的库；其中，os用于处理文件路径和目录操作，random用于随机选择像素，matplotlib.pyplot则用于绘制图像。
  随后，我们定义函数load_image(image_path)；这个函数接收一个影像文件路径image_path作为输入参数。随后，在函数内使用gdal库打开该影像文件，然后提取其第一个和第二个波段的数据，并分别存储在band1和band2中。最后，函数返回这两个波段的数据。
  接下来，我们定义函数plot_time_series(image_folder, pic_folder, num_pixels)；这个函数接收三个输入参数，分别为image_folder、pic_folder和num_pixels。其中，image_folder为包含多个.tif格式的影像文件的文件夹路径，pic_folder是保存生成的时间序列图像的文件夹路径，而num_pixels则指定了随机选择的像素数量，用于绘制时间序列图——这个参数设置为几，我们最后就会得到几张结果图像。
  在这个函数的内部，我们通过os.listdir函数获取image_folder中所有以.tif结尾的影像文件，并将这些文件名存储在image_files列表中。然后，我们创建两个空列表band1_merge和band2_merge，用于存储所有影像文件的2个波段数据。接下来，我们遍历所有影像文件，逐个加载每个影像文件的全部波段数据，并将它们添加到对应的列表中。其次，使用random.sample函数从像素索引的范围中随机选择num_pixels个像素的索引，并保存在pixel_indices列表中。接下来，我们遍历并恢复pixel_indices中的每个像素索引，计算该像素在每个影像中的每个波段的时间序列数据，并存储在band_list_1、band_list_2列表中。
  随后，我们即可绘制两个时间序列图，分别表示2个波段在不同影像日期上的数值。最后，我们将图像保存到指定的文件夹pic_folder中，命名规则为x_y，其中x与y分别代表像素的横、纵坐标。
  执行上述代码，我们即可在指定的文件夹路径下看到我们生成的多张曲线图；如下图所示。

  其中，每1张图像都表示了我们2个波段在这段时间内的数值走势；如下图所示。

  至此，大功告成。

来源:https://www.cnblogs.com/fkxxgis/p/18340048
免责声明：由于采集信息均来自互联网，如果侵犯了您的权益，请联系我们【E-Mail:cb@itdo.tech】 我们会及时删除侵权内容，谢谢合作！