1.Fmask软件对1C级产品进行处理，识别像素类别

不知道Fmask是什么可以先去百度一下
软件下载,链接到github地址
在这里插入图片描述
我下载的是4.5版本，无脑安装即可。
双击打开软件（需要等一会），长这样

路径选择E:\S2\S2A_MSIL1C_20220806T032531_N0400_R018_T49RBQ_20220806T054540\S2A_MSIL1C_20220806T032531_N0400_R018_T49RBQ_20220806T054540.SAFE\GRANULE\L1C_T49RBQ_A037195_20220806T033824
在这里插入图片描述
点击run

状态会发生变化

当然也可以采用python实现，详见这个链接
Python会运行快一些，但是生成的是.img格式，转成tif好像多了些奇怪的值。故还是用软件吧！哪个简单用哪个

显示finished表示运行完毕
在这里插入图片描述
可以看到里面多了一个t文件夹，里面是Fmask.tif文件。

DN值所代表的含义如下：
0 => clear land pixel

1 => clear water pixel

2 => cloud shadow

3 => snow

4 => cloud

255 => no observation

运行完之后，我们会得到一个tif文件，因为我们只需要云掩膜也就是保留0和1的像元（纯净的陆地像素与纯净的水像素），其他和云相关的像元删掉，也就是将云所在的像素去掉。这里采用python实现，代码如下：

from osgeo import gdal, ogr, osr
import os
import datetime

path = r"cloud.tif"

if __name__ == '__main__':
    start_time = datetime.datetime.now()

    inraster = gdal.Open(path)
    inband = inraster.GetRasterBand(1)
    prj = osr.SpatialReference()
    prj.ImportFromWkt(inraster.GetProjection())

    outshp ="mask_cloud.shp"#记得改一下文件名路径
    drv = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outshp):
        drv.DeleteDataSource(outshp)
    polygon = drv.CreateDataSource(outshp)
    poly_layer = polygon.CreateLayer(path[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)
    newfield = ogr.FieldDefn('value', ogr.OFTReal)
    poly_layer.CreateField(newfield)

    gdal.Polygonize(inband, None, poly_layer, 0)
    polygon.SyncToDisk()
    polygon = None

    # 打开生成的 Shapefile 文件
    shp_datasource = ogr.Open(outshp, 1)
    shp_layer = shp_datasource.GetLayer()

    # 遍历要素并删除 value 不为 0 或 1 的要素
    shp_layer_def = shp_layer.GetLayerDefn()
    delete_indices = []
    for i in range(shp_layer.GetFeatureCount()):
        feature = shp_layer.GetFeature(i)
        value = feature.GetField("value")
        if value != 0 and value != 1:
            delete_indices.append(i)
        feature = None

    # 从后往前删除要素，避免索引错位
    for idx in reversed(delete_indices):
        shp_layer.DeleteFeature(idx)

    # 重新同步文件
    shp_layer.ResetReading()
    shp_datasource.SyncToDisk()

    end_time = datetime.datetime.now()
    print("Succeeded at", end_time)
    print("Elapsed time:", end_time - start_time)

注意需要有gdal库，推荐本地安装，这里不会可以看这篇文章。

运行完这段代码后会得到一个掩膜文件，用ArcGIS打开长这样，1值代表有效像素
在这里插入图片描述

2.1C级产品生成云掩膜后再批量进行sen2cor大气校正

这我相信你会的！
然后把校正后的2A级产品用SNAP打开

3.批量重采样与掩膜操作

1.导入2A级产品

你应该知道吧，先重采样到你先要的分辨率，这里我们选择10m.

2.导入我们在第1步生成的云掩膜文件

这时候就要进行我们的掩膜操作了（有什么疑问的去看官方文档，点一下“帮助”认真读一下）
在这里插入图片描述
找到vector Data文件夹，选中之后点击上方菜单栏Vector-import-shp

记得这里选否不然每个多边形都生成一个文件

可以看导在Vector Data下面和Masks下面多了一个mask
双击打开可以看到它的属性值，1代表有效像素，0代表无效像素，我们需要抠掉的像素
在这里插入图片描述

这里最好先重采样到10m，我找了一景重采样之后的影像来做演示
看看效果，mask是否完全覆盖了云

然后进行掩膜操作，说白了就是把有云的地方抠掉，然后再进行后续的镶嵌操作补上来，云多真没办法，哎！误差大。

3.SNAP掩膜操作

在这里插入图片描述
输入输出参数，选好路径即可

处理参数选择使用矢量掩膜，拉到最后选择我们之前用Python生成的那个shp文件

波段那里全选。点击run运行即可。
可以看到掩膜之后相当于把云抠掉了

然后再选取同一传感器，同一地点，相近时间的影像在SNAP软件中做同样的处理
然后把他们镶嵌在一起
在这里插入图片描述
注意是选择dim文件，对话框里面有三个选项卡，I/O Parameters(输入输出参数)，Map Projection Definition
(定义投影)和Variables&Coniditions(变量和条件)。首先我们在输入输出参数选项卡中指定输入的影像，点击下面箭头指向的"加号”即可选择所需影像文件。在本选项卡内还需要对Nme(文件名)、文件输出类型和目录进行设置。
在这里插入图片描述