本文介绍基于Python语言，针对一个含有大量遥感影像栅格文件的文件夹，从其中第2景遥感影像开始，分别用每一景影像减去其前一景影像，从而求取二者的差值，并将每一个所得到的差值结果保存为新的一景遥感影像文件的方法。
首先，来看一下我们具体的需求。我们现在有一个文件夹，其中存放着不同成像时间的单波段遥感影像文件（多波段遥感影像文件也可以用本文的代码，只不过就是在代码读取遥感影像数据的时候，先指定一下具体要读取的波段序号即可），其文件名就表示遥感影像的成像时间，且我们已经按照文件名，对这些遥感影像文件加以排序了；如下图所示。其中，每一景遥感影像的空间范围、地理参考信息、栅格行数与列数等都是一致的。
我们希望其中每一景遥感影像之间的差值。例如，首先用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函数，传入待处理影像的文件夹路径，即可按要求开始处理影像。
运行上述代码。首先，在运行过程中，会显示当前已经计算好的差值结果文件的文件名，如下图所示；从而方便我们获取代码的执行进度。
其次，执行代码完毕后，我们即可在结果文件夹中看到对应的多个差值结果文件；如下图所示。
我们随机打开几景原始遥感影像，及其对应的差值结果文件；如下图所示，可以看到，我们差值结果文件中的像素，就是原始遥感影像文件之间像素的差值。
至此，大功告成。
欢迎关注：疯狂学习GIS
