企业官网建设流程全解析

• GIS
• 遥感
• 数据工程

本文是 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 name=value] <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 <name>=<value>] [--if <format>] [--oo <name>=<value>] [-f <format>] [--overwrite] <input> <output>

• <input>：输入栅格数据集；
• <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传 4，true传 8。选择 8-连通通常会使多边形更容易“长大”，小图斑更易被合并。

3.4 --mask：有效性掩膜

--mask <mask>

使用指定文件的第一波段作为有效性掩膜（validity mask）：掩膜波段中值不为零的像素才被认为“适合纳入多边形”，即参与筛分；值为零的像素视为无效/NoData 区域，既不会成为多边形的一部分，也不会被修改、不会影响多边形面积。

源码实现中掩膜波段取自m_maskDataset.GetDatasetRef()->GetRasterBand(1)（即掩膜文件的第一波段），如获取失败会报错Cannot get mask band.。这在处理带真实背景/无效区域的栅格时非常有用——例如掩膜外的海面、云区被排除在统计之外，从而避免背景像素形成“巨型多边形”把边缘小面误合并掉。

四、标准选项说明

文档中以 collapse 形式内联引用了gdal_options/目录下的通用选项，本节逐条展开（对应文件均位于 gdal_options 目录）：

其中--overwrite行为在测试 test_gdalalg_raster_sieve.py 中有明确验证：不传--overwrite再次运行会抛出already exists异常，传--overwrite后即可成功重跑。

五、返回值（Return status code）

执行成功返回状态码0；出错时返回非零状态码。需要注意：以警告形式发出的非阻塞错误（non-blocking errors）仍被视为成功执行。该约定在 return_code.rst 中统一定义，适用于所有gdal子命令。

六、数据读取语义：整数化与浮点问题

文档明确指出：输入数据集按整数数据读取，浮点值会被四舍五入为整数。这意味着：

• 对浮点栅格执行筛分时，原始小数信息在算法内部被取整；
• 某些场景可能需要预先重缩放（re-scaling）源数据，文档给出的典型例子是32 位浮点数据、取值范围 min=0 ~ max=1——此时绝大多数像素在取整后都变成 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”）。

从源码可以梳理出算法执行的完整流程：

1. 第一遍扫描（枚举多边形并统计面积）：逐行读取，用GDALRasterPolygonEnumerator按连通性（4/8）给每个连通区域分配多边形 ID，同时累计每个多边形的像素数（gdalsievefilter.cpp）；
2. 合并 ID 映射收尾：通过CompleteMerges()把第一遍扫描中产生中间编号的多边形片段合并到最终 ID，并合并对应面积（gdalsievefilter.cpp）；
3. 第二遍扫描（寻找最大邻居）：重新枚举多边形，对每个像素与上下左右（8 连通时还包括对角）的邻居比较，为每个小多边形记录“最大相邻多边形”anBigNeighbour（gdalsievefilter.cpp）；
4. 链式追认：若某小多边形的最大邻居也小于阈值，则沿“最大邻居链”继续向上查找，直到找到一个不小于阈值的多边形作为最终归宿；若找不到（如被 NoData 包围的孤立小面），则保持原值不动（gdalsievefilter.cpp）。这一行为与文档说明一致：“小于阈值但没有任何达到阈值大小的邻居的多边形不会被改变，被 NoData 包围的多边形因此不会被改动”；
5. 第三遍扫描（回写结果）：再次逐行处理，把需要合并的像素值改写为其最终目标多边形的像素值并写出（gdalsievefilter.cpp）。

7.1 时间复杂度与内存占用特征

该算法对输入文件做3 遍完整扫描。内存占用与多边形数量成正比（源码注释估算约每多边形 24 字节），而与栅格总大小无直接关系。因此：

• 大片平滑的栅格即使文件很大也能高效处理；
• 极“嘈杂”、充满大量单像素多边形的栅格会因多边形数量激增而消耗大量内存。

测试 sieve.py 中的性能用例验证了“处理时间随像素数大致线性增长”的特性。

7.2 掩膜与 NoData 的边界行为

• 被掩膜排除（掩膜值为 0）的像素无论原始值是什么都不属于任何多边形；
• 算法把 NoData 标记的多边形（GP_NODATA_MARKER）直接忽略，不参与合并（gdalsievefilter.cpp）；
• 若所有像素都被掩膜（全掩膜情形），输入输出相同时直接成功返回，不同时则退化为整波段拷贝（gdalsievefilter.cpp），测试 sieve.py 覆盖了这一边界。

八、非原生流式输出（GDALG）

GDAL 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-doc=sieve可查看该步骤的专属帮助，参数与独立使用时完全一致。例如，把“裁剪 → 筛分 → 写回”串成一条管线：

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 COMPRESS=LZW --overwrite input.tif output.tif

10.3 用测试数据亲手验证

仓库自带筛分测试数据 sieve_src.grd（AAIGrid 文本格式，5×7、NODATA=132、含 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、创建选项取{}/TILED=YES/COMPRESS=LZW）逐一验证了 4/8 连通的输出校验和（364 vs 370）以及 TILED/压缩创建选项是否生效——这也印证了-c开关与--co选项在真实输出中可观测的差异。

十一、实践建议与注意事项小结

1. 先确认数据是整数类型：浮点栅格（尤其 0~1 范围）会因取整而失真，建议先用gdal raster scale重缩放；
2. 合理选择连通性：-c（8 连通）合并更激进，适合碎斑密集的影像；默认 4 连通更保守；
3. 善用掩膜隔离无效区：用--mask排除背景，避免无效值像素参与多边形统计，也避免孤立小面“无处可并”；
4. 控制多边形总数：极嘈杂数据会因多边形数量巨大而内存占用飙升，必要时先做中值/众数滤波类预处理；
5. 输出覆盖需显式声明：目标已存在时默认报错，需要--overwrite；
6. 在 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
• 遥感
• 数据工程