GDAL `gdal raster sieve` 详解：用大小阈值合并清除栅格小面（4/8 连通、掩膜与流水线用法）

GIS遥感数据工程【免费下载链接】gdalGDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.项目地址https://gitcode.com/gh_mirrors/gd/gdal点击查看免费下载本文是 GDAL 官方命令文档 gdal_raster_sieve.rst 的深度实战指南。gdal raster sieve是 GDAL 3.11 起加入新版gdal命令家族的栅格子命令用于把面积小于指定阈值像素数的栅格“小面”polygon合并到其最大相邻面中是栅格分类结果去噪、地类图斑破碎化消除的常用工具。读完本文你将掌握该子命令的全部参数语义、与底层GDALSieveFilter()算法的对应关系、掩膜与连通性的实际影响以及在gdal raster pipeline中把它作为流水线步骤的方法。一、命令总览与适用场景gdal raster sieve的作用一句话概括移除小于给定阈值以像素计的栅格多边形并用其最大相邻多边形的像素值替换它们。这里的“多边形”并非矢量多边形而是指栅格中由像素值相同且彼此连通connected的区域。典型应用场景包括分类影像如土地覆盖分类结果中的“椒盐”噪声斑块清理将破碎的小图斑合并到周围主导地类中简化分析对象在gdal raster pipeline中作为一步自动化处理串联在其他栅格步骤之后。该子命令自GDAL 3.11引入versionadded:: 3.11从GDAL 3.12起它还可以作为 gdal raster pipeline 的一个潜在步骤使用。1.1 与旧版 Python 工具 gdal_sieve 的关系在引入本子命令之前GDAL 一直提供独立的 Python 脚本gdal_sieve见旧版命令文档 gdal_sieve.rst其用法为gdal_sieve [--help] [--help-general] [-q] [-st threshold] [-4] [-8] [-o namevalue] srcfile [-nomask] [-mask filename] [-of format] [dstfile]该脚本依赖 GDAL Python 绑定才能运行。官方文档明确指出gdal raster sieve是它在新版gdal命令行界面中的等价实现两者共享同一个底层算法。二、Synopsis命令行语法运行gdal raster sieve --help-doc可获得完整语法原文档通过program-output指令在构建时自动嵌入。结合源码 gdalalg_raster_sieve.h 中定义的三类参数完整的调用形态为gdal raster sieve [--help] [--help-doc] [--version] [-b band] [-s size-threshold] [-c] [--mask mask] [--append] [--co namevalue] [--if format] [--oo namevalue] [-f format] [--overwrite] input outputinput输入栅格数据集output输出栅格数据集GDALG 非流式输出。新版gdal命令遵循统一帮助约定--help展示概览--help-doc展示完整文档--version输出版本信息。三、程序专属选项详解文档定义了 4 个程序专属选项它们在实现层面对应 gdalalg_raster_sieve.cpp 中注册的参数3.1 -b, --band 选择输入波段-b band指定参与筛分的输入波段1 起始索引。默认值为1源码中int m_band 1。当多波段数据只希望清理某一个波段时使用。注意输出数据集只包含被处理的这一个波段。3.2 -s, --size-threshold 大小阈值-s size-threshold保留多边形的最小面积像素数默认值为 2源码中int m_sizeThreshold 2。小于该阈值的多边形将并入其最大相邻多边形。3.3 -c, --connect-diagonal-pixels对角像素连通-c控制多边形连通性的判定方式默认不加-c只把边接触touching the edges的像素视为连通等价于4-连通4-connectivity加-c后额外把角接触corners的像素也视为连通等价于8-连通8-connectivity。在实现层该布尔值被直接映射为连通度参数传给底层算法源码 gdalalg_raster_sieve.cpp 中m_connectDiagonalPixels ? 8 : 4一目了然——false传 4true传 8。选择 8-连通通常会使多边形更容易“长大”小图斑更易被合并。3.4 --mask有效性掩膜--mask mask使用指定文件的第一波段作为有效性掩膜validity mask掩膜波段中值不为零的像素才被认为“适合纳入多边形”即参与筛分值为零的像素视为无效/NoData 区域既不会成为多边形的一部分也不会被修改、不会影响多边形面积。源码实现中掩膜波段取自m_maskDataset.GetDatasetRef()-GetRasterBand(1)即掩膜文件的第一波段如获取失败会报错Cannot get mask band.。这在处理带真实背景/无效区域的栅格时非常有用——例如掩膜外的海面、云区被排除在统计之外从而避免背景像素形成“巨型多边形”把边缘小面误合并掉。四、标准选项说明文档中以 collapse 形式内联引用了gdal_options/目录下的通用选项本节逐条展开对应文件均位于 gdal_options 目录选项语法说明--append--append将输入栅格作为新子数据集subdataset追加到已有输出文件中仅对支持追加子数据集的驱动有效如 GeoTIFF、GPKG输出文件不存在时会先创建append_raster.rst。--co--co NAMEVALUE输出创建选项可多次重复。不同格式驱动的创建选项各不相同例如 GeoTIFF 支持COMPRESS、TILED等可用gdal --formats | grep raster | grep rw | sort查看候选驱动co.rst。--if--if format指定用于打开输入文件的格式/驱动名可重复指定多个候选驱动用于跳过自动驱动探测if.rst。--oo--oo NAMEVALUE数据集打开选项格式相关可重复oo.rst。-f/--of/--format/--output-format-f OUTPUT-FORMAT指定输出栅格格式of_raster_create_copy.rst。--overwrite--overwrite允许覆盖已存在的目标文件/数据集默认情况下若目标已存在命令会直接报错退出overwrite.rst。其中--overwrite行为在测试 test_gdalalg_raster_sieve.py 中有明确验证不传--overwrite再次运行会抛出already exists异常传--overwrite后即可成功重跑。五、返回值Return status code执行成功返回状态码0出错时返回非零状态码。需要注意以警告形式发出的非阻塞错误non-blocking errors仍被视为成功执行。该约定在 return_code.rst 中统一定义适用于所有gdal子命令。六、数据读取语义整数化与浮点问题文档明确指出输入数据集按整数数据读取浮点值会被四舍五入为整数。这意味着对浮点栅格执行筛分时原始小数信息在算法内部被取整某些场景可能需要预先重缩放re-scaling源数据文档给出的典型例子是32 位浮点数据、取值范围 min0 ~ max1——此时绝大多数像素在取整后都变成 0/1 两个值栅格细节几乎全部丢失必须先用gdal raster scale或gdal raster convert之类步骤把数据放大到有意义的整数量级。这一行为在底层算法中体现得更为彻底gdalsievefilter.cpp 在读取输入时即以GDT_Int64类型执行GDALRasterIO即无论原始波段是什么类型都会按 64 位整数语义读入参与多边形编号与合并。这也是为什么文档特别提醒浮点数据要谨慎、必要时先重缩放。七、底层算法原理GDALSieveFilter 的三遍扫描gdal raster sieve最终调用 C 层算法接口GDALSieveFilter()其声明位于 gdal_alg.h实现位于 gdalsievefilter.cpp。该接口同样可以脱离命令行被 C/C 程序直接调用CPLErr GDALSieveFilter(GDALRasterBandH hSrcBand, GDALRasterBandH hMaskBand, GDALRasterBandH hDstBand, int nSizeThreshold, int nConnectedness, char **papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg);其中nConnectedness只能取 4 或 8对应命令行的-c开关papszOptions当前为空文档注释“None currently supported”。从源码可以梳理出算法执行的完整流程第一遍扫描枚举多边形并统计面积逐行读取用GDALRasterPolygonEnumerator按连通性4/8给每个连通区域分配多边形 ID同时累计每个多边形的像素数gdalsievefilter.cpp合并 ID 映射收尾通过CompleteMerges()把第一遍扫描中产生中间编号的多边形片段合并到最终 ID并合并对应面积gdalsievefilter.cpp第二遍扫描寻找最大邻居重新枚举多边形对每个像素与上下左右8 连通时还包括对角的邻居比较为每个小多边形记录“最大相邻多边形”anBigNeighbourgdalsievefilter.cpp链式追认若某小多边形的最大邻居也小于阈值则沿“最大邻居链”继续向上查找直到找到一个不小于阈值的多边形作为最终归宿若找不到如被 NoData 包围的孤立小面则保持原值不动gdalsievefilter.cpp。这一行为与文档说明一致“小于阈值但没有任何达到阈值大小的邻居的多边形不会被改变被 NoData 包围的多边形因此不会被改动”第三遍扫描回写结果再次逐行处理把需要合并的像素值改写为其最终目标多边形的像素值并写出gdalsievefilter.cpp。7.1 时间复杂度与内存占用特征该算法对输入文件做3 遍完整扫描。内存占用与多边形数量成正比源码注释估算约每多边形 24 字节而与栅格总大小无直接关系。因此大片平滑的栅格即使文件很大也能高效处理极“嘈杂”、充满大量单像素多边形的栅格会因多边形数量激增而消耗大量内存。测试 sieve.py 中的性能用例验证了“处理时间随像素数大致线性增长”的特性。7.2 掩膜与 NoData 的边界行为被掩膜排除掩膜值为 0的像素无论原始值是什么都不属于任何多边形算法把 NoData 标记的多边形GP_NODATA_MARKER直接忽略不参与合并gdalsievefilter.cpp若所有像素都被掩膜全掩膜情形输入输出相同时直接成功返回不同时则退化为整波段拷贝gdalsievefilter.cpp测试 sieve.py 覆盖了这一边界。八、非原生流式输出GDALGGDAL 3.12 起本命令支持通过GDALG输出格式把整条命令行序列化为 JSON 文件随后可用 gdalg 驱动以栅格数据集方式打开该 JSON实现按需on-the-fly/流式执行同一处理管线见 gdalg_raster_compatible_non_natively_streamable.rst。但需要特别注意sieve 算法并非原生流式兼容源码中它继承自 GDALRasterPipelineNonNativelyStreamingAlgorithm并覆写了IsNativelyStreamingCompatible()返回不兼容。因此在实际执行时会先生成一个临时数据集——从源码看RunStep() 首先调用CreateTemporaryCopy()把输入复制为临时栅格再在临时栅格上原地执行GDALSieveFilter()。这意味着打开 GDALG 输出时可能产生显著的处理时间因为必须先把数据物化到临时文件临时文件的创建可通过配置项控制测试 test_gdalalg_raster_sieve.py 展示了GDAL_RASTER_PIPELINE_USE_GTIFF_FOR_TEMP_DATASET与CPL_TMPDIR对临时文件行为的约束。九、在 gdal raster pipeline 中作为步骤自 GDAL 3.12 起sieve 可作为 gdal raster pipeline 的一步gdal_raster_pipeline.rst 中可见其被列为 pipeline 支持的步骤之一gdal raster pipeline --help-docsieve可查看该步骤的专属帮助参数与独立使用时完全一致。例如把“裁剪 → 筛分 → 写回”串成一条管线gdal raster pipeline \ gdal raster read --input /path/to/input.tif \ gdal raster clip ... \ gdal raster sieve --size-threshold 10 --connect-diagonal-pixels \ gdal raster write --output /path/to/output.tif注意 pipeline 中需要显式指定read与write步骤来界定数据流的起止sieve 作为非原生流式步骤会触发临时数据集物化。十、示例与验证10.1 文档示例清理波段 2 中小于 10 像素的多边形$ gdal raster sieve -b 2 -s 10 input.tif output.tif该命令读取input.tif的第 2 波段把面积小于 10 像素的多边形并入其最大相邻多边形结果写入output.tif。此例完整对应文档示例gdal_raster_sieve.rst。10.2 常见变体# 默认波段、默认阈值 2输出为 GTiff $ gdal raster sieve input.grd output.tif # 8 连通 阈值 5 掩膜文件 $ gdal raster sieve -c -s 5 --mask mask.tif input.tif output.tif # 输出 LZW 压缩的 GeoTIFF并允许覆盖已存在文件 $ gdal raster sieve -s 10 --co COMPRESSLZW --overwrite input.tif output.tif10.3 用测试数据亲手验证仓库自带筛分测试数据 sieve_src.grdAAIGrid 文本格式5×7、NODATA132、含 107/115/123/140/148/156/100/101/102/103 等多值混合可直接复现官方算法测试 sieve.py 的预期结果gdal raster sieve -s 2 input.grd out_4conn.tif # 4 连通校验和 364 gdal raster sieve -c -s 2 input.grd out_8conn.tif # 8 连通校验和 370自动化测试 test_gdalalg_raster_sieve.py 以参数化方式connect_diagonal_pixels取 False/True、创建选项取{}/TILEDYES/COMPRESSLZW逐一验证了 4/8 连通的输出校验和364 vs 370以及 TILED/压缩创建选项是否生效——这也印证了-c开关与--co选项在真实输出中可观测的差异。十一、实践建议与注意事项小结先确认数据是整数类型浮点栅格尤其 0~1 范围会因取整而失真建议先用gdal raster scale重缩放合理选择连通性-c8 连通合并更激进适合碎斑密集的影像默认 4 连通更保守善用掩膜隔离无效区用--mask排除背景避免无效值像素参与多边形统计也避免孤立小面“无处可并”控制多边形总数极嘈杂数据会因多边形数量巨大而内存占用飙升必要时先做中值/众数滤波类预处理输出覆盖需显式声明目标已存在时默认报错需要--overwrite在 pipeline 中注意物化开销sieve 非流式兼容作为 pipeline 步骤时会生成临时数据集。十二、参考文档与源码索引命令文档gdal_raster_sieve.rst算法 C 接口GDALSieveFilter 声明、实现 gdalsievefilter.cpp子命令实现gdalalg_raster_sieve.cpp、gdalalg_raster_sieve.h自动化测试test_gdalalg_raster_sieve.py、sieve.py、test_gdal_sieve.py旧版 Python 工具文档gdal_sieve.rst通用选项定义gdal_options返回值约定return_code.rst流水线文档gdal_raster_pipeline.rst赞分享GIS遥感数据工程【免费下载链接】gdalGDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.项目地址https://gitcode.com/gh_mirrors/gd/gdal点击查看免费下载相关推荐GDAL gdal raster select 命令详解栅格波段选择、重排与掩膜处理实战指南GDAL gdal raster select 命令详解栅格波段选择、重排与掩膜处理实战指南 gdal raster select 是 GDAL 3.11 引GIS遥感数据工程GDAL 栅格概览删除命令详解gdal raster overview delete 的用法与底层实现GDAL 栅格概览删除命令详解 gdal raster overview delete 的用法与底层实现 导读 gdal raster overview deGIS遥感数据工程GDAL gdal raster contour 等值线提取完全指南从 DEM 栅格生成矢量等值线与等值面GDAL gdal raster contour 等值线提取完全指南从 DEM 栅格生成矢量等值线与等值面 本篇技术指南围绕 GDAL 3.11 引入的新一代GIS遥感数据工程上一篇CKEditor 5 Widget 内部机制深度剖析type-around 功能的禁用之道与 data-cke-ignore-events 事件隔离下一篇es-toolkit asyncNoop 详解异步空操作函数的实现原理与工程实践创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考