今天我们介绍一下Landsat的质量标记波段，也就是QA波段，这个波段标记了Landsat云及其阴影的信息，可以轻松的帮助我们实现云及其阴影的去除。其中，Landsat的云及阴影检测使用的是Fmask算法。 首先我们要了解Landsat的QA波段像元值都代表什么意思。这里我们以Landsat-8为例。
上图是GEE里面对QA波段不同比特代表含义讲解，这里我们截取其中前9比特的含义。可以看出3比特和4比特分别代表云及其阴影，那么比特代表什么意思呢？其实这个比特就是我们把QA波段的值从十进制转化为二进制，3比特就是从右到左第4个数（因为是从0开始的）。这就意味着，我们把QA波段转化为二进制后，如果第四位数是1，那么这个位置对应的像元就是云，如果第五个位置是1，这个像元就是云阴影。
numpy.bitwise_and
numpy.bitwise_and这个函数是处理QA波段非常好用的函数。我们需要识别云像元，就需要把二进制值在第四个位置为1的像元找出来。bitwise_and可以把输入的两个数或者数组都自动转换为二进制，并进行比较，只有两个二进制数某个位置都为1，才会输出1.比如说两个数转换为二进制以后分别为100100100和100，那么输出的值就会是100，如果是100100100和100100000，输出值会是100100000。
那么去云的话，我们就用QA和8（二进制为1000）作为bitwise_and的输入参数，输出值如果是8（二进制为1000）话，就是云，云阴影就把8替换为16。
下面我们看一下具体代码：
代码语言：javascript复制import numpy as np
from osgeo import gdal
out_path='D:DATAcloud_mask.tif'
QA_band = gdal.Open(r'D:DATALC09_L2SP_055016_20220703_20220705_02_T1LC09_L2SP_055016_20220703_20220705_02_T1_QA_PIXEL.TIF')
red_band = gdal.Open(r'D:DATALC09_L2SP_055016_20220703_20220705_02_T1LC09_L2SP_055016_20220703_20220705_02_T1_SR_B4.TIF')
in_band = red_band.GetRasterBand(1)
x_size = in_band.XSize
y_size = in_band.YSize
QA_array = QA_band.ReadAsArray().astype(np.float)
red_array = red_band.ReadAsArray().astype(np.float)
#识别云及阴影位置
cloud_mask=(np.bitwise_and(QA_array,8)==8)
cloud_shadow_mask=(np.bitwise_and(QA_array,16)==16)
#对波段进行去云及阴影
out_red_array=np.where((cloud_mask | cloud_shadow_mask),np.nan,red_array)
driver = red_band.GetDriver()
outimg = driver.Create(out_path, x_size, y_size, 1, gdal.GDT_Float32)
outimg.SetGeoTransform(red_band.GetGeoTransform())
outimg.SetProjection(red_band.GetProjection())
outBand = outimg.GetRasterBand(1)
outBand.WriteArray(out_red_array)
outBand.FlushCache
outimg = None
下面是红波段去云的效果，
原始图像
去云后的红波段
这样去云和阴影就完成拉！
