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

去云结果对比图:

热门相关:战神小农民   至尊凰妃   人族镇守使   无敌天下   隋唐君子演义