我有一个DEM,我想对其进行平滑或概括以消除地形的极端情况(切掉峰和填谷)。理想情况下,我还希望控制“模糊性”的半径或水平。最后,我将需要一组从稍微模糊到真正模糊的栅格。 (理论上,最模糊的是所有值的算术平均值的恒定栅格。)

我是否可以使用任何工具或方法(基于Esri,GDAL和GRASS)?

我需要在家烘烤自己的高斯模糊例程吗?

我可以使用低通滤镜(例如ArcGIS的滤镜),如果是这样,我是否需要运行一遍才能获得大半径效果? br />

评论

仅将栅格导出到更大的像元大小该怎么办?这还会导致极端的静音吗?

是的,这也将减少极端情况(假设隐式重采样涉及某种形式的平均),但这是平滑DEM的一种糟糕方法:您将创建少量大块。顺便说一句,通常不需要导出栅格来执行此操作;聚合以及重新采样到不同的像元大小是通常在基于光栅的软件中发现的基本操作。

#1 楼

高斯模糊只是加权焦距平均值。您可以使用短距离圆形邻域序列(未加权)来高精度地重新创建它:这是中心极限定理的应用。

您有很多选择。 “过滤器”太有限-仅适用于3 x 3邻域-因此请不要打扰。大型DEM的最佳选择是将计算在ArcGIS之外带入使用快速傅里叶变换的环境中:它们进行相同的焦点计算,但是(相比之下)它们执行起来非常快。 (GRASS具有FFT模块。该模块用于图像处理,但如果可以将其以合理的精度重新缩放到0..255范围,则可以将其按入DEM的服务。)除非至少有两个解决方案值得考虑的是:



创建一组邻域权重,以近似高斯模糊,以获得较大的邻域。使用此模糊的连续遍历来创建更平滑的DEM序列。

(权重以exp(-d ^ 2 /(2r)计算),其中d是距离(如果愿意,可以在单元格中显示),r是有效半径(也可以在单元格中)。必须在至少延伸到3r的圆内计算。这样做之后,将每个权重除以它们的总和,这样最后它们的总和就是1。 ;只需反复执行循环焦点均值。我正是为了研究派生网格(如坡度和坡向)如何随DEM的分辨率变化而做的。
两种方法都可以很好地工作,并且在经过前几遍之后,几乎没有什么选择,但是收益递减:n个连续聚焦装置(均使用相同的邻域大小)的有效半径仅为n乘以焦点均值半径的平方根。因此,对于大量的模糊,您将要从大半径的邻域重新开始。如果使用未加权的均值,则在DEM上运行5-6次。如果您使用的权重接近高斯,则只需要执行一次即可:但是您必须创建权重矩阵。

这种方法确实将DEM的算术平均值作为极限值。 >

评论


如果您的数据出现峰值,您可以先尝试使用中值过滤器(en.wikipedia.org/wiki/Median_filter),然后再根据Whuber的建议应用更一般的模糊处理。

– MerseyViking
2011年5月9日在8:44

@Mersey这是一个很好的建议。我从未见过具有局部异常值的DEM,但是我再也不必处理原始DEM(例如原始LIDAR结果)。您无法使用FFT进行中值滤波,但是(通常)仅需要3 x 3的邻域,因此无论如何这都是一种快速的操作。

– hu
2011年5月9日14:49

谢谢胡伯。我必须承认,我只使用过预处理的LiDAR数据,但是SRTM数据中存在一些明显的峰值,这将受益于中值滤波器。但是,它们的确倾向于2或3个样本宽,因此需要一个更大的中值滤波器。

– MerseyViking
2011年5月9日15:19

@Mersey您仍然可以使用较大的5 x 5或7 x 7的中值滤镜。但是,如果您打算(例如)使用101 x 101滤镜,则请等待!您还提出了一个值得阐述的重要观点:在做任何事情之前对DEM进行探索性分析是一个非常好的主意。这将包括识别峰值(局部异常值)并表征其大小和程度。在使用过滤器清除它们之前,您需要确保它们确实是伪影(而不是某些实际现象)!

– hu
2011年5月9日17:12

+1用于高程数据的FFT。实际上,我已经针对32位NED数据在草丛中进行了这项工作,以消除双向条带化。最后,这也是有问题的,因为它重新引入了困扰许多其他轮廓派生DEM的梯田效果。

–杰伊·瓜纳里(Jay Guarneri)
13年5月1日在21:45

#2 楼

我一直在探索SciPy的signal.convolve方法(基于本食谱),并通过以下代码段取得了一些非常不错的成功: import numpy as np from scipy.signal import fftconvolve def gaussian_blur(in_array, size): # expand in_array to fit edge of kernel padded_array = np.pad(in_array, size, 'symmetric') # build kernel x, y = np.mgrid[-size:size + 1, -size:size + 1] g = np.exp(-(x**2 / float(size) + y**2 / float(size))) g = (g / g.sum()).astype(in_array.dtype) # do the Gaussian blur return fftconvolve(padded_array, g, mode='valid')

我在另一个函数中使用了此函数,该函数通过GDAL读取/写入float32 GeoTIFF(无需重新缩放为0-255字节进行图像处理),并且我一直在尝试使用像素大小(例如2、5、20),并且输出效果非常好(在ArcGIS中以1:1像素和恒定的最小/最大范围可视化):
注意:此答案已更新为使用更快的基于FFT的信号。fft卷积处理功能。

评论


+1不错的解决方案!我不确定,但是signal.convolve使用FFT是一个很好的选择。

– hu
2011年6月14日在20:28

我一直在寻找一些我正在编写的自动缝合工具的模糊代码,偶然发现了这一点。干得好@MikeToews!

–拉吉·亚瑟(Ragi Yaser Burhum)
13年1月15日在18:07

@RagiYaserBurhum很想听听更多关于您的工具的信息。 MikeToews很棒的答案和非常感谢的代码片段。

–杰伊·劳拉(Jay Laura)
13年1月16日在2:55

@JayLaura没什么特别的,只是编写一个工具来自动缝制我和一些朋友拿着气球拍摄的图像。使用Orfeo工具箱类orfeo-toolbox.org/SoftwareGuide/…

–拉吉·亚瑟(Ragi Yaser Burhum)
13年1月16日下午5:37

@whuber在修改此例程后,并未使用FFT,但现在使用FFT,并且速度更快。

– Mike T
2014年5月12日在1:28

#3 楼

如果时间不是太长且太复杂,这可能是对MikeT出色答案的评论。我已经玩了很多,并根据他的功能制作了一个名为FFT卷积滤波器(尚未进入“实验”阶段)的QGIS插件。除了平滑外,该插件还可以通过从原始像素中减去平滑后的栅格来锐化边缘。

我在此过程中对Mike的功能进行了一些升级: />
有效性检查是不言而喻的,但重要的是强制转换为浮动并返回。在此之前,由于将值除以值的总和(g / g.sum()),该函数使整数数组变为黑色(仅零)。

#4 楼

在QGIS中,使用Orfeo Toolbox图像过滤可以轻松获得良好的结果。这是合理的快速,批处理模式可以正常工作。可以使用高斯,均值或各向异性扩散。

请注意,Radius是指像元数,而不是距离。

下面是使用平滑(高斯)的示例:
原始文件:









#5 楼

高斯模糊和炫酷动画的不错解决方案。关于上述Esri过滤器工具,基本上就是硬编码为3x3大小的Esri“焦点统计”工具。焦点统计工具为您提供了更多有关移动滤镜的形状,大小和要运行的统计信息的选项。
http://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/focal-statistics.htm

您还可以制作“不规则”过滤器您可以在其中传递自己的文本文件,其中包含要用于每个单元格的权重。文本文件在过滤器区域中具有所需的任意多行,各列以空格分隔。我猜您应该始终使用奇数行和列,因此目标单元格位于中间。

我创建了一个excel电子表格,可以使用不同的权重进行播放,然后将其复制/粘贴到此文件中。如果您调整公式,它将获得与以上相同的结果。