转载自：https://blog.csdn.net/lidahuilidahui/article/details/102765420

03-SNAP处理Sentinel-2 L2A级数据（一）

• 前言
• 云检测
• • 云检测方法
  • 去云方法
  • • 掩膜法
    • • 再谈Sen2Cor插件
      • 叠加掩膜图层
    • 合成法
    • • Sen2Three插件
• 重采样操作
• • 法一：使用常规的栅格重采样工具
  • 法二：使用专门的Sentinel-2重采样工具
• 超分辨率合成操作
• • 大杀器----Sen2Res
  • • Sen2Res安装
  • Sen2Res超分辨率合成
  • Sen2Res超分辨合成与重采样对比：
  • • 原始60m的B1波段对比
    • 原始20m的B8A波段对比
    • 假彩色图对比
• 裁剪操作
• • 空间子集裁剪
  • • 通过像素坐标范围确定
    • 通过地理坐标范围确定
  • 通过波段选择裁剪
  • 元数据裁剪
• 镶嵌操作
• • 普通镶嵌操作
  • 多尺度镶嵌
• 去云掩膜操作
• • 波段运算创建掩膜
  • 掩膜操作
• 后语
• 参考文献

（原创文章，转载请注明来源！）

 

前言

本文为《SNAP与snappy》第3篇博文。本文主要介绍Sentinel-2 L2A数据（也就是Sen2cor大气校正后的产品数据）的一些基本处理，包括去云、裁剪、镶嵌、重采样、超分辨率合成等操作。先说一下，这篇博文的操作不是一天能后完成的（我写这篇博客画了两天时间）。

博主的电脑配置如下：
运行内存：16G;
CPU: 8 代 i7
操作系统：Win 10,64位系统；
SNAP版本：V6.0.10
中上配置的笔记本电脑。

如果你想SNAP处理的更快一些的话，可以采用以下措施：尽量使运行内存更大；尽可能使用流程图处理（避免多次的读写操作）；一般地使用SNAP的gpt（命令行批处理工具）会比使用GUI界面操作更快一些；使用固态硬盘；如果你有大内存的情况（或者分块操作也是可以）下可以设置设置多进程或多线程操作（snanppy可能更方便），并行处理数据等；使用Linux系统会比Windows快一些。

云检测

云检测方法

按云量划分光学影像的单个像素，可分为三种情况：

1. 无云像素；
2. 部分有云的像素（混合像元）；
3. 全为云的像素。

云（及其形成的阴影）检测分类与去除也是光学遥感的一个热点方向。其难点主要在混合像元的识别，云与高亮地物的识别（如白色的建筑物等），云阴影（厚云下比较明显）等。

云检测（分类）可分为：

• 阈值法：被归到物理方法中，如晴空下的白云，在影像中反射率很高，利用多光谱影像丰富的光谱信息（根据光谱反射率、亮度、纹理特征等特性）设置一个阈值（也有动态阈值的）即可将其识别出来（高于阈值的部分）。目前动态阈值法还在用。
• 模式识别法：也就是常说的无监督和监督分类法，只不过是从地物分类变为了云分类。分类方法有：神经网络、支持向量机、聚类、卷积神经网络等机器学习、深度学习算法等。这是目前主流的云检测（分类）方法。
• 多时相、源分析：前面两种方法主要针对单源单时相数据进行，但实际上可以利用多时相多源数据数据进行，往往能获得更好的监测（分类）效果，因为同一地区不太可能长期都有云覆盖，无云的数据往往能提供更好的特征信息。
• 综合法：结合前面所提的两种或以上的综合方法监测方法。

去云方法

掩膜法

其原理为利用云检测（分类）算法的结果创建一个掩膜文件，将光学影像中的云覆盖部分直接掩膜掉，干脆利落，避免烦恼(效果好坏取决于云检测的精度)。

再谈Sen2Cor插件

在本专栏02篇中介绍了用Sen2Cor对Sentinel-2 L1级做大气校正批处理。Sen2Cor做大气校正，已经有了许多博客文章。

但是，实际上Sen2Cor在大气校正时，会对云、云阴影、卷云、雪、水体等进行场景分类（可以看作地物分类，只是云是大气层的物体），这是绝大多数博文未曾提到的，未曾发现这些信息，意味这些信息很难被利用上。

Sen2Cor的主要场景分类流程如下（更多细节请查看Sen2Cor官方文档）：

实际上，如果你仔细对比过Sen2Cor大气校正前后的Sentinel-2 L1C和L2A级数据，会发现Sentinel-2 L2A的数据文件多了许多文件。博主这里打开的是Sen2Cor处理前后的数据集(L1C数据按原始分辨率做大气正)：

L2A(大气校正后的)数据与L1C数据的差别在于：

1. 多了Index Coding文件夹：
   里面是quality_scene_classification属性表文件，记录的是场景分类结果的编码值、数据类型等描述性信息。
   
2. Bands文件夹(这个文件是存放栅格数据的)多了quality文件夹： 里面存放的是质量控制文件: （都可以通过右击，点击property查看属性信息）
   （这里列出栅格尺寸大小是为了让你留意栅格文件的分辨率）

• quality_aot栅格文件（尺寸：10980X10980）: 气溶胶厚度（AOT, Aerosol Optical Thickness）质量控制文件（无单位）,栅格大小为10980X10980，和B2(波段2)栅格大小相同，表明其分辨率为10m;
• quality_wvp栅格文件（尺寸：10980X10980）：水蒸汽(Water Vapour, WVP)质量控制文件（单位:cm)，分辨率为10m。
• quality_cloud_confidence栅格文件（尺寸：5290X5290）：云质量控制的自信度文件（单位：百分比%），栅格大小为5290X5290，和B5(波段5)栅格大小相同,表明其分辨率为20m;
• quality_snow_confidence栅格文件（尺寸：5290X5290）雪质量控制的自信度文件（单位：百分比%）,分辨率为20m；
• quality_scene_classification栅格文件（尺寸：5290X5290）场景质量分类图，含12个类别，对于去云操作而言，需要关注的是DARK_FEATURE_SHADOW(地形等)，CLOUD_SHADOW（云阴影）、CLOUD_MEDIUM_PROBA（中概率云）、CLOUD_HIGH_PROBA（高概率云）、THIN_CIRRUS（薄卷云）这几个类别。
  
  之前，有人问我是否可以直接可以用这个场景分类图（尽管很好看）当作地物分类图，答案肯定是不能，一个是因为地物类别是很少的，另一个，这里主要针对的对象不同，不是陆地地表地物，精度很低。不过，你非要用的话，那我也办法。

1. Bands文件夹下少了B10波段栅格文件:
   这主要是因为波段10是卷云波段，需要的大气顶部（TOA，Top of Atmosphere）的反射率而不是地表反射率(BOA, Bottom of Atmosphere)。因此，Sen2Cor不会对这个波段做处理。如果你想用这个波段，你需要从L1C数据中提取。
2. Marks文件夹（栅格掩膜文件）下多了一个scl文件：
   里面存放的是10类掩膜文件，其实是从场景分类图中创建的：
   共10个掩膜文件，文件名是很容易理解的，我不再解释。注意设这里的掩膜文件大小为5490X5490,也就是20m分辨率。

叠加掩膜图层

如果你想将掩膜图像添加到栅格影像上的话，特别注意的是两者的分辨率应是相同的。

由于scl文件夹下大的掩膜文件都是20m分辨率的，所以我们需要打开一个20m大的波段，以便掩膜图层叠加20m的灰度影像上。我选择是波段B8A（20m分辨率的近红外波段），先双击L2A数据的B8A栅格数据文件。

如果你的视图和我的不一样，请点击菜单栏的Window，下拉菜单栏Reset Windows。

点击右侧栏的Mask Manager,将滚动条拉到最后，勾选scl_dark_feature_shadow（可能是云阴影、山体、高楼等形成的阴影，很难说，不过这部分的像素很难利用，分类时很可能产生误差，因此，考虑到实际需要，也可以当作云阴影一起掩膜掉）scl_cloud_shadow, scl_cloud_high_proba, scl_thin_cirrus这四个掩膜文件（不选scl_cloud_medium_proba，这个通常会把许多不是云的像素囊括进来）。

默认的颜色显示在图像上是不明显的，可以自己设置更明显的颜色（尽量避免白色、灰色、黑色，因为这是灰度影像，否则很难看到掩膜图像和原影像的差别）。也可以重新设置透明度的数值（0.5，相当于50%，这里我没有修改）（云设置为浅棕红色，云阴影设置为浅湖蓝色）。
结果如下，还可以，但是：

或许你还想真彩色影像对比看看，看下Sen2Cor的云及其阴影掩膜文件大的精度如何，能否用于掩膜。这很好办！

在弹出窗口中保持默认选项, 点击OK：

结果如下：

点击Window—>Tile Horizontally(横向分屏显示):

默认结果：

查看全图：在导航窗口（Navigation）点击下图中的红色圈的组件

结果如下：

可以看到右上角的云识别效果很好，但是左下角的云似乎识别稍差（特别是云阴影)。需要说明的是可能不同的数据识别效果是不一样。总体上还可以。

抱歉！我忘记告诉你，其实，L1C级数据也是带有掩膜文件的（不过很遗憾的是，在生成Sentinel-2 L1C数据时并没有进行云阴影的检测）。好吧，在这一起展示一下，其实，其带有的掩膜文件有许多。不过，只需关注这二个就好：
Mask文件夹下的opaque_clouds（不透明云层，也就是有一定厚度的云层）、cirrus_clouds（卷云）
展开后可以看到：
opaque_clouds，cirrus_clouds文件下有3种分辨率(10，20，60m)下的掩膜文件

我想做的事是把不透明云和卷云叠加到右侧的真彩色影像上，由于其分辨率是10m，所以，这里选择的10m的掩膜文件（不过实际上，当的影像分辨率为10m时,只能看10m的掩膜文件）。先点击一下右侧的真彩色影像，以确保后续的掩膜文件的确是叠加在其影像上的。
同样的操作，我想我不必重复说了：

结果如下：

可以明显地看到，右上角的云L1C数据的自带掩膜效果不是很好，漏了许多范围，但是左下角部分要优于Sen2cor的掩膜结果。当然，光针对一个数据来说算法好坏，是欠妥的。但是，实际上，不同检测算法在单一数据上的表现确实可以确定的，因为你的的确确看到了其在单一数据上结果的好坏。

不过，Sen2Cor的掩膜也有吻合得不错的情况（这是另一个数据集）：

而生成的L1C自带的检测算法效果很差。总体来说Sen2cor生成的掩膜文件要好于L1C集数据自带的掩膜。
由于不同的波段含有不同的分辨率且掩膜文件也有不同的分辨率，因此，现在还不能应用掩膜，需要在重采样或超分辨率合成之后才能其进行掩膜

合成法

其原理利用多时相的影像合成一幅无云影像，例如：最大NDVI像素法、最好像元法、平均像素值法、加权平均合成法等。视数据多少、影像好坏等条件可以合成周、月、年无云影像。

Sen2Three插件

这是一个第三方插件，我希望，你能把SNAP官方提供的插件全部安装上，其实插件的大小都不大，没必要计较这一点储存空间了。其实，还有许多的第三方插件的（linux系统的选择可能更多一些），甚至你可以自己编写相应的插件（如果你的代码能力很强的话）。SNAP的第三插件功能强大，你不应该置若罔闻！

由于该插件需要多个Sen2Cor的大气校正后的影像做输入，需要做大量运算，并且目前处理近年来的数据2018-2019年可能有时有点bug。所以，后文没有对其进行测试。更多的操作过程，请看Sen2Three官方文档。其主要流程如下：

重采样操作

由于Sentinel-2各波段之间存在三种不同的分辨率，如果不统一为同一分辨率的数据，许多后续操作都无法进行。例如裁剪、一些指数计算、分类等。注意：需要确保重采样后数据的输出路径有足够的硬盘存储空间。

法一：使用常规的栅格重采样工具

点击Raster—>Geometric Operations—>Resampling:简单介绍一下参数，就不做结果演示了。

选项①：按定义的源数据参考波段的分辨率确定输出的分辨率（注意默认选择的B1是60m的分辨率，如果你保持默认，所有的10m, 20m的波段都会被重采样为60m，或许你应该查下Sentinel-2各波段的分辨率）

选项②：按照定义栅格尺寸大小设置分辨率，10m分辨率对应的栅格大小为10980X10980, 20m分辨率对应的栅格大小为5490X5490, 60m分辨率对应的栅格大小为1830X1830

选项③：按照定义像素分辨率输出值，单位为m,常用的10，20，60这三个值。

选项④：上采样（提高分辨率，像素分辨率的数值变小），有最邻近（Nearest）、双线性（Bilinear）、三次卷积（Bicubic）三种重采样方法。

选项⑤：降采样（降低分辨率，像素分辨率的数值变大），有首个像素值法、最小像素值法，最大像素值法，平均像素值法，中位数像素值法。（例如将10m的分辨率的4个栅格降采样为20m分辨率的1个栅格，就取第一个栅格的值；4个栅格中的最小像素值；4个栅格中的最大像素值；4个栅格中的平均像素值；4个栅格中的中位数像素值）

选项⑥：后面这个辅助标记选项（只针对降采样）较为复杂，有First,Min,Max, Minmedian, Maxmedian，这个标志有点像像素索引，我想不到它有何作用，也许别的操作可能需要用到它，一般是不要用到它的，保持默认就行。

选项⑦：创建重采样的金字塔文件，有利于加快图像显示，缺点会增加一些额外的存储空间，保持默认就行。
更详细的介绍，请参看帮助文档Help。

法二：使用专门的Sentinel-2重采样工具

我更推荐这种方法，这个官方专门为Sentinel-2重采样设置的模块，重采样的速度会快一些。（博主的电脑采样默认设置，将1个数据的所有波段（除10m波段外）重采样为10m，花了1个半小时左右，处理后，一个数据由原来的近2G增大到了17多G, 后来我发现有些辅助的栅格数据也被重采样了，所以数据量大了许多，所以需要确保的你硬盘有足够的存储空间）
操作流程：optical—>Geometric—>S2 Resampling Processor
其参数和普通的栅格模块重采样是基本一样的，不再赘述，点击run，执行，然后耐心等待就好：

结果：
可以在重采样的数据集上右击打开一个RGB图（Open RGB Image Window）：
你会发现可以合成系统定义好的RGB图像变多了，这是因为它们各波段具有相同的分辨率，可以合成的缘故，当然，你可以自定义用于合成RGB的R、G、B通道文件，用于合成更多的假彩色合成图；
这里打开两个系统定义的RBG图（一个可以指示农作物生长情况的假彩色合成图，一个真彩色合成图，你可以选择别的RGB合成图研究一下）

超分辨率合成操作

也许重采样能能够解决Sentinel-2波段之间分辨率不同的问题，但是这未必是最优的解决方案，想象一下，你直接将60m的数据重采样为10m，将原本的1个栅格变为36个栅格数据，肯定存在不少偏差。

许多光学影像提供全色影像，用于和多光谱波段融合，以提高其分辨率，然而，Sentinel-2影像并没有全色波段，这一想法，无法实现。好在，欧空局提供一个第三方插件Sen2Res，可以实现20m,60m的各波段的超分辨率合成，合成为10m的波段，合成的效果大大优于重采样的结果。

然而，使用Sen2Res插件的代价也是颇大的，博主使用Sen2Res处理一个Sentinel-2数据集，花了近1天半的时间，花的时间太多。建议只针对云量很少近无云的影像进行或者在需要精细分类时使用这个模块。不过，超分辨率合成效果确实不错。

大杀器----Sen2Res

Sen2Res安装

请参照欧空局STEP第三方插件Sen2Res介绍网站的说明，安装到SNAP上。

我相信，你可以看懂安装步骤的，建议安装文件到路径：C:\Users\用户名.snap\auxdata 下，这个路径是存放是辅助文件例如DEM、坐标系定义文件以及插件模块等，以便统一管理SNAP的插件。（当然别的文件夹路径也是可以的）。后续按上图步骤安装即可。

安装好的结果如下图所示：

Sen2Res插件是由印度理工学院（该学院在印度的地位相当于国内的清华大学，在印度非常牛逼,感兴趣大的话，可以百度一下）的GeoStat项目的成果，得到了欧空局相关专家的认可，同时作为SNAP的开源第三方插件，可以免费下载。核心代码是用C++写的，只是集成为Java的NetBeans模块文件（因为SNAP的模块也是这样的模块设置）。如果你想了解更多关于Sen2Res模块的细节，可以到网址http://nicolas.brodu.net/recherche/sen2res/查看更详细的介绍。

Sen2Res超分辨率合成

导入原始的Sen2Cor校正后得到的L2A级数据，点击Optical—>Sentinel-2 Super Resulotion, 打开Sen2Res超分辨率合成操作面板：

只需要设置输出目录就行（注意需要确保输出路径具有足够的存储空间11G左右）：
文件名后缀为_super_resolved

尽管提供了高级参数设置，但是如果你不懂多线程、多进程操作，Sen2Res超分辨率合成算法的话，请保持默认设置就行：

所以，这背后的操作，仅仅是需要设置一下设置的路径就行，然后便是漫长的等待（1天半呀，好在，博主有两个电脑，要不这一半天就没电脑用了，因为这个插件默认尽可能会占用所有的内存，在其运行后别的进程程序很难再使用）。。。。。。

Sen2Res超分辨合成与重采样对比：

需要放大才能看到其效果，加载影像时可能还需要耐心等待一下,好让SNAP反应一下。

原始60m的B1波段对比

所以你看到明显的区别了吗？（左边的采样的B1波段；右边是超分辨率合成的SRB1波段,前缀SR表示Super Resolution,超分辨率之意）

原始20m的B8A波段对比

假彩色图对比

R: B11;G:B8;B:B2;
普通重采样数据：

超分辨合成数据：
效果比较：

超分辨率合成后的波段其纹理等细节大大优于原来的波段重采样的结果，对比度明显增强（注意，默认打开放大的显示效果，没有调整拉伸参数，所以，你不必怀疑是否拉伸的情况造成的），如果你担心精度的话，可以到Sen2Cor官网查看他们论文上的验证结果。

裁剪操作

接下来会裁剪两个Sentinel-2 L2A数据集，然后将这两个数据集镶嵌起来。我想要的效果是得到两个上海市Sentinel-2 L2A数据集中崇明岛部分部分。
如果你需要的是不规则的裁剪，那么需要借助掩膜工具才能实现，不过这种不规则裁剪最后放出最后出图时才做，在此之前尽可能地做矩形裁剪。

由于各波段都具有相同的操作了，可以进行裁剪操作了。

点击Ratser—>Subset(subset, 子集的意思，裁剪即是取子集):
注意，必须先点击一个数据集，才能点亮子菜单栏的Subset操作。

也许我该说一下它的一些参数：

空间子集裁剪

尽管裁剪的操作很简单，但是仍有许多细节值得探讨的，不排除一些细节很可能影响效果。

通过像素坐标范围确定

直接控制输出的栅格尺寸范围，可以设置X、Y轴方向栅格的起点、步长和终点。此外，你还应该知道图像坐标系是以左上角为原点开始的，向右为x轴，向下为y轴。
（注意到博主的裁剪操作了吗，见上图原始两个数据的栅格范围在x轴向是一致的，因此博主并没有动x轴的栅格范围，避免镶嵌后在x轴生成黑色背景的缺失值。）

通过地理坐标范围确定

如果你知道裁剪范围的经纬度边界的话，也可以使用这种方式。（注意：单位是度，而不是度分秒，如果是度分秒的话，可以转化一下单位）

如果不想自己转化坐标单位为度的话，可以参照下图设置显示带小数点以度为单位大的经纬度坐标。在右下角的红色方框可以查看影像中对应的经纬度。从而确定边界范围。

通过波段选择裁剪

原始的Sentinel-2数据除了B1-B12（无B10波段栅格数据）外,还含有许多辅助的栅格数据，根据需要选择波段，通常只需要选择原始B1-B12(无B10)这些波段就行，其它的都可以不选以减少数据量。但是，如果你想保留原始的掩膜数据的话，必须勾选quality_scene_classification栅格波段，因为Masks文件夹的scl文件下得掩膜数据依赖于这个场景分类栅格图件（博主这里多选了一个云质量自信度栅格文件，因为如果是掩膜去云的话，这个波段也是可以用来创建掩膜的）。

元数据裁剪

抱歉这个参数不能改的。
确认后点OK就行。还有一点必须注意点的是，生成裁剪数据集是在内存中生成的，而不是在硬盘上。
如果你右击生成的裁剪数据集，再点Colse Product,那么是会丢失裁剪后的数据。确保数据保存下载来，你需要右击生成的裁剪数据集，点击Save。

镶嵌操作

由于博主目前只处理了一副超分辨率影像，只能用两幅普通重采样的数据来做测试。
这是博主裁剪好的同一天的两个数据Sentinel-2 L2A从采样后的数据集。

普通镶嵌操作

点击Raster—>Geometric Operations—>Mosaicing:

输入输出参数（I/O Prameters）
点击“+”号添加数据，输入镶嵌后的数据文件名和输出目录：
地图投影定义（Map Projection Definition）

坐标系定义（Custom CRS）：输入的Sentinel-2 L2A数据为正射校正的UTM/WGS投影坐标系，为了以原坐标系一致，最好选择UTM/WGS（Automatic），会自动确定投影带。

正射校正参数(orthorectification):不需要在选正射校正参数,因为Sentinel-2 L2A已经经过正射校正了。
镶嵌边界（Mosaic Border）：原始的边界经纬度不用动，像素分辨率设置为10m(需手动设置)
其它，可以点放大镜看看源数据的范围和位置。

变量和条件（Variales&Conditions）

变量（Variables）：指的是选择的波段，第一个圆圈添加的是源数据的波段（可以选择多个），第二个圆圈添加的是变量，可以通过波段运算表达式（Expression）确定，比如新建一个变量B2_add_B3,Expression可设置为：B2 + B3，这样镶嵌后的数据集就会包含这个变量波段B2_add_B3。

条件(conditions）:这个是逻辑条件设置，可以设置某一波段的逻辑条件，例如B8>3000等（只会对B8波段大于3000。为什么可以是3000？Sentinel-2 L2A大气校正后反射率放大了10000倍，反射率为0-10000的整数），某种意义上这可以实现掩膜操作。

这里没有创建变量和设置条件：

镶嵌结果：

可以看到一个很好的镶嵌效果。注意，如果是不同时间或不同源（Sentinel-2A和Sentinel-2B）的影像镶嵌可能会有色差（这点许多软件都很难解决，尽管从理论上说可以做到镶嵌后无色差，但是，这样会改动原始数据的像素值，某种程度上会影响后面进一步制图、波段运算、反演等操作的结果，或许更好的办法是在制图、波段运算、反演等后面再镶嵌会好点）。

这里是同一天同源的数据，因此，镶嵌效果很好（尽管我这里使用有点云覆盖的影像，不过不重要了，用于测试而已）。

你会发现有每个栅格波段都多了一个加了_count后缀的波段，这个波段用于记录每个栅格位置上的像素数值（范围重叠的位置就会加1，在某个栅格上，三个数据重叠，那么是该栅格数值为3）。
打开B1_count，你会发现范围重叠部分比较亮，因为较大，为2，而重叠部分范围为1，边界黑色部分为0：

如果你打开Pixel Info查看像素值，移动鼠标到对应位置即可显示为看到B1_count各栅格的值。

对比：

多尺度镶嵌

这是SNAP为Sentinel-2等存在不同分辨率的数据而开发的镶嵌模块，这个操作允许你对存在不同分辨率的数据进行镶嵌而无需重采样。意味着你可以直接对大气校正后的Sentinel-2数据进行镶嵌。
点击Raster—>Geometric Operations—>Multi-size Mosaic:
它的参数和普通的Mosaic参数基本一样的，我不再重复。

输入输出参数（I/O Prameters）

地图投影定义（Map Projection Definition）

变量（Variables）

镶嵌后有一个缺点：元数据会丢失，Masks里面的掩膜文件会丢失。不过好在，我们保留了quality_scene_classification波段，可以用它创建云掩膜。

去云掩膜操作

这里去云操作放在了裁剪和镶嵌之后，去云的掩膜操作其是可以放在重采样或超分辨率合成之后进行的，或许你可以试下操作顺序造成的差别。

打开镶嵌后的quality_scene_classification波段，点击颜色操作板（Colour Manipulation）,你会发现为了灰度图（其实，可以从裁剪后的数据集quality_scene_classification波段复制颜色表（或者说是颜色模板也可以））:

打开任意一个裁剪后的数据quality_scene_classification波段保存颜色表(scene_classification.cpd,建议保存在路径：C:\Users\用户名.snap\auxdata\color_palettes，该路径是颜色表存放的目录)
再回到镶嵌后的quality_scene_classification波段，点击颜色操作板（Colour Manipulation），导入刚保存的颜色表文件句就好。

不过没有label值，很遗憾只好查看原裁剪数据颜色表中的编码值（见上一幅图）：我们需要利用DARK_FEATURE_SHADOW，CLOUD_SHADOW，CLOUD_HIGH_PROBA，THIN_CIRRUS，编码值分别为2， 3， 9， 10，来创建栅格掩膜文件，不过，镶嵌后的quality_scene_classification波段的像素值为浮点值，对应的编码值为：2.0，3.0，9.0， 10.0。

波段运算创建掩膜

回到镶嵌后的quality_scene_classification波段，右击，选择Band Maths(或者选择Raster—>Band Maths也可以)

设置波段名（Name）,把Virtual(…)前面的勾去掉（这个选项说的是是否创建虚拟波段，虚拟波段不会存储下来），点击Edit Expression选项设置表达式。

可以看到波段名（Data sources,可以用于选择某个波段）；Constants(表示重要的常数,例如Pi, E等)；Operators(包含算术(加减乘除等运算)、关系（大于、小于等运算）、逻辑运算操作符（与，非运算等）)；Functions（表示函数，包含基本的数学函数，例如绝对值函数、三角函数等）；@为占位符，需用波段名或数字（常数Constants）代替。

从Operators从选择if then else条件语句创建0-1掩膜图像，滚动Data sources的滑动滚动条，以便选择波段quality_scene_classification代替占位符@。

完整的波段运算表达式为（注意红色的提示语句，如果绿色表示语法没问题，红色，则会出现错误）：

还可以是别的波段运算表达式，例如通过大于小于运算创建。

确认无误点击OK。返回的界面再点击OK。

结果展示：
尽管颜色板显示的值不对（3个），但是事实上打开的cloud_shadow_mask波段为只有0和1值得二值图像（就是栅格掩膜文件，0值黑色，1值白色，黑色的部分会被掩膜掉）。不信的话，你可以点下Pixel Info窗口，移动鼠标看看。

掩膜操作

由于SNAP中并没有使用栅格进行掩膜的功能（SNAP提供了矢量掩膜的功能，我迟点再谈谈这个功能），不过，SNAP每个波段都可以设置显示有效像素值。利用这个特点也是可以实现栅格掩膜的。

对于镶嵌后的每个波段（例如波段1），右击选择Properties:

将Valid-Pixel Expression(有效像素表达式)设置为（点击右边红圈的…，可以修改表达式）：cloud_shadow_mask > 0.0
（可以复制一个写好的表达式，以便应用到别的波段上）

可以利用别的波段来创建这个有效像素表达式：

点击Close 关闭即可。
其它波段同理可以这样设置，抱歉，我没有找到更方便的修改属性的方法。

来看下B1对比效果
B1
未修改表达式前（原始有效像素表达式为：B1_count >0.0 , 修改后cloud_shadow_mask > 0.0）

掩膜前：

掩膜后：

白色的强反射的云变少了许多。

看看掩膜后的真彩色合成图：

SNAP中还有一个检测云和掩膜的工具----Idxpix，也许后面还会谈到它及去云的处理。你想了解更多的话，请到SNAP的官方论坛，搜索Idxpix的话题。

后语

本来打算每个月写一篇《SNAP与snappy》专栏的文章，奈何总是各种事情干扰，主要忙于导师项目以及学院的各种杂事。虽然我也想一下子，把所有的东西完成，但是无可奈何，SNAP包含的东西太多了，远非数人数月所能完成或学会。事情总要一件一件地来，路需要一步一步走。

SNAP是一款强大的开源遥感处理软件，有强大的社区、论坛资源，任何一个追求开源、自由科学精神的RSer、GISer等地学人都应关注，学习一下。SNAP软件下载次数已超过30万，SNAP论坛有5300多个讨论主题，以及近40000的贴子，包括Google Engine,等商业公司，DIAS等政府部门单位以及遍布各国的科学家以及学生等都在使用。很遗憾的是，国内对SNAP大的关注以及利用不够，SNAP的好多强大功能以及源码、算法等待挖掘。

如果有兴趣的话，可以加入博主创建的欧空局SNAP处理交流群：665903216（这个群已满人），欧空SNAP处理交流群（二）：1102493346。

抱歉，我又写得的太多了，最后，希望你能有所收获。

参考文献

[1] 侯舒维, 孙文方, 郑小松. 遥感图像云检测方法综述[J]. 空间电子技术, 2014, 11(3):68-76.
[2] http://step.esa.int/main/third-party-plugins-2/sen2cor/sen2cor_v2-8/
[3] http://step.esa.int/thirdparties/sen2three/1.1.0/sen2three-1.1.0.htmldoc/sum.html#
[4] http://nicolas.brodu.net/recherche/superres/index.html
[5] https://forum.step.esa.int/t/cloud-mask-snap-sentinel-2/5298/47