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

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

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

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

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…