IDL根据Landsat QA波段去云处理【代码】

IDL根据Landsat QA波段去云处理【代码】

​ landsat QA波段（质量评估波段）是Landsat卫星影像数据中的一个特殊波段，他在Landsat5-9的每个产品中都存在。虽然我们常用的Landsat影像数据有B1-B7波段，但QA波段并不是其中之一。它可以反映出云、云阴影、雪等类别的像素，常常应用在影像处理中对云像素去除。

​ 最近有在写landsat像素去云处理，查了网上许多QA波段值解释说明，发现都是基于二进制的，但IDL不同于GEE的算法，没有>>这种的按位运算符，只能先转成二进制，再自己写算法处理。算法写好后，为了发博客就去查了官网，又发现官网更新的QA波段值解释说明已经更新到了十进制，于是又写了一下根据十进制的去云处理（真的大哭）。

方法一：根据QA给定的二进制值解释进行处理

​ 上面的图片列出了QA波段的每一位所代表的含义，该含义为二进制存储的信息。

​ QA波段的存储方式为十进制，所以转换为二进制值进行判断，下图为某一像素二进制值说明。该像素为云的可能性很大。

代码思路：

1. 读取图像，将十进制的数据转换为二进制格式
2. 云像素识别，并标记，例如（只去除云像素和云阴影像素），为了方便，只使用了bit为3和4的两个为参考，并未加入置信值（confidence)
3. 创建掩膜，对原图像进行掩膜

PRO LANDSAT_MASK_CLOUD   COMPILE_OPT IDL2   e = ENVI()      raster = e.OpenRaster('F:gbytempLC09_L2SP_127031_20220509_20220511_02_T1LC09_L2SP_127031_20220509_20220511_02_T1_MTL.xml')   qaPixelRaster = e.OpenRaster('F:gbytempLC09_L2SP_127031_20220509_20220511_02_T1LC09_L2SP_127031_20220509_20220511_02_T1_QA_PIXEL.TIF')   data = qaPixelRaster.GetData()   dimensions = SIZE(data, /DIMENSIONS)   dataBit = data.toBits() ;  QA Bit    Description                values ;     0      Fill ;     1      Dilated Cloud              1 ;     2      Cirrus                     1 ;     3      Cloud                      1 ;     4      Cloud Shadow               1 ;     5      Snow                       1 ;     8-9    Cloud Confidence           01Low 10Reserved 11 High ;     10-11  Cloud Shadow Confidence    01Low 10Reserved 11 High ;     12-14  Snow/Ice Confidence        01Low 10Reserved 11 High ;     14-15  Cirrus Confidence          01Low 10Reserved 11 High    stop   mask = MAKE_ARRAY(dimensions, VALUE=1, /INTEGER)   FOR N = 0, dimensions[0]-1 DO BEGIN     FOR M = 0, dimensions[1]-1 DO BEGIN ;      本文只用到bit 3（云）、bit 5（云阴影）进行去云操作 ;      其中3和4表示二进制的位置，从右往左数（0开始）所以3和4的索引位置为-4和-5       IF dataBit[-4, N, M] EQ 1 OR dataBit[-5, N, M] EQ 1 THEN BEGIN         mask[N, M] = 0       ENDIF     ENDFOR   ENDFOR       file = e.GetTemporaryFilename()   maskRaster = ENVIRaster(mask, URI=file)   maskRaster.Save   maskedRaster = ENVIMaskRaster(raster[0], maskRaster)    e.Data.Add, maskedRaster   view=e.GetView()   layer=view.CreateLayer(maskedRaster)   stop END 

​ 去云结果对比图：

方法二：根据QA给定的十进制值解释进行处理

​ 十进制值解释含义如下：

代码思路：

1. 读取图像
2. 云像素识别，并标记，例如（只去除云像素和云阴影像素），为了方便，只使用了高置信值云22280、和高置信值云阴影23888为参考，
3. 创建掩膜，对原图像进行掩膜

PRO LANDSAT_MASK_CLOUD   COMPILE_OPT IDL2   e = ENVI()      raster = e.OpenRaster('F:gbytempLC09_L2SP_127031_20220509_20220511_02_T1LC09_L2SP_127031_20220509_20220511_02_T1_MTL.xml')   qaPixelRaster = e.OpenRaster('F:gbytempLC09_L2SP_127031_20220509_20220511_02_T1LC09_L2SP_127031_20220509_20220511_02_T1_QA_PIXEL.TIF')   data = qaPixelRaster.GetData()   dimensions = SIZE(data, /DIMENSIONS)    stop   mask = MAKE_ARRAY(dimensions, VALUE=1, /INTEGER)   FOR N = 0, dimensions[0]-1 DO BEGIN     FOR M = 0, dimensions[1]-1 DO BEGIN       IF data[N, M] EQ 55052 OR data[N, M] EQ 23888 THEN BEGIN         mask[N, M] = 0       ENDIF     ENDFOR   ENDFOR    file = e.GetTemporaryFilename()   maskRaster = ENVIRaster(mask, URI=file)   maskRaster.Save   maskedRaster = ENVIMaskRaster(raster[0], maskRaster)    e.Data.Add, maskedRaster   view=e.GetView()   layer=view.CreateLayer(maskedRaster)   stop END 

​ 去云结果对比图：