ArcPy批量对大量遥感影像相减做差

今天药忘吃喽~ 2024-03-17 16:18 49阅读 0赞

本文介绍基于**Python**中**ArcPy**模块，对大量栅格遥感影像文件**批量**进行**相减做差**的方法。

首先，我们来明确一下本文的具体需求。现有一个存储有多张`.tif`格式遥感影像的文件夹，其中每一个遥感影像的文件名中都包含有该图像的成像年份，且每一个遥感影像的空间范围、像元大小等都是一致的，可以直接进行栅格相减；且文件夹内除了`.tif`格式的遥感影像文件外，还具有其它格式的文件；如下图所示。

![c5129ec251b5bad8899c07eec107abe5.png][]

我们希望，对于**同一年**成像的两景遥感影像分别进行做差处理。例如，将上图中的`2001.tif`文件减去`2001_N.tif`文件，将`2005.tif`文件减去`2005_N.tif`文件，以此类推。

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

# -*- coding: utf-8 -*-
    """
    Created on Sun Apr 24 11:12:37 2022
    
    @author: fkxxgis
    """
    
    import arcpy
    
    tif_file_path="E:/LST/Data/MODIS/16_True/"
    dif_file_path="E:/LST/Data/MODIS/17_Difference/"
    arcpy.env.workspace=tif_file_path
    
    tif_file_name=arcpy.ListRasters("*","tif")
    tif_file_year=tif_file_name[0][0:4]
    one_year_tif_list=[]
    
    for tif_file in tif_file_name:
        if tif_file[0:4]==tif_file_year:
            one_year_tif_list.append(tif_file)
            if tif_file==tif_file_name[len(tif_file_name)-1]:
                arcpy.gp.Minus_sa(one_year_tif_list[0],
                                  one_year_tif_list[1],
                                  dif_file_path+tif_file_year+"_Dif.tif")
        else:
            arcpy.gp.Minus_sa(one_year_tif_list[0],
                              one_year_tif_list[1],
                              dif_file_path+tif_file_year+"_Dif.tif")
            one_year_tif_list=[]
            one_year_tif_list.append(tif_file)
            tif_file_year=tif_file[0:4]

其中，`tif_file_path`是原有计算平均值前遥感图像的保存路径，`dif_file_path`是我们新生成的求取平均值后遥感影像的保存路径，也就是结果保存路径。

在这里，和我们前期的博客[Python ArcPy批量拼接长时间序列栅格图像][Python ArcPy]类似，需要首先在**资源管理器中**，将`tif_file_path`路径下的各文件以“**名称**”排序的方式进行排序；随后，利用`arcpy.ListRasters()`函数，获取路径下原有的全部`.tif`格式的图像文件，并截取第一个文件的部分文件名，从而获取其成像时间的具体年份。

接下来，遍历`tif_file_path`路径下全部`.tif`格式图像文件。其中，我们通过一个简单的判断语句`if tif_file[0:4]==tif_file_year:`，来确定某一年的遥感影像是否已经读取完毕——如果已经读取完毕，例如假如`2001`年成像的`2`幅遥感影像都已经遍历过了，那么就对这`2`景遥感影像做差，并开始对下一个年份（即`2005`年）成像的`2`景遥感影像继续加以计算；如果还没有读取完毕，例如假如`2001`年成像的`2`幅遥感影像目前仅遍历了第`1`幅，那么就不做差，继续往下遍历，直到遍历完`2001`年成像的`2`幅遥感影像。

这里相信大家也看到了为什么我们要在前期先将文件夹中的文件按照“**名称**”排序——首先，是为了保证**同一年成像的`2`景遥感影像都排列在一起**，遍历时只要遇到一个**新的年份**，程序就知道**上一个年份**的`2`张图像都已经遍历完毕了，就可以将**上一个年份**的`2`张栅格图像加以做差；其次，是为了保证我们的**被减数**（例如`2005.tif`文件）排在**减数**（例如`2005_N.tif`文件）的前面，从而方便我们进行做差运算。

在这里，我们实现两张栅格遥感影像相减操作的函数是`arcpy.gp.Minus_sa()`函数，其第一个参数是被减数，第二个参数是减数，第三个参数是结果保存路径与名称。

最后，通过`if tif_file==tif_file_name[len(tif_file_name)-1]:`这个判断，来确认是否目前已经遍历到文件夹中的**最后一个**图像文件。如果是的话，就需要将**当前成像年份**的`2`景图像进行差值的求取，并宣告代码完成运行。

在 **IDLE (Python GUI)** 中运行代码。代码运行完毕后，我们可以看到求取差值之后的遥感影像已经存在于我们的结果保存路径中了。

![e1f201df4c0793c825a86117dae96c12.png][]

至此，大功告成。

[c5129ec251b5bad8899c07eec107abe5.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/15/a31f84e5ee5e4dc3857c339083b6b4f0.png
[Python ArcPy]: https://www.cnblogs.com/fkxxgis/p/17218043.html
[e1f201df4c0793c825a86117dae96c12.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/15/cadedb66eed544d294bb398baa7f90bd.png