Python忽略NoData计算多张遥感影像的像元平均值：whitebox库

偏执的太偏执、 2024-03-22 22:37 88阅读 0赞

在文章[Python ArcPy批量计算多时相遥感影像的各像元平均值][Python ArcPy]中，我们介绍了基于**Python**中**Arcpy**模块实现多时相遥感影像数据的平均值求取方法。但是这一方法具有一个问题，即对于任意一个像元，只要该像元在**任意一个时相**的图像中是无效值（即为**NoData**），那么该像元在最终求出的平均值**结果图中**像素值也将会是无效值**NoData**。这就导致在我们最终计算得到的平均值结果图层中，具有很多空白区域（像素值为**NoData**的区域）。

为了解决这一问题，这里我们再介绍一种基于**Python**中另一个地理空间数据分析库——`whitebox`，实现多时像遥感影像数据逐像元平均值的求取方法。

首先，需要下载并安装`whitebox`这一模块。如果大家电脑中已经有了**Anaconda**环境，就可以直接按照[Python地理分析库whitebox在Anaconda中的配置][Python_whitebox_Anaconda]这篇文章中介绍的方法下载、安装`whitebox`。

本文要实现的需求和文章[Python ArcPy批量计算多时相遥感影像的各像元平均值][Python ArcPy]中的一致，这里就不再赘述。本文所需用到的代码如下。

# -*- coding: utf-8 -*-
    """
    Created on Sun Apr 17 15:04:29 2022
    
    @author: fkxxgis
    """
    
    import glob
    from whitebox import WhiteboxTools
    
    tif_file_path="E:/LST/Data/MODIS/test/"
    average_file_path="E:/LST/Data/MODIS/06_Average/"
    
    wbt=WhiteboxTools()
    wbt.work_dir=tif_file_path
    
    tif_file_name=glob.glob(tif_file_path+"*.tif")
    tif_file_year=tif_file_name[0][-18:-14]
    one_year_tif_list=[]
    
    for tif_file in tif_file_name:
        if tif_file[-18:-14]==tif_file_year:
            one_year_tif_list.append(tif_file)
            tif_file_temp=tif_file
            if tif_file==tif_file_name[len(tif_file_name)-1]:
                wbt.average_overlay(inputs=';'.join(one_year_tif_list),
                                    output=average_file_path+tif_file_year+"_Ave.tif")
        else:
            wbt.average_overlay(inputs=';'.join(one_year_tif_list),
                                output=average_file_path+tif_file_year+"_Ave.tif")
            one_year_tif_list=[]
            one_year_tif_list.append(tif_file)
            tif_file_year=tif_file[-18:-14]

其中，`tif_file_path`是原有计算平均值前遥感图像的保存路径，`average_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`年成像的`8`幅遥感影像都已经遍历过了，那么就对这`8`景遥感影像加以逐像元的平均值求取，并开始对下一个年份（即`2005`年）成像的遥感影像继续加以计算；如果还没有读取完毕，例如假如`2001`年成像的`8`幅遥感影像目前仅遍历到了第`5`幅，那么就不求平均值，继续往下遍历，直到遍历完`2001`年成像的`8`幅遥感影像。

这里相信大家也看到了为什么我们要在前期先将文件夹中的文件按照“**名称**”排序——是为了保证**同一年成像的所有遥感影像都排列在一起**，遍历时只要遇到一个**新的年份**，程序就知道**上一个年份**的所有图像都已经遍历完毕了，就可以将**上一个年份**的所有栅格图像加以平均值求取。

本文代码和前期博客中代码不一样的部分就在于，这里是用到`whitebox`模块而非`arcpy`模块来实现同一年份遥感影像的逐像元平均值求取。在这里，`wbt.average_overlay()`函数就是我们实现这一步骤的关键，其中`inputs`参数表示需要进行平均值计算的同一年份的所有遥感影像，`output`表示求取平均值后得到的结果图像。

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

这里需要注意，由于我们在此没有用到`arcpy`模块，因此代码也就不一定非要在 **IDLE (Python GUI)** 中运行了，常见的编译器都可以运行。在代码运行过程中，还可以看到具体运行情况与进度。

![8eca7307048b7aaf5161033bd6fd0d6d.png][]

代码运行完毕后，即可得到求解平均值后的结果图层。

最后还有一个问题——在我用这一代码进行实践后发现，如果**计算平均值前的图层**具有**两个或两个以上**的波段，那么得到的结果图层整体看还好，如下图所示。

![4757ba176e19fc3c2f621d7b5a9cfa75.png][]

但放大后，会发现得到的结果呈现出如下所示的条带状。

![798ca752510e8ddfef6c53dd3d07459e.png][]

而如果**计算平均值前的图层**仅具有**一个**波段的话，就不会出现这种问题；如下图所示。

![f801eb194b5a1dc2981d02b7459bf518.png][]

因此，大家在使用本文的代码对大量**长时间序列**栅格遥感影像的**每一个像元**进行**忽略Nodata值**的**多时序平均值**求取时，一定注意输入图层要**仅含有一个波段**；否则结果就会出现条带状的错误。

[Python ArcPy]: https://www.cnblogs.com/fkxxgis/p/17330411.html
[Python_whitebox_Anaconda]: https://www.cnblogs.com/fkxxgis/p/17227352.html
[8eca7307048b7aaf5161033bd6fd0d6d.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/14/e8d63ba57c2540e897a52ef5b4874db0.png
[4757ba176e19fc3c2f621d7b5a9cfa75.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/14/31858f5507f54948822f891b6df60886.png
[798ca752510e8ddfef6c53dd3d07459e.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/14/222cff63cc8340c2bb6607d4852dff9a.png
[f801eb194b5a1dc2981d02b7459bf518.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/14/f59ecb460ede46c7a5ad574001064272.png