从卫星图像创建无云图像合成

通常我们在处理光学图像时经常遇到的问题之一是云和阴影的存在， 因为它们会阻挡图像中的地面信息。 对于小区域， 可以等待下一张清晰的图像。 然而，对于被多个足迹覆盖的大面积区域， 尤其是海岛或热带等阴天频繁的地区， 几乎不可能获得清晰的图像。 在实践中， 我们使用在一个时间窗口内收集的多个场景来创建无云图像合成。

ArcGIS Pro 3.1 中的新几何中线是用于图像合成工作流的栅格函数。 它通过最小化到所有重叠像素的距离之和来计算光谱空间中的最佳像素近似值， 有助于减少异常像素, 并在保持波段之间光谱关系的同时产生更好的图像合成。 本文将向大家展示如何应用此新功能从堆叠的 Landsat 图像创建无云图像合成。

研究区域和数据

研究区域是澳大利亚最南端的塔斯马尼亚岛， 大部分岛屿都被森林覆盖， 北部海岸和贯穿塔斯马尼亚东部的主要河谷系统中有一些混合农业用途。 我们的目标是使用 Landsat 场景绘制 2021 年的岛屿地图。 该岛被 12 个 Landsat 足迹覆盖（左下方）， 从 Landsat Collection 2 中的 Landsat 8 和 Landsat 9 中总共选择了 121 张可用图像， 使用云量标准小于 60%， 以确保每个位置都被清晰的像素覆盖。 右图显示了一些示例。

在 ArcPy 工作流中使用新的几何中值栅格函数从这些图像中创建一个合成图。 最后还提供了用于创建无云图像立方体的代码片段。 对于那些喜欢手动工作流的人， 还将介绍该工作流中使用的相应工具和功能链。

第 1 步：创建图像集

第一步是从图像堆栈创建栅格集合。 尽管可以使用 ArcPy 从源文件直接创建栅格集合对象， 但很多人喜欢先使用地理处理工具创建镶嵌数据集， 即持久栅格集合， 因为这对 ArcPy 工作流和手动工作流都有好处。 具体来说，为添加栅格到镶嵌数据集工具使用多波段处理模板。 此多波段模板会将表面反射波段和 QA 波段组合到一个栅格中， 允许定义和应用云遮罩。 现在将从导入模块和定义 Image Analyst 许可证开始：

接下来，将创建一个镶嵌数据集并将源图像中的数据添加到一个文件夹中。 它将添加 Landsat 8 和 Landsat 9 图像的表面反射率和 QA 波段， 然后从镶嵌数据集创建一个栅格集合对象， 现在包含 121 个栅格的栅格集合。

有关创建镶嵌数据集的手动工作流，请参阅此链接（https://pro.arcgis.com/en/pro-app/latest/help/data/imagery/creating-mosaic-datasets-wf.htm）。

第 2 步：遮盖云层和阴影

为了创建更好的图像合成， 将使用 Landsat 产品中提供的 QA 波段移除云层和阴影。 QA 类别使用位索引定义， 其中 0 表示填充，1 表示膨胀云，2 表示卷云，3 表示云，4 表示云影。 这些可以在 Landsat Collection 2 (详情参见 https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/media/files/LSDS-1619_Landsat-8-9-C2-L2-ScienceProductGuide-v4.pdf)产品文档（表 6-2，第 13 页）中找到， 可使用转置位栅格函数解释 QA 位以创建 QA 掩码。 下面的代码定义了一个 RemoveCloud 函数， 它创建一个 QA 掩码并将其应用于数据波段。 该函数将用于使用栅格集合对象的 Map 方法处理栅格集合中的每个项目。

对于手动工作流， 下面是用于移除云层和阴影的相应栅格函数链， 可使用栅格函数批处理编辑器将其应用于镶嵌数据集， 这将屏蔽镶嵌数据集中每个栅格的云和阴影。

第 3 步：生成无云图像合成

移除云之后，集合中的栅格将有很多 NoData 漏洞。 下一步是从这些栅格创建一个合成。 新的几何中 Median 栅格函数是一种数值方法， 可根据定义的条件迭代运行以近似获得最佳结果。 参数 epsilon 定义精度， 例如，当连续两次运行的结果接近或超过最大迭代次数时， 迭代将停止，将使用默认值， 然后将其保存为云栅格格式。

现在拥有从 121 张图像生成的图像合成， 无缝且无云。

对于手动工作流， 可使用几何中值函数定义函数模板， 并将其应用于从栅格函数生成栅格地理处理工具以创建无云图像合成。

创建无云图像立方体

如果需要将数据处理成每月一个合成或每年一个合成， 则可以使用组字段， 处理每个组内的图像， 并创建无云图像立方体。 在步骤 2 中完成移除栅格集合的云和阴影后， 请完成以下步骤：

1.添加一个名为 year 的新字段， 使用从 "StdTime" 字段中提取的年份值填充此字段，并按此新字段分组：

2.使用 keys 方法从组字段中获取唯一值列表， 并将几何中值函数应用于每个组。 在此处理过程中， 准备将在下一步中使用的栅格、名称和年份列表。

3.从三个列表中创建栅格集合， 将其转换为多维栅格， 并保存为多维 CRF， 这将是分析就绪图像立方体。

ArcGIS 还提供了其他组合像素的方法， 例如计算平均值、计算重叠像素的中值， 或者使用 NDVI 来寻找绿色像素。 这些计算密集度较低的方法可以是满足各种要求的替代方法， 但几何中值可以保持波段之间的光谱关系， 这对分析工作流程是有利的。