ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

ArcGIS栅格Nodata填补全指南:从焦点统计到空间插值的实用方案

ArcGIS栅格Nodata填补全指南:从焦点统计到空间插值的实用方案 做水文分析的时候DEM数据里有一小片Nodata流向计算出来直接黑了一块后面所有汇流累积量全部作废做遥感影像镶嵌两景数据重叠区之外留了一圈空值想算个NDVI都算不了还有一次做坡度分级图因为原始栅格有空洞输出结果大范围空白被甲方直接打回。这些场景应该每个做GIS的都不陌生——栅格数据里的Nodata空洞可以说是实际项目里最烦人、最绕不开的问题之一。这篇文章就围绕ArcGIS里“填补栅格空缺值Nodata”这件事展开把几种主流填补方法的原理、操作步骤、选型逻辑和踩坑经验一次性说清楚。文章适合正在做地形分析、遥感影像处理、科研出图或者项目交付时被Nodata卡住的同学不管你是刚入门的小白还是已经做了几年GIS项目的老手都能从里面找到能直接用的方案。1. 栅格里的NoData到底是什么为什么它会成为一道坎1.1 NoData不是0是“压根没测到”很多新手最容易犯的错就是把NoData当成0来处理。这两者在ArcGIS里的意义完全不同。0是一个有效测量值比如某个位置的海拔就是0米或者某个像元的地表温度就是0摄氏度而NoData表示这个位置“什么都没有”既没有测量值也没有推算值它是栅格数据里一种特殊的状态标记。在ArcGIS的浮点型栅格里NoData通常用一个极其特殊的数值来表示比如-9999但ArcGIS在读取时并不会把它当普通数值参与计算而是将其识别为“无效像元”在统计、渲染、分析时自动跳过。在整型栅格里NoData往往用255或者其他预先约定的值来标记。如果你在属性表里看到某个值被标记为“NoData”那就是这个像元被排除在所有运算之外了。这个区别带来了一个很有意思的后果如果你在栅格计算器里直接写一个dem * 2的表达式遇到NoData的位置输出仍然是NoData而0的位置会变成0。也就是说NoData会像病毒一样在计算链里传播。这就是为什么它必须被处理掉否则它会影响后续所有分析环节。1.2 空洞是怎么产生的想填空洞先搞清楚空洞从哪来。根据我的实际项目经验栅格里的Nodata区域大致有这么几个来源遥感影像的云和阴影遮挡。光学遥感影像在拍摄时被云层、云影、山体阴影遮挡传感器在这些位置没有收到有效信号记录的DN值就是无效的转成反射率或辐射亮度后自然形成NoData区域。Landsat、Sentinel-2这类影像尤其常见一景影像里5%~10%的NoData都不稀奇。传感器硬件问题。有些传感器的探测器阵列在某个波段有坏道或坏线扫描出的影像会有一条或多条NoData条带。虽然现在主流传感器在轨定标后会在L1级产品里做坏线修复但在某些原始数据或第三方处理产品里还是能看到。DEM数据里的测量盲区。雷达遥感如SRTM、ALOS在水体、冰川、陡峭峡谷等区域信号反射要么太弱要么太乱反演出来的高程就是无效值。所以你会看到很多公开DEM在水库、河流、冰川区域是整片NoData。镶嵌拼接的接缝外区域。多景影像拼到一起重叠区之外的像元会因为没有数据来源而被设为NoData。这是最常见也最让人头疼的情况因为你辛苦镶出来的影像四周总是带一圈“黑边”。重投影和坐标转换。栅格从WGS84转到CGCS2000或者其他投影坐标系时像元位置会被重新采样边缘区域因为原数据范围不够新坐标系下无法计算到有效值就会形成空洞。这个在最初的标题相关热词里有“栅格影像wgs84坐标系转cgcs2000”这种情况所以特别说明一下。1.3 哪些场景必须填哪些场景千万别填区分“必须填”和“不能填”很重要。下面我用一个表格来梳理场景类型是否建议填补原因水文分析填洼、流向、累积量必须填流向计算遇到NoData会产生逻辑中断下游结果全部出错坡度坡向、曲率等地形因子提取必须填邻域算子遇到NoData像元窗口内计算失效输出也变NoData多时相影像代数运算如NDVI必须填NoData在运算中传播结果区域会缩小甚至断裂插值分析、机器学习建模的输入数据必须填算法通常要求完整网格输入制图出图时想隐藏空洞不建议填单纯显示问题用掩膜或透明色处理即可不需要造数据研究区外的天然空缺绝对不能填把研究区外的空白区域插值补上等于制造了虚假数据精度要求极高的专业成果需要谨慎填填补值只是“合理估计”不是实测值必须记录清楚填补过程和不确定性一句话总结如果空洞影响到后续分析链条那就必须处理如果只是显示层的问题就别画蛇添足。还有一个很关键的点填补Nodata本质上是在用数学模型“猜”那些位置的数值相当于给数据注入了不确定性所以任何正经项目里都建议把填补过程写到处理报告中。2. 动手前的准备摸清空洞的分布和成因2.1 用栅格计算器快速定位Nodata在动手填补之前一定要先搞清楚空洞的分布情况。闭着眼睛就上焦点统计填了半天发现还有大面积空洞没填上这种情况我见得太多了。第一步先用栅格计算器Spatial Analyst Toolbox → Map Algebra → Raster Calculator统计一下NoData像元的分布范围表达式很简单IsNull(dem)这个表达式会生成一个新的栅格NoData区域取值1有效区域取值0。加载到地图里你可以直接目视看到空洞到底在哪儿。然后打开属性表看像是元数量统计一下空洞占整体像元的比例。如果想更精确地知道空洞像元数和面积可以用GetRasterProperties工具查看栅格的“SUMM”等属性或者直接用Zonal Histogram分区统计。这一步花两分钟做一下能帮你判断后面该用哪种填补策略。2.2 判断空洞的类型目视检查之后接下来要判断空洞形态。根据我处理过的数据空洞大致分为三类每类的处理思路完全不同孤立小洞。面积小、数量少的离散NoData像元通常三五成群直径不超过几个像元。这种情况用焦点统计或Nibble工具就能快速填好操作简单、速度快、效果好。大面积连片空洞。像云遮挡区域、水体区域、镶嵌边缘那种几十甚至上百像元宽度的空洞。这种空洞用窗口小的焦点统计填不动因为统计窗口内有效像元太少用Nibble会产生明显的块状感最稳妥的方案是空间插值。边缘条带型空缺。沿着栅格四周分布的一圈或几行NoData。这种空缺的填补策略和小洞、大洞都不同因为边缘外侧没有任何参考数据焦点统计和Nibble都会失效通常需要用边缘向外推的趋势面插值或者干脆裁剪掉这一圈。2.3 一个必须留意的处理顺序问题先投影还是先填补这个细节很容易被忽略但在重投影场景里真的能坑人。如果你拿到的是WGS84坐标系的栅格后续要做CGCS2000或者高斯投影那就要先想清楚是先重投影再填补还是先填补再重投影我的建议几乎永远是先重投影再填补。原因在于重投影的过程会做重采样像元位置、大小、值都发生了变化原来NoData区域的形状和位置也会跟着变。如果你先在原始坐标系里把空洞填好了再重投影重采样过程可能因为像元对齐问题产生新的NoData或者改变已经填好的值。反过来如果先重投影在新的坐标系里统一处理空洞那么后续分析链路上就不会再被重投影影响。当然有一种例外如果重投影时NoData区域特别大、占比特别高重采样后空洞会变得非常不规则填补难度更大。这种情况下我会先做一次粗略的填充让重投影后的数据更“完整”然后再在新坐标系里做一次精细填补。但这种情况不多见一般先投影后填补就够了。3. 最常用的三种填补方法原理与实操3.1 栅格计算器焦点统计小空洞的快速填充这个方法是处理孤立小空洞最经典、最快速的方案核心思路就是把NoData像元的值替换为其周边有效像元值的统计值通常是平均值或中位数。在栅格计算器里输入的表达式如下Con(IsNull(dem), FocalStatistics(dem, NbrRectangle(3,3,CELL), MEAN, NODATA), dem)我来拆解一下这个表达式到底在干什么。IsNull(dem)生成一个布尔掩膜标记出空洞位置FocalStatistics(dem, NbrRectangle(3,3,CELL), MEAN, NODATA)对每个像元计算其3x3邻域的平均值而且特意设置了NODATA参数表示统计时忽略邻域内的NoData像元只对有效像元求平均最后用Con()条件函数将空洞位置替换为邻域平均值有效位置保留原值。极易出错的一个点就是FocalStatistics第四个参数NODATA。如果不设置这个参数ArcGIS默认会把邻域内的NoData当作特殊值单独处理统计结果会变成NoData那整个表达式就废了。网上很多教程没提这个参数照抄下来的结果就是填不上空洞还以为工具出问题了。如果数据里有一些较大的空洞3x3窗口填不上可以改成5x5、7x7甚至更大的窗口。但要注意窗口越大填补结果越平滑细节损失越严重尤其是地形数据窗口过大会把山谷、山脊的形态磨平。这里还要多说一句迭代问题。一次ConFocalStatistics处理之后如果空洞区域很大最中心位置可能还是NoData。解决办法很简单把输出结果再次代入表达式也就是做多次迭代。比如先3x3填一遍再5x5填一遍直到IsNull的结果里没有值为1的像元为止。手动嵌套几层也可以但我更推荐用后面会讲的Python脚本做自动迭代。3.2 Nibble工具用最近的有效像元“啃”掉空洞Nibble是Spatial Analyst工具箱里的一个专用工具位置在Spatial Analyst Tools → Map Algebra → Nibble它做的事情很直观把被掩膜标记为空洞的像元替换成离它最近的有效像元的值像一个“蚕食”的过程。Nibble有两个关键输入输入栅格。需要填补的原始栅格。掩膜栅格。用来标记哪些位置需要被替换。掩膜栅格里值为1的位置是空洞要被替换值为0或NoData的位置是有效区保留原值。同时有一个常用可选参数“使用NoData值进行NibbleNibble NoData values”勾选后Nibble会忽略输入栅格中NoData像元的可用性直接把NoData像元替换为最近非NoData像元的值。这个参数在ArcMap和ArcGIS Pro的对话框里都有记得勾上。实际操作的流程一般是这样先用栅格计算器把空洞位置生成掩膜Con(IsNull(dem), 1, 0)然后运行Nibble工具将输入栅格设为原始DEM掩膜栅格设为刚才生成的掩膜勾选Nibble NoData values输出新栅格。运行完成后空洞区域会被周围最近的有效像元值填满。Nibble的优点在于速度非常快而且不会产生插值方法那种平滑过渡的效果填补结果和周围数据的纹理一致性比较强。它的缺点是空洞越大填补出的区域会越长条状看起来像一块一块的“补丁”视觉上有点怪另外它只会用最近像元的值不会考虑周围数据的变化趋势所以在大面积空洞中填充出的值可能会和周围数据衔接不自然。我的使用经验是对于小于10像元的孤立小洞Nibble非常好用对于大片云遮挡区域或者大面积水体空洞效果就不太行了建议用空间插值。3.3 空间插值法大范围空洞的平滑填补当空洞范围大到几十上百像元时焦点统计和Nibble都难以满足要求这时候最可靠的做法就是空间插值。核心思路是利用空洞周围的有效像元作为样本点通过插值算法模拟出整个区域的连续表面再把插值结果填入原来的空洞位置。具体操作流程如下第一步把有效像元转成点。使用工具箱里的Raster to Point工具输入原始栅格输出点要素类。转换过程中NoData像元会被自动跳过所以得到的点集就是“有效测量值”的离散样本。第二步用插值工具生成完整表面。ArcGIS里最常用的插值工具是IDW反距离权重、克里金Kriging和样条函数Spline。以IDW为例IDW适用于数据变化比较平缓的场景参数较少运算快。常用幂指数power为2搜索半径可以用“变量”模式默认采样点数为12。克里金适用于有明确空间相关性的数据计算更精细但参数多半变异函数模型选择球面、指数、高斯等需要一定的统计学基础。如果数据平滑性较好选“球面”如果数据波动较大可以试“指数”。样条函数能生成非常平滑的表面但可能在边界处出现“过冲”导致产生超出原始数据范围的极值这一点对地形数据尤其要注意。第三步用Con函数把插值结果只填入空洞。表达式Con(IsNull(dem), interpolated_surface, dem)这里有几条经验提醒一下。插值工具默认的输出范围是输入样本点的整体外接矩形也就是说插值结果会覆盖到原始栅格没有数据的外侧区域。使用Con合并时这些外侧区域原本不是NoData但因为原始栅格在那里也是NoDataCon会把这些区域也填上插值结果。如果你的原始栅格范围本来就比研究区大这没问题如果你的NoData区域仅限于内部空洞而研究区外的NoData是“天然空缺”那你得先用Extract by Mask或者SetNull把研究区外的插值结果重新设回NoData。解决方法是把研究区范围作为掩膜让插值只在有效范围内输出ExtractByMask(interpolated_surface, study_area_boundary)再用Con合并。这个方法稍显麻烦但能确保“研究区内的空洞被填好研究区外的天然空缺保持不变”。空间插值法的最大优势是填补结果平滑、连续视觉效果和数学意义上的合理性都更好。缺点是参数设置不当容易出现极端值而且运算量比前两种方法大得多。如果只是填补几个孤立小洞杀鸡用牛刀不仅慢还可能引入不必要的“平滑假象”所以要根据空洞大小选择方法。4. ArcGIS Pro的“填充NoData”栅格函数和Python批处理4.1 栅格函数填空洞也有“黑盒”方案ArcGIS Pro里有一个专门的“填充NoDataFill NoData”栅格函数位置在“分析”选项卡 → “栅格函数” → “填充NoData”。这个函数不需要写表达式对话框里设置几个参数就能直接把NoData填掉。它的原理其实结合了趋势面拟合和焦点统计内部通过多次迭代估算空洞区域的值直到填补完成或达到预设的最大迭代次数。和ArcMap里手动组合工具不同这个栅格函数对大面积空洞的支持要好很多。参数方面主要设置最大迭代次数默认是10和邻域窗口大小默认是3x3。我在实际使用中如果空洞宽度在10像元以内默认参数就能填得很好如果空洞宽度超过20像元我会把最大迭代次数调到20~30窗口大小调到5x5效果会明显提升。这个函数的缺点是比较“黑盒”你很难精确控制每一个空洞的填补逻辑对做科研或者需要高精度成果的场合来说可解释性弱了一些。不过如果只是项目交付前“让数据好看一点”或者做快速预处理这个工具的性价比很高。4.2 Python脚本自动迭代万无一失手动在栅格计算器里嵌套表达式处理几次迭代效率太低而且容易出错。但我写了下面这段Python脚本之后迭代填补这件事就变成了一键操作。这个脚本结合了Con、IsNull和FocalStatistics自动迭代直到没有NoData为止import arcpy from arcpy.sa import * arcpy.CheckOutExtension(Spatial) arcpy.env.workspace rD:\your_workspace arcpy.env.overwriteOutput True in_raster dem.tif out_raster dem_filled.tif # 如果栅格路径中有空格要加上引号 temp in_raster max_iterations 10 # 设置最大迭代次数防止死循环 for i in range(max_iterations): # 统计当前NoData数量 mask IsNull(temp) nodata_count_result arcpy.GetRasterProperties_management(mask, SUM) nodata_count float(nodata_count_result.getOutput(0)) print(f第 {i1} 次迭代前 NoData 像元数: {int(nodata_count)}) if nodata_count 0: print(NoData已全部填充完毕) break # 当前窗口大小随迭代次数增加而增大 # 这样既填小洞也逐步扩大范围去补大洞 window_size 3 i * 2 # 3, 5, 7, 9... filled Con(IsNull(temp), FocalStatistics(temp, NbrRectangle(window_size, window_size, CELL), MEAN, NODATA), temp) temp filled temp.save(out_raster) print(f处理完成结果保存在: {out_raster}) arcpy.CheckInExtension(Spatial)这个脚本做的事情是每轮迭代检查当前栅格里是否还有NoData如果有就根据当前迭代次数自动调整窗口大小继续填充。因为窗口从3x3开始逐步增加到5x5、7x9、9x9……所以就算是大面积空洞随着迭代次数增加也能被慢慢填满。GetRasterProperties_management(mask, SUM)这段用来统计IsNull结果中值为1的像元数也就是NoData像元总数精确又方便。我最开始写这段脚本时犯过一个错误窗口大小固定为3x3结果大空洞填了20轮都还有NoData程序跑了好久也没结束。后来把窗口改成动态递增效率立刻提升了。这个经验也分享出来省得大家走弯路。如果有多幅栅格需要批量处理可以在外层再加一层For循环遍历文件夹里的所有文件把in_raster换成循环变量即可。批处理时要注意输出路径不能和输入路径相同否则会覆盖原数据处理顺序乱了就找不回原始数据了。5. 填补之后的验证和那些绕不开的坑5.1 怎么判断填补效果到底好不好填补完成以后别急着高兴先花几分钟做一下效果验证。我见过不少人在填补后直接拿去做分析结果因为某些位置补出来的值离谱分析结果还是很奇怪走了一大圈弯路才发现是填充环节出了问题。验证主要看几点一是统计指标变化。打开填充前后栅格的属性表对比均值、标准差、最小值和最大值。如果填补之后最小值突然变成-9999或者最大值变成了几百上千的异常值那说明填充过程把某些NoData值错误地带入了有效数据或者插值过程产生了过冲极值。二是视觉检查。把填充区域用单独的符号色渲染出来反复放大缩小查看看是否存在明显的“补丁感”、条带状伪影、或者与周边数据格格不入的突兀值。地形的自然变化通常是有连续性的如果填补区域出现一圈圈的年轮状纹理或者突然出现一个山尖或坑洞就要检查参数设置了。三是边界连续性。对DEM这类连续表面数据重点检查空洞边界附近有没有明显的台阶状跳跃。一个简单方法是在ArcScene或ArcGIS Pro的3D视图里把栅格叠加显示旋转视角从侧面观察看填补区域和周边区域是否形成自然过渡。5.2 常见问题排查我根据自己的踩坑经历把填补Nodata最常见的几个问题整理了一下现象可能原因解决方案焦点统计填充后仍是NoData窗口太小邻域内全是NoData或FocalStatistics的NoData参数没设为“NODATA”增大窗口修改参数并重新计算填补区域出现块状补丁空洞面积偏大Nibble产生“近值复制”效应改用空间插值方法插值结果出现极端值Spline过冲IDW幂指数设置不当搜索半径太小换克里金或降低幂指数填补后范围扩大到研究区外Con合并时没有裁剪插值表面覆盖了外圈先用ExtractByMask或SetNull处理插值结果栅格计算器报错说表达式无效栅格名称带空格或路径为中文没勾选处理扩展模块栅格名称加双引号CheckOutExtension(Spatial)填补结果整体偏移或模糊重投影/重采样时方法选的是双线性或三次卷积值本身被平滑对分类或离散数据用最近邻对连续数据用双线性但接受一定平滑NoData被当成0参与运算栅格属性中NoData值未正确识别用“构建栅格属性表”或复制栅格工具重新定义NoData值这里再展开一个特别容易踩的坑NoData值的定义出了问题。比如你用某些工具导出的TIFF文件NoData可能被写在文件元数据里但ArcGIS有时不一定正确读取导致NoData被当成普通数值参与运算。这时先右键栅格打开属性看“NoData值”一栏是否正确如果不正确最简单的办法是使用“Raster Calculator”里的SetNull(dem 特殊值, dem)把那个特殊值重新定义回NoData再处理。5.3 个人经验和建议最后分享几条我在实际项目里总结出来的经验。第一操作前永远备份一份原始栅格。听起来是废话但我真的见过有人直接覆盖原文件填完发现参数没调对想回头找原始数据已经来不及了。ArcGIS的overwrite默认关闭还好如果开了覆盖模式一不小心保存到原路径那就只能重新下载或者重新处理原始数据了。第二在成果说明里记录填补过程。如果是给科研论文或者工程报告用一定要写清楚用的什么方法、什么参数、什么时候处理的。尤其你是用插值法的时候那部分数据已经带了模型的不确定性。很多评审专家都会盯这一块。第三填补和数据分析尽量在同一坐标系下完成。坐标系统一能减少重采样带来的额外误差。之前遇到过有人先在WGS84下填补再转到CGCS2000结果填补好的位置因为重采样又出现了新的NoData来回折腾了很久。第四不要迷信某一个工具。“一招鲜吃遍天”在GIS数据处理里是不存在的。大洞、小洞、边缘空缺每一种情况都有最适合的方案。我自己在小空洞场景用焦点统计大面积空洞用空间插值边缘条带用趋势面外推加裁剪这算是一个比较成熟的组合策略。第五填补完以后顺手检查一下下游分析的输入。比如你填的是DEM填完之后跑一遍填洼工具和流向工具看看是否还有报错或者大面积NoData输出如果你填的是遥感影像填完以后算NDVI看看是否有异常的饱和区。这种复盘检查能帮你及时发现填补阶段遗留的问题而不是等整个流程走完了再回头找。处理栅格Nodata这件事说难不难说简单也不简单。难的是选对方法、调好参数、控制误差简单的是一旦把各种方法吃透了整个过程其实就是几分钟的操作。希望在看完这篇文章之后你遇到Nodata不用再头皮发麻从容打开ArcGIS按部就班把它处理得干干净净。
返回列表