Cloud cover calculation for subset of a Landsat scene

[cullen-molitor]

I would like to calculate the cloud cover over a subset of Landsat scenes. The eo:cloud_cover property is good for the full scene, but I am only taking a small subset (~1 km²). I would like to calculate a new percentage of cloud cover over this small subset of the scene, and possibly to use this new statistic to filter which images I run through our model. This would hopefully allow me to use images with a high cloud cover overall, but clear skies on the subset of the scene I am using. It would also allow me to evaluate how cloud cover impacts the model we are running imagery through. For example, right now a scene might only have 5 percent cloud cover, but the subset we take could be much higher and it would be nice to keep track of this information.

I have tried digging into the cloud_qa but it appears to only be for Landsat 4-7. I then looked into using the qa_pixel band but I am unsure how to decode the bit index values. This notebook is the closest I could find for guidance, but I haven't been able to get the attributes.

Any recommendation on how to calculate a cloud cover percentage on a subset of an scene? With a little upfront guidance, I would be willing to draft a tutorial notebook for the planetary computer examples repository. Any help would be greatly appreciated!

Here is a small example of what I have done so far to visualize the qa_pixel band:

import stackstac
import pystac_client
import pyproj
import planetary_computer as pc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

lat, lon =  -5, 20

catalog = pystac_client.Client.open('https://planetarycomputer.microsoft.com/api/stac/v1')

search = catalog.search(
    collections=['landsat-c2-l2'],
    intersects=dict(type="Point", coordinates=[lon, lat]),
    datetime="2020-12-01/2021-01-01",
    query={"eo:cloud_cover": {"lt": 25}}
)

items = pc.sign(search)

stack = stackstac.stack(items, assets=["qa_pixel"])

x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)

buffer = 2000  # meters

aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]

image = aoi.data.squeeze().compute()

plt.figure(figsize=(10,10))
im = plt.imshow(image, interpolation='none')
values = np.unique(image.ravel())
colors = [im.cmap(im.norm(value)) for value in values]
patches = [mpatches.Patch(color=colors[i], label=f"{values[i]}") for i in range(len(values))]
plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0. )
plt.show()

[TomAugspurger]

https://stackstac.readthedocs.io/en/latest/examples/gif.html#Mask-cloudy-pixels-using-the-QA-band has a nice example of this.

I think the key sections are to

1. Load the qa_pixel as unit16 values

qa = stackstac.stack([items[0].to_dict()], assets=["qa_pixel"], dtype="uint16", fill_value=0).isel(band=0)
qa

2. Make a cloud mask for the values you care about

mask_bitfields = [1, 2, 3, 4]  # dilated cloud, cirrus, cloud, cloud shadow
bitmask = 0
for field in mask_bitfields:
    bitmask |= 1 << field


bad = qa & bitmask  # just look at those 4 bits
mask = bad == 0

3. Use that mask on your other bands

rgb = stackstac.stack([items[0].to_dict()], assets=["red", "green", "blue"])
good = rgb.where(good)

This would make for a nice addition to our example notebook.

[cullen-molitor]

Thanks, this is exactly the type of documentation I was looking for!