Making QA masks from a Landsat 8 scene with unpackqa

import rasterio
from matplotlib import pyplot as plt
import numpy as np
import pooch

import unpackqa

This is a Landsat 8 Collection 2 Level 2 product for Path 14, Row 41, May 3, 2013. It's from the example datasets here: https://www.usgs.gov/core-science-systems/nli/landsat/landsat-sample-products

L8_SCENE = 'https://landsat.usgs.gov/sites/default/files/C2_Sample_Data/LC08_L1TP_140041_20130503_20200912_02_T1.zip'
# This is the product code internal to unpackqa only.
L8_qa_product = 'LANDSAT_8_C2_L2_QAPixel'

file_paths = pooch.retrieve(
    L8_SCENE,
    known_hash=None,
    processor=pooch.Unzip()
)
L8_QA_pixel_file = [f for f in file_paths if 'QA_PIXEL' in f][0]

# Load the QA raster as a numpy array.
with rasterio.open(L8_QA_pixel_file) as src:
    img_qa_data = src.read(1)
    
    
def quick_raster_plot(img_data, title='', legend=False, unique_values=[0,1]):
    plt.figure(figsize = (8,8))
    plt.title(title)
    plt.imshow(img_data, interpolation='none')
    if legend:
        plt.colorbar(ticks=unique_values,values=unique_values, fraction=0.03)

First we'll make a mask for a single flag, the "Cloud" flag (bit 4). See the flag details here https://sdtaylor.github.io/unpackqa/Landsat.html

cloud_only_mask = unpackqa.unpack_to_array(img_qa_data, 
                                             product=L8_qa_product, 
                                             flags=['Cloud'])

quick_raster_plot(cloud_only_mask, title='Cloud Only Mask')

Along with clouds you probably also want to mask cloud shadows. It's possible to pull both those masks at once and combine them for a single mask.

Here the cloud_and_shadow_mask array has an added axis at the end, with length 2, representing the 2 masks. Next those 2 axis are summed, so that any value > 0 represents a pixel with either a cloud or cloud shadow

cloud_and_shadow_mask = unpackqa.unpack_to_array(img_qa_data, 
                                                   product=L8_qa_product, 
                                                   flags=['Cloud','Cloud_Shadow'])
print('Original Scene Shape: {}'.format(img_qa_data.shape))
print('Cloud and Shadow mask shape: {}'.format(cloud_and_shadow_mask.shape))

Original Scene Shape: (7351, 7551)
Cloud and Shadow mask shape: (7351, 7551, 2)

cloud_and_shadow_mask = cloud_and_shadow_mask.sum(axis=-1) > 0
quick_raster_plot(cloud_and_shadow_mask, title='Cloud and Cloud Shadow mask')

Most QA bands are binary (either 0 or 1) but some have more values. Here the Cloud Confidence QA flag can take the following values.

• 0 : for no confidence level set
• 1 : Low confidence
• 2 : Medium confidence
• 3 : High confidence

These can be plotted as well to show the variation in cloud confidence across the scene.

cloud_confidence = unpackqa.unpack_to_array(img_qa_data, 
                                              product=L8_qa_product, 
                                              flags=['Cloud_Confidence'])

quick_raster_plot(cloud_confidence, title = 'Cloud Confidence Values', unique_values=[0,1,2,3], legend=True)

If you need several masks for different purposes it would be useful to be able to reference them directly.
For example, say you want to mask out cloudy pixels, but you also want the water mask to highlight rivers or lakes.
The unpack_to_dict function allows for easier separation to use flags later.

masks = unpackqa.unpack_to_dict(img_qa_data, 
                                  product=L8_qa_product, 
                                  flags=['Cloud','Water'])

quick_raster_plot(masks['Cloud'], title='Cloud Mask from dictionary')

plt.figure(figsize = (10,10))
quick_raster_plot(masks['Water'], title='Water Mask from dictionary')

<Figure size 720x720 with 0 Axes>