Python中ArcPy实现忽略无效值NoData的长时间序列栅格遥感影像批量逐像元求取像素平均值

  本文介绍基于Python中whitebox模块，对大量长时间序列栅格遥感影像的每一个像元进行忽略NoData值的多时序平均值求取。

  在Python中ArcPy实现对大量长时间序列栅格遥感影像批量逐像元求取像素平均值（https://blog.csdn.net/zhebushibiaoshifu/article/details/124344520）中，我们介绍了基于Python中Arcpy模块实现多时相遥感影像数据的平均值求取方法。但是这一方法具有一个问题，即对于任意一个像元，只要该像元在任意一个时相的图像中是无效值（即为NoData），那么该像元在最终求出的平均值结果图中像素值也将会是无效值NoData。这就导致在我们最终计算得到的平均值结果图层中，具有很多空白区域（像素值为NoData的区域）。

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

  首先，需要下载并安装whitebox这一模块。如果大家电脑中已经有了Anaconda环境，就可以直接按照Anaconda下Python中whitebox模块的下载与安装方法（https://blog.csdn.net/zhebushibiaoshifu/article/details/124229346）这篇文章中介绍的方法下载、安装whitebox。

  本文要实现的需求和Python中ArcPy实现对大量长时间序列栅格遥感影像批量逐像元求取像素平均值（https://blog.csdn.net/zhebushibiaoshifu/article/details/124344520）中的一致，这里就不再赘述。本文所需用到的代码如下。

# -*- 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实现对大量长时间序列栅格遥感影像批量逐像元求取像素平均值（https://blog.csdn.net/zhebushibiaoshifu/article/details/124344520）这篇文章是非常类似的。首先，同样需要在资源管理器中，将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) 中运行了，常见的编译器都可以运行。在代码运行过程中，还可以看到具体运行情况与进度。

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

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

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

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

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

欢迎关注公众号/CSDN/知乎/微博：疯狂学习GIS