Python批量绘制多张遥感影像并分别设定子图标题

  本文介绍基于Python语言的matplotlib模块与gdal模块，读取大量长时间序列遥感影像，分别将其不同时相的图像作为子图，绘制在1个完整的大图中，并分别为每1个子图构建、显示标题的方法。

  首先，我们明确一下本文的需求。现有一个文件夹，其中包含大量.tif格式的遥感影像文件，如下图所示；其中，每1景遥感影像的文件名称中都含有其成像时间，且每1景图像成像时间都相差8天——从2021年的第001天开始，到第361天结束，一共是46景图像。

  我们希望实现的是，绘制每1景遥感影像，并将绘制得到的46景图像放在1个大图中，并且标注每1景图像具体的成像时间；最终的预期结果如下图所示。其中，下图仅展示了46景图像比较靠后的几景图像，大家理解意思即可。

  此外，从上图还可以看到，在结果图的末尾处，还有一个图例，表示图像中像元的数值。

  明确了需求，接下来即可开始代码的撰写。本文所用代码如下。

# -*- coding: utf-8 -*-
"""
Created on Wed Mar 13 21:04:48 2024

@author: fkxxgis
"""

import os
import fnmatch
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
from matplotlib.colors import LinearSegmentedColormap

file_path_tif = r"E:\04_Reconstruction\96_TIFF\GF_Original"
file_path_plot = r"E:\04_Reconstruction\96_TIFF"
file_extension = "*.tif"
band_idx = 4;
row = 12
col = 4

file_all = os.listdir(file_path_tif)
file_all_tif = [file_tif for file_tif in file_all if fnmatch.fnmatch(file_tif, file_extension)]
file_all_tif.sort()

if band_idx == 1:
    cmap_colors = [(0, "white"), (0.1, "powderblue"), (0.2, "lightblue"), (0.3, "lightskyblue"), (0.4, "skyblue"),
                   (0.5, "deepskyblue"), (0.6, "cornflowerblue"), (0.7, "royalblue"), (0.8, "blue"),
                   (0.9, "mediumblue"), (1, "darkblue")]
    cmap = LinearSegmentedColormap.from_list("", cmap_colors)
    file_pic_name = "blue_ori.png"
elif band_idx == 2:
    cmap_colors = [(0, "white"), (0.1, "palegreen"), (0.2, "lightgreen"), (0.3, "greenyellow"), (0.4, "lawngreen"),
                   (0.5, "springgreen"), (0.6, "lime"), (0.7, "limegreen"), (0.8, "forestgreen"),
                   (0.9, "green"), (1, "darkgreen")]
    cmap = LinearSegmentedColormap.from_list("", cmap_colors)
    file_pic_name = "green_ori.png"
elif band_idx == 3:
    cmap_colors = [(0, "white"), (0.1, "darksalmon"), (0.2, "salmon"), (0.3, "tomato"), (0.4, "orangered"),
                   (0.5, "red"), (0.6, "indianred"), (0.7, "firebrick"), (0.8, "brown"),
                   (0.9, "maroon"), (1, "darkred")]
    cmap = LinearSegmentedColormap.from_list("", cmap_colors)
    file_pic_name = "red_ori.png"
elif band_idx == 4:
    cmap_colors = [(0, "white"), (0.1, "ivory"), (0.2, "lemonchiffon"), (0.3, "papayawhip"), (0.4, "navajowhite"),
                   (0.5, "sandybrown"), (0.6, "orange"), (0.7, "darkorange"), (0.8, "peru"),
                   (0.9, "chocolate"), (1, "saddlebrown")]
    # cmap = LinearSegmentedColormap.from_list("", cmap_colors)
    cmap = "YlOrRd"
    file_pic_name = "NIR_ori_2.png"
else:
    cmap = "Greens"
    file_pic_name = "NDVI_MODIS.png"

idx = 1
plt.rc("font", family = "Times New Roman")
fig = plt.figure(figsize = (10, 40))

for file_tif in file_all_tif:
    dataset = gdal.Open(os.path.join(file_path_tif, file_tif))
    band = dataset.GetRasterBand(band_idx)
    array_band = band.ReadAsArray()
    array_band = np.nan_to_num(array_band, nan=0)
    if np.all(array_band == 0):
        img = plt.imshow(np.zeros((1, 1)), cmap = cmap)
        plt.axis("off")
        plt.title("2021" + str(idx * 8 - 7).zfill(3), y = -0.2)
        idx += 1
        continue
    array_band[(array_band > 10000) | (array_band < 0)] = 0
    array_band = array_band / 10000
    dataset = None
    
    ax = plt.subplot(row, col, idx)
    img = plt.imshow(array_band, cmap = cmap)
    plt.axis("off")
    plt.title("2021" + str(idx * 8 - 7).zfill(3), y = -0.2)
    print(idx, "finished!")
    idx += 1

cax = plt.axes([0.66, 0.15, 0.24, 0.015])
colorbar = plt.colorbar(cax=cax, orientation="horizontal")
colorbar.outline.set_visible(False)
plt.savefig(file_path_plot + "/" + file_pic_name, dpi = 300, bbox_inches = "tight", pad_inches = 0)
plt.show()

  其中，代码的具体含义如下。

  首先，我们导入所需的模块和库。其中，os用于处理文件和路径操作，fnmatch用于文件名匹配，numpy用于处理数组数据，matplotlib.pyplot用于生成图像和子图，osgeo.gdal用于读取遥感图像数据，matplotlib.colors.LinearSegmentedColormap则用于定义自定义颜色映射。

  其次，我们定义一些变量。其中，file_path_tif为包含待绘图的遥感图像文件的路径，file_path_plot为生成的图像文件保存的路径，file_extension为待处理的图像文件的扩展名，band_idx为要处理的图像数据的波段索引（因为我这里每次只绘制1个波段，所以增加了这个参数）；而row和col则分别为生成的子图的行数和列数。

  接下来，我们获取符合条件的图像文件。其中，首先使用os.listdir()函数获取文件路径中的所有文件；随后，使用fnmatch()函数根据扩展名筛选出符合条件的文件，并对文件名按照首字母升序的顺序进行排序。

  随后，我们根据波段索引选择颜色映射。在我这里，是希望对具有蓝光、绿光、红光与近红外等4个波段的遥感影像加以绘制，每次绘制其中的1个波段；因此，对于不同的波段索引，定义不同的颜色映射；其中，对于不同波段，我们创建LinearSegmentedColormap对象，或使用预定义的颜色映射。其中，我们重点介绍一下本部分中的如下代码。

    cmap_colors = [(0, "white"), (0.1, "powderblue"), (0.2, "lightblue"), (0.3, "lightskyblue"), (0.4, "skyblue"),
                   (0.5, "deepskyblue"), (0.6, "cornflowerblue"), (0.7, "royalblue"), (0.8, "blue"),
                   (0.9, "mediumblue"), (1, "darkblue")]
    cmap = LinearSegmentedColormap.from_list("", cmap_colors)

  其中，我们定义了一个颜色映射，用于将图像数据的像素值映射到具体的颜色。cmap_colors是一个包含颜色映射信息的列表，其每个元素都是一个元组，包含两个值：第1个值是图像数据的归一化像素值（范围为0到1之间），第2个值是与归一化像素值对应的颜色。例如，(0, "white")表示归一化像素值为0时对应的颜色是白色，(0.1, "powderblue")表示归一化像素值为0.1时对应的颜色是粉蓝色。而随后的这行代码，使用LinearSegmentedColormap类从颜色映射列表cmap_colors创建一个颜色映射对象cmap。LinearSegmentedColormap是matplotlib.colors模块中的一个类，用于创建线性分段的颜色映射；其接受一个颜色映射列表作为输入，并根据列表中定义的归一化像素值和对应颜色的关系来生成颜色映射。

  回到前述完整的代码中。随后，即可开始循环处理每个图像文件。首先，打开图像文件并读取指定波段的数据；同时，对数据进行预处理，将无效值（NoData值）替换为0，并进行范围缩放——将异常数据都改为0，并将非异常数据缩小10000倍。接下来，创建子图并显示图像数据，按照图像编号设置成像时间并作为该子图的title，同时打印处理进度。

  最后，我们即可添加颜色条并保存图像。首先，在图像中添加颜色条；其次，保存生成的图像文件，并显示生成的图像。

  运行上述代码，可以得到如下所示的图片结果。

  图片比较模糊，但是可以看到，46景图像已经按照成像时间的顺序加以可视化了，且配置了图例。

  至此，大功告成。