IDL根据Landsat QA波段去云处理【附代码】
IDL根据Landsat QA波段去云处理【附代码】
landsat QA波段（质量评估波段）是Landsat卫星影像数据中的一个特殊波段，他在Landsat5-9的每个产品中都存在。虽然我们常用的Landsat影像数据有B1-B7波段，但QA波段并不是其中之一。它可以反映出云、云阴影、雪等类别的像素，常常应用在影像处理中对云像素去除。
最近有在写landsat像素去云处理，查了网上许多QA波段值解释说明，发现都是基于二进制的，但IDL不同于GEE的算法，没有>>这种的按位运算符，只能先转成二进制，再自己写算法处理。算法写好后，为了发博客就去查了官网，又发现官网更新的QA波段值解释说明已经更新到了十进制，于是又写了一下根据十进制的去云处理（真的大哭）。
方法一：根据QA给定的二进制值解释进行处理
上面的图片列出了QA波段的每一位所代表的含义，该含义为二进制存储的信息。
QA波段的存储方式为十进制，所以转换为二进制值进行判断，下图为某一像素二进制值说明。该像素为云的可能性很大。
代码思路：
- 读取图像，将十进制的数据转换为二进制格式
- 云像素识别，并标记，例如（只去除云像素和云阴影像素），为了方便，只使用了bit为3和4的两个为参考，并未加入置信值（confidence)
- 创建掩膜，对原图像进行掩膜
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给定的十进制值解释进行处理
十进制值解释含义如下：
代码思路：
- 读取图像
- 云像素识别，并标记，例如（只去除云像素和云阴影像素），为了方便，只使用了高置信值云22280、和高置信值云阴影23888为参考，
- 创建掩膜，对原图像进行掩膜
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
去云结果对比图：