ArcPy学习心得系列（5）遥感数据中值滤波与均值滤波实现方法(不计入NoDATA)

在数据处理与应用的过程中，我们难免会遇到一些低质量的遥感数据，低质量遥感数据一般是由于天气因素导致的，在云量较多时，卫星传感器所采集到的地面信息被云层所遮挡，导致遥感图像成像过程中产生了较多噪声，对遥感数据的精确度造成了一定影响，如果不解决这些数据中的噪声，会导致我们从遥感数据中所提取到的信息也不够精确。

为更好的处理遥感数据中的噪声，我们一般会使用均值滤波或中值滤波对遥感图像进行去噪，去除遥感图像噪声的目的在于最大程度上还原原始数据，以便于我们更好的提取遥感图像中的有效信息。

那么，在进行滤波代码的实现以前，让我们一起来看一下滤波的实现原理吧！

 

其中，滤波算法的主要思路为，设置一个n*n(n=2x+1,x>0)大小的窗口,随后对这个窗口的中心点像元进行重新赋值，以5*5的窗口以及均值滤波算法为例，在上图中，窗口A的中心像元(0,0)的值为19，而经过了均值滤波后，也就是在窗口B中，原窗口A的中心像元(0,0)的值变成了11，而这个值，是由原窗口A中所有像元的值加起来的总和除以25得到的。

通过上述算法，即可对一个窗口内的遥感数据进行滤波，基于上述思路，对整个遥感图像进行遍历，以窗口的中心像元为基准，取到N个相同大小的窗口(窗口大小会随应用场景的变化而变化)，对其中的每一个窗口中心像元进行滤波，以下，为滤波过程的核心代码，其中，self.w_1是人为设置的窗口大小，height和width分别表示遥感数据的行数和列数，result表示滤波后的新数据。

edge = int((self.w_l - 1) / 2)#计算遍历的起始位置与终点
height = data.shape[0]#获取数据的行数
width = data.shape[1]#获取数据的列数
for i in range(edge,height-edge):#设定遍历起始行与终点行
    for j in range(edge,width-edge):#设定遍历起始列与终点列
        newdata =data[i - edge:i + edge + 1, j - edge:j + edge + 1]#取一个n*n的窗口
        if self.no_data_value in newdata:
            continue
        else:
            result[i, j] = np.median(newdata)#调用np.median求中值函数

基于以上代码，即可完成中值滤波(不对NoDATA及数组边缘进行计算)。

除了滤波部分的代码以外，我们还需了解，ArcPy是如何对遥感数据进行读取以及对计算结果进行输出的，首先，我们需要利用ArcPy的Raster方法对数据进行读取，随后利用ArcPy的RasterToNumPyArray方法，将数据转换为array数组，并利用以上滤波算法对数组进行滤波。在滤波完成后，我们还需要将滤波后的结果保存为TIFF文件。

以下部分为数据读取与成果输出部分的代码。

arcpy.env.outputCoordinateSystem = self.Path#设置输出坐标系
arcpy.env.overwriteOutput = True#可以覆盖数据
tif = arcpy.Raster(self.Path)#读取数据
ExtentXmin = tif.extent.XMin#取x坐标最小值
ExtentYmin = tif.extent.YMin#取y坐标最小值
lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin)#取得数据起始点坐标(左下角坐标)
Widthcellsize=tif.meanCellWidth#像元宽度
Heightcellsize=tif.meanCellHeight#像元高度
self.no_data_value = tif.noDataValue#获取数据NoDATA
#将Raster转换为Array数组
self.data = arcpy.RasterToNumPyArray(tif, nodata_to_value=self.no_data_value)
newdata = self.Filter()
if newdata != None:
    #将滤波后的array数组转换为Raster，其中需要指定其左下角坐标，像元大小及NoDATA等参数
    newRaster = arcpy.NumPyArrayToRaster(newdata,lowerLeft,Widthcellsize,Heightcellsize,value_to_nodata=self.no_data_value)
    newRaster.save(self.Save)#将Raster保存到存储目录下
    arcpy.BuildPyramids_management(self.Save)#建立金字塔

至此，本次讲解到此结束，以下为此算法的全部代码。

import arcpy
import numpy as np

class RasterFilter:
    def __init__(self,DataPath,windows_length,method,SavePath):
        self.Path = DataPath
        self.w_l = windows_length
        self.method = method
        self.Save = SavePath

    def Start(self,method="Median"):
        arcpy.env.outputCoordinateSystem = self.Path#设置输出坐标系
        arcpy.env.overwriteOutput = True#可以覆盖数据
        tif = arcpy.Raster(self.Path)#读取数据
        ExtentXmin = tif.extent.XMin#取x坐标最小值
        ExtentYmin = tif.extent.YMin#取y坐标最小值
        lowerLeft = arcpy.Point(ExtentXmin,ExtentYmin)#取得数据起始点坐标(左下角坐标)
        Widthcellsize=tif.meanCellWidth#像元宽度
        Heightcellsize=tif.meanCellHeight#像元高度
        self.no_data_value = tif.noDataValue#获取数据NoDATA
        #将Raster转换为Array
        self.data = arcpy.RasterToNumPyArray(tif, nodata_to_value=self.no_data_value)
        newdata = self.Filter()
        if newdata != None:
            #将滤波后的array数组转换为Raster，其中需要指定其左下角坐标，像元大小及NoDATA等参数
            newRaster = arcpy.NumPyArrayToRaster(newdata,lowerLeft,Widthcellsize,Heightcellsize,value_to_nodata=self.no_data_value)
            newRaster.save(self.Save)#将Raster保存到存储目录下
            arcpy.BuildPyramids_management(self.Save)#建立金字塔
    def Filter(self):
        data = self.data#获取array数据
        result = data.copy()#复制array到result
        edge = int((self.w_l - 1) / 2)#计算遍历起始位置与终点
        height = data.shape[0]#获取数据的行数
        width = data.shape[1]#获取数据的列数
        if height - 1 - edge <= edge or width - 1 - edge <= edge:#如何数据行列数小于设置的窗口大小，则返回并抛出异常
            print("The parameter k is to large.")
            return None
        if self.method == "Median":
            for i in range(edge,height-edge):#设定遍历起始行与终点行
                for j in range(edge,width-edge):#设定遍历起始列与终点列
                    newdata =data[i - edge:i + edge + 1, j - edge:j + edge + 1]#取一个n*n的窗口
                    if self.no_data_value in newdata:
                        continue
                    else:
                        result[i, j] = np.median(newdata)#调用np.median求中值函数

        elif self.method == "Average":
            for i in range(edge,height-edge):
                for j in range(edge,width-edge):
                    newdata =data[i - edge:i + edge + 1, j - edge:j + edge + 1]
                    if self.no_data_value in newdata:
                        continue
                    else:
                        result[i, j] = np.average(newdata)#调用np.average求均值函数
        return result
if __name__=='__main__':
    rFilter = RasterFilter(r"F:\NPPOFNDVI\ASTERGDEMV2DEM_fujian\ASTGTM2_N23E116_dem.tif",5,"Median",r"F:\NPPOFNDVI\ASTERGDEMV2DEM_fujian\DEM2.tif")
rFilter.Start()

以下为DEM数据滤波前后的对比图。