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:\gbytemp\LC09_L2SP_127031_20220509_20220511_02_T1\LC09_L2SP_127031_20220509_20220511_02_T1_MTL.xml')
  qaPixelRaster = e.OpenRaster('F:\gbytemp\LC09_L2SP_127031_20220509_20220511_02_T1\LC09_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:\gbytemp\LC09_L2SP_127031_20220509_20220511_02_T1\LC09_L2SP_127031_20220509_20220511_02_T1_MTL.xml')
  qaPixelRaster = e.OpenRaster('F:\gbytemp\LC09_L2SP_127031_20220509_20220511_02_T1\LC09_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

​ 去云结果对比图：