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

news/2024/11/10 16:17:54/文章来源:https://www.cnblogs.com/gaogao-web/p/18336580

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_CLOUDCOMPILE_OPT IDL2e = 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 Highstopmask = MAKE_ARRAY(dimensions, VALUE=1, /INTEGER)FOR N = 0, dimensions[0]-1 DO BEGINFOR M = 0, dimensions[1]-1 DO BEGIN
;      本文只用到bit 3（云）、bit 5（云阴影）进行去云操作
;      其中3和4表示二进制的位置，从右往左数（0开始）所以3和4的索引位置为-4和-5IF dataBit[-4, N, M] EQ 1 OR dataBit[-5, N, M] EQ 1 THEN BEGINmask[N, M] = 0ENDIFENDFORENDFORfile = e.GetTemporaryFilename()maskRaster = ENVIRaster(mask, URI=file)maskRaster.SavemaskedRaster = ENVIMaskRaster(raster[0], maskRaster)e.Data.Add, maskedRasterview=e.GetView()layer=view.CreateLayer(maskedRaster)stop
END

去云结果对比图：

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

十进制值解释含义如下：

代码思路：

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

PRO LANDSAT_MASK_CLOUDCOMPILE_OPT IDL2e = 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)stopmask = MAKE_ARRAY(dimensions, VALUE=1, /INTEGER)FOR N = 0, dimensions[0]-1 DO BEGINFOR M = 0, dimensions[1]-1 DO BEGINIF data[N, M] EQ 55052 OR data[N, M] EQ 23888 THEN BEGINmask[N, M] = 0ENDIFENDFORENDFORfile = e.GetTemporaryFilename()maskRaster = ENVIRaster(mask, URI=file)maskRaster.SavemaskedRaster = ENVIMaskRaster(raster[0], maskRaster)e.Data.Add, maskedRasterview=e.GetView()layer=view.CreateLayer(maskedRaster)stop
END

去云结果对比图：