Unpacking Landsat bit-packed QA flags

openEO Platform
1

Has anyone ever unpacked the Landsat Collection2 bit-packed QA layers in openEO (bandname BQA)?

SentinelHub provides some utility functionality for this (see below)…
Could be useful to have this also in openEO as a UDP?
@daniel.thiex

2

ok then … lets make a feature request :slight_smile:
@jeroen.dries @benjamin.schumacher

3

Hi Patrick,
as it happens, we started integrating Landsat-8 on CDSE due to a project that needs it. So it might be that we arrive at the very same question quite soon. We’ll then aim for a generic solution.

best regards,
Jeroen

4

Hi Patrick,
I worked with the BQA layer by converting the integer number to a binary representation and then reclassify based on the (for me) important Bits.

udf_reclassify = openeo.UDF("""
from openeo.udf import XarrayDataCube
from openeo.udf.debug import inspect
import numpy as np

def reclassify_number(number):
    # Convert the 16-bit integer to binary representation
    binary_representation = np.binary_repr(number, width=16)
    # Check state of Bit 6 (first of two Cloud Confidence Bits)
    # based on https://www.usgs.gov/landsat-missions/landsat-sr-derived-spectral-indices-pixel-quality-band
    if binary_representation[6] == '1':
        return 1
    else:
        return 0

# Vectorize the function
vectorized_reclassify_number = np.vectorize(reclassify_number)

def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:
    array = cube.get_array()
    temp_array = array.values.astype(np.uint16)
    
    # Vectorize the function to work with NumPy arrays
    results_array = vectorized_reclassify_number(temp_array)

    array.values = results_array
    return cube
""")

ls_bqa = connection.load_collection(
    "LANDSAT8-9_L2",
    temporal_extent=temporal_extent,
    bands=["BQA"],
)
ls_bqa_chip = ls_bqa.filter_spatial(geometries=aoi)
reclassified_bqa = ls_bqa_chip.apply(process=udf_reclassify)

In this case I was looking at Bit 6 of the BQA layer. All pixels with medium or high cloud confidence are set to 1, all others are NA (at least I think so).
Using aggregate_spatial and reducer="mean", I then calculated the share of cloud covered pixels for every date in temporal_extent.

I hope this helps you.
Best regards,
Felix

1 Like
5

Dear Felix,
many thanks for posting your solution. THis is exactly what I was looking for. Hope to find the time to covert this into a UDP sometime soon…