使用arcpy移除遥感影像云层-编程知识

先讲思路，然后上代码：

去除云层

思路1：

如果同一地理区域的多个图像，其中一些部分有丰富的云，而另一些部分没有云，则可以将它们组合起来，以便无云的部分替代多云的部分。这种方法很简单，但没有考虑季节性变化或随时间变化的差异；因此，相同的地方通常看起来可能会有所不同。在这种情况下，可以逐一合并来自相近日期的高频场景以生成合成图像并减少图像上的季节性差异。

思路2：

考虑阴影或云层覆盖；这种技术通常称为像素排序。首先，根据颜色和云层覆盖百分比预先选择多个图像。然后，可以使用机器学习算法重新排序和合并所有图像，或者对每个像素使用一些马赛克运算符（例如均值、中值和混合等）。

思路3：

此外，假设可以使用高性能计算处理器。在这种情况下，可以遵循自动化方法或机器学习分类，例如随机森林或支持向量机。可以拍摄同一场景的多张图像，并根据一段时间内像素的模式观察（即比较像素及其周围环境的相似性）来合并它们，从而给出一个没有云的位置的外观。

数据下载与去除云层代码

在这篇文章中，我将分享一个基于 ArcGIS Pro 的 ArcPy 站点包的代码示例。

使用 Landsat 9

这个基本示例使用来自陆地卫星的图像，这些卫星是被动遥感系统（即，它们仅测量由太阳传输并由目标物体反射或发射的辐射）。Landsat地面数据接收、数据处理、数据归档、产品生成和产品分发请参阅官网。Landsat 9 是 Landsat 系列中最新的卫星。它于 2021 年 9 月发射，很大程度上复制了其前身 Landsat 8 的属性。它携带两台具有中等空间分辨率传感器的仪器（15米、30米和100米，取决于光谱带）。

https://earthexplorer.usgs.gov/

主要步骤：

- 在哪里下载 Landsat 图像？我使用USGS EarthExplorer 工具（可以免费注册）。它是一个广泛用于下载多个卫星图像或地理空间数据的便捷平台。要消除云层覆盖，通常需要同一地理位置的一组重要图像（例如，共享轨道）。
- 因为国外的数据全一些，经过对比选择，我确定了哥伦比亚 (74°28'44''W 4°21'51''N) 不同日期的先验六幅图像（来自 Landsat 8/9 Collection 2 Level 2），云量低于 40%，轨道路径=008，行=057。
- 我使用像素质量评估带（QA_PIXEL 栅格文件，包含在选定的 Landsat 数据产品中），其中包含从图像数据和云掩模信息收集的质量统计数据。用于保持清澈区域和水域的高置信度像素值分别为 21824 和 21952。目前，ArcGIS 只允许使用具有相同波段数的栅格分配条件空值；因此我使用“复合波段”地理处理工具。
- 我使用“镶嵌到新栅格”地理处理工具将文件合并到新栅格中。在这种情况下，重叠图像的输出像元值将是重叠像元的平均值，而其剩余点（仍然有云）现在将具有“NoData”值。

在本例中，我使用软件ArcGIS Pro 3.0.2。该脚本是我开发的一个较大脚本的简化版本，用于处理大量文件夹和文件；因此，该代码的某些步骤可以轻松调整。强烈建议创建和使用地理数据库文件。在此示例中，该文件用于保存输出并包含在“Pathname2”的末尾。下面是不同时期的5副影像及其云量！

哥伦比亚 (74°28'44''W 4°21'51''N)，2020 年 8 月 29 日 — 云量：27.6%

哥伦比亚 (74°28'44''W 4°21'51''N)，2020 年 9 月 14 日 — 云量：32.0%

哥伦比亚 (74°28'44''W 4°21'51''N)，2021 年 1 月 4 日 — 云量：34.6%

哥伦比亚 (74°28'44''W 4°21'51''N)，2021 年 9 月 17 日 — 云量：31.0%

哥伦比亚 (74°28'44''W 4°21'51''N)，2021 年 10 月 3 日 — 云量：28.7%

哥伦比亚 (74°28'44''W 4°21'51''N)，2021 年 12 月 14 日 — 云量：28.8%

# Import arcpy module          
import arcpy          # Overwrite          
arcpy.env.overwriteOutput = True  # Overwriting saved files          # Local variables          
landsat_folders = "Bulk Order 008057" # Bulk folders downloaded from USGS Earth Explorer          
frame = ["LC08_L2SP_008057_20200829_20200906_02_T1", "LC08_L2SP_008057_20200914_20200919_02_T1",          "LC08_L2SP_008057_20210104_20210309_02_T1", "LC08_L2SP_008057_20210917_20210925_02_T1",          "LC08_L2SP_008057_20211003_20211013_02_T1", "LC09_L2SP_008057_20211214_20220120_02_T1"]          # Workspaces: Pathname1=where you will include the Landsat folders, Pathname2=where you would extract the output          
workspaces = ["Pathname1", "Pathname2"]          
arcpy.env.workspace = workspaces[0]          for j in range(len(frame)):          qapixel = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_QA_PIXEL.tif"          band1 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B1.tif"          band2 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B2.tif"          band3 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B3.tif"          band4 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B4.tif"          band5 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B5.tif"          band6 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B6.tif"          band7 = landsat_folders + "/Landsat 8-9 OLI_TIRS C2 L2/" + frame[j] + "/" + frame[j] + "_SR_B7.tif"                arcpy.management.CompositeBands([qapixel, qapixel, qapixel, qapixel, qapixel, qapixel, qapixel], workspaces[1] + "comp_qapixel")          arcpy.management.CompositeBands([band1, band2, band3, band4, band5, band6, band7], workspaces[1] + "comp_bands")          
outSetNull = arcpy.sa.SetNull(workspaces[1] + "comp_qapixel", workspaces[1] + "comp_bands", "VALUE <> 21824 AND VALUE <> 21952")          
outSetNull.save("Pathname2" + "c_" + frame[j][17:25])          # Mosaic to New Raster. This step allows you to create the final image with no/less clouds and delete unnecessary outputs          
arcpy.management.MosaicToNewRaster([workspaces[1] + "c_" + frame[0][17:25], workspaces[1] + "c_" + frame[1][17:25], workspaces[1] + "c_" + frame[2][17:25], workspaces[1] + "c_" + frame[3][17:25], workspaces[1] + "c_" + frame[4][17:25]], workspaces[1] + "c_" + frame[5][17:25]], workspaces[1], "new" + frame[j][10:16], pixel_type = "16_BIT_UNSIGNED", number_of_bands = 7, mosaic_method = "MEAN")          
arcpy.management.Delete([workspaces[1] + "c_" + frame[0][17:25], workspaces[1] + "c_" + frame[1][17:25], workspaces[1] + "c_" + frame[2][17:25], workspaces[1] + "c_" + frame[3][17:25], workspaces[1] + "c_" + frame[4][17:25], workspaces[1] + "c_" + frame[5][17:25], workspaces[1] + "comp_qapixel", workspaces[1] + "comp_bands"])