Python GDAL求取不同栅格文件两两之差

  本文介绍基于Python语言，针对一个含有大量遥感影像栅格文件的文件夹，从其中第2景遥感影像开始，分别用每一景影像减去其前一景影像，从而求取二者的差值，并将每一个所得到的差值结果保存为新的一景遥感影像文件的方法。

  其中，本文所需实现的需求，和我们之前的文章ArcPy对多张栅格影像批量相减求取差值非常类似；但是在上述文章中，我们是基于Python中ArcPy模块实现需求的。而在本文中，我们将通过另一个Python模块——gdal库，来实现这一需求；大家基于实际需要，选择这两篇文章中的代码即可。

  首先，来看一下我们具体的需求。我们现在有一个文件夹，其中存放着不同成像时间的单波段遥感影像文件（多波段遥感影像文件也可以用本文的代码，只不过就是在代码读取遥感影像数据的时候，先指定一下具体要读取的波段序号即可），其文件名就表示遥感影像的成像时间，且我们已经按照文件名，对这些遥感影像文件加以排序了；如下图所示。其中，每一景遥感影像的空间范围、地理参考信息、栅格行数与列数等都是一致的。

  我们希望其中每一景遥感影像之间的差值。例如，首先用2020009.tif这个文件，减去2020001.tif这个文件，得到一个差值结果文件（本文选择将这个结果图像命名为2020009_diff.tif）；随后用2020017.tif这个文件，减去2020009.tif这个文件，得到一个差值结果文件（这个结果图像命名为2020017_diff.tif）；以此类推，直到文件夹内的最后一个遥感影像文件。

  明确了需求后，我们就可以开始具体的操作。首先，本文所需用到的代码如下。

# -*- coding: utf-8 -*-
"""
Created on Sun Dec 31 23:47:43 2023

@author: fkxxgis
"""

import os
from osgeo import gdal

def subtract_images(image1_path, image2_path, output_path):
    image1 = gdal.Open(image1_path)
    band1 = image1.GetRasterBand(1)

    image2 = gdal.Open(image2_path)
    band2 = image2.GetRasterBand(1)

    data1 = band1.ReadAsArray()
    data2 = band2.ReadAsArray()

    diff = data2 - data1

    driver = gdal.GetDriverByName("GTiff")
    output = driver.Create(output_path, image1.RasterXSize, image1.RasterYSize, 1, band1.DataType)

    output_band = output.GetRasterBand(1)
    output_band.WriteArray(diff)

    output.SetGeoTransform(image1.GetGeoTransform())
    output.SetProjection(image1.GetProjection())

    image1 = None
    image2 = None
    output = None

def process_images(folder_path):
    file_names = os.listdir(folder_path)
    file_names.sort()
    
    for i in range(len(file_names) - 1):
        image1_path = os.path.join(folder_path, file_names[i])
        image2_path = os.path.join(folder_path, file_names[i + 1])

        output_name = file_names[i + 1].replace(".tif", "_diff.tif")
        output_path = os.path.join(result_path, output_name)

        subtract_images(image1_path, image2_path, output_path)
        print(output_name)

folder_path = "H:/Data_Reflectance_Rec/NDVI/NDVI_2020"
result_path = "H:/Data_Reflectance_Rec/NDVI/NDVI_2020_Dif"
process_images(folder_path)

  其中，我们首先导入所需的模块。在这里，os模块用于处理文件和文件夹路径，gdal模块则用于读取和处理遥感影像数据。

  接下来，我们定义了一个subtract_images函数，用于计算两幅影像之间的差异。这个函数简单的流程如下：首先，打开影像image1_path和image2_path，并读取两幅影像的第一个波段的数据（如果大家有多个波段需要计算，那么就可以通过循环，分别对每一个波段加以处理），随后直接计算两幅影像数据的差异；接下来，创建一个新的影像文件output_path，并将差异数据写入其中；同时，设置新影像的地理转换和投影信息。最后释放打开的影像对象。

  其次，我们定义了一个process_images函数，用于处理一个文件夹中的一系列影像，我们前面的subtract_images函数，就是在这个process_images函数中调用的。这个函数简单的流程如下：首先，获取文件夹中的文件名，并按升序进行排序；其次，遍历文件名列表，对每对相邻的影像文件进行差值计算（调用subtract_images函数）；接下来，将输出影像保存到指定的结果文件夹中，并以原始文件名为基础生成新的文件名。同时，在上述处理的过程中，打印结果文件名。

  最后，我们定义了一个文件夹路径folder_path，指定待处理的影像所在的文件夹路径；同时，定义了一个结果文件夹路径result_path，指定差值影像输出的文件夹路径。接下来，调用我们前面的process_images函数，传入待处理影像的文件夹路径，即可按要求开始处理影像。

  运行上述代码。首先，在运行过程中，会显示当前已经计算好的差值结果文件的文件名，如下图所示；从而方便我们获取代码的执行进度。

  其次，执行代码完毕后，我们即可在结果文件夹中看到对应的多个差值结果文件；如下图所示。

  我们随机打开几景原始遥感影像，及其对应的差值结果文件；如下图所示，可以看到，我们差值结果文件中的像素，就是原始遥感影像文件之间像素的差值。

  至此，大功告成。