How to retrieve Sentinel images without significant cloud coverage using Python

Question

Background

I am trying to retrieve images of certain areas. My aim is to get images without much cloud on them. I would like to specify the percentage of pixels with clouds on them.

---

What I am doing to tackle the problem

I am using the pystac_client to retrieve images. Install pystac_client and rasterio:

%%shell
pip install pystac-client
pip install rasterio

Import useful modules:

from pystac_client import Client
import rasterio
import matplotlib.pyplot as plt

Specify a bounding box (two points defining a rectangle) on the surface of Earth (ie latitudes and longitudes of the southeasternmost and the northwesternmost points of a rectangle, sides are constant latitudes & longitudes). Let's use a point in Singapore for that:

lat, lon = 1.2987, 103.6549 # a random point in Singapore
bbox={'lonLower':lon-0.1,'latLower':lat-0.1,'lonHigher':lon+0.1,'latHigher':lat+0.1}

Search for Sentinel images, similar to what I have done here:

catalog = Client.open("https://earth-search.aws.element84.com/v0")

mysearch = catalog.search(collections=['sentinel-s2-l2a-cogs'],
                          bbox=[bbox['lonLower'],bbox['latLower'],bbox['lonHigher'],bbox['latHigher']],
                          query =  {"eo:cloud_cover":{"lt":1}},
                          datetime='2020-01-01/2020-12-31', 
                          max_items=10)   

resdict = mysearch.get_all_items_as_dict()

I believe the line query = {"eo:cloud_cover":{"lt":1}} should make sure I am only retrieving images with cloud cover less than 1 percent. I haven't found the exact place in the documentation where the usage of this argument is specified, but in another place, the cloud_cover variable is described as:

The estimate of cloud cover as a percentage (0-100) of the entire scene.

To make sure I am not getting this wrong, I have also tried changing 1 to 0.01, it produced the same result as I present here.

---

Plot a retrieved image

Select the first image, use the B08 band (just for example), plot:

url = resdict['features'][0]['assets']['B08']['href']
src = rasterio.open(url)
plt.imshow(src.read(1), cmap='pink')

Singaporean coastline is clearly recognizable. It seems however, that the cloudy parts of the image is way more than 1 percent.

---

Check the Scene Classification Layer

For the same image, I get the Scene Classification Layer by doing:

url = resdict['features'][0]['assets']['SCL']['href']
src = rasterio.open(url)

Each pixel is assigned an integer value, based on what the layer is according to Sentinel's team. The key is:

(Souce: this, "Classification Mask Generation" section.)

Let's plot the values on a histogram:

plt.hist(src.read(1).flatten())

It is clear that pixels classified as CLOUD_MEDIUM_PROBABILITY (Label 8) and CLOUD_HIGH_PROBABILITY (Label 9) are way more prominent than 1 percent.

---

Possible workaround

It would be possible to retrieve many images, and using the SCL, select the ones which have low cloud coverage. However, this method is pretty slow, and would not scale well if I use it in my real-world application.

---

Question

It seems that the above method did not work to retrieve nearly cloudless images. If there is an easy fix to the above method, I'd be glad to hear, if there isn't, I am open to a completely new solution as well.

How do I retrieve Sentinel images, if I want ot specify the maximum ratio of pixels with clouds?

---

---

EDIT: about radouxju's answer:

I tried the Sentinelsat API. I set the api up:

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from datetime import date

api = SentinelAPI('my_username', 'my_password', 'https://apihub.copernicus.eu/apihub')

Then:

# search by polygon, time, and SciHub query keywords
products = api.query("POLYGON ((-3.25 54.5, -3.25 54.7, -3.45 54.7, -3.25 54.5))",
                     date=('20151219', date(2015, 12, 29)),
                     platformname='Sentinel-2',
                     cloudcoverpercentage=(0, 30))

# download all results from the search
api.download_all(products)

As suggested by the docs. The polygon I gave is just an example. The putput is:

Downloading products: 0% 0/1 [00:00<?, ?product/s] LTA retrieval: 0% 0/1 [00:00<?, ?product/s]

Then nothing happens. Upon closer examination, the found that the website says:

Copernicus Open Access Hub no longer stores all products online for immediate retrieval. Offline products can be requested from the Long Term Archive (LTA) and should become available within 24 hours. Copernicus Open Access Hub’s quota currently permits users to request an offline product every 30 minutes.

So I think this solution does not work for near-instant data retrievel, unfortunately.

Answer 1

I think the thing that is tripping you up, is the sentinel:valid_cloud_cover property. Unfortunately I didn't find any documentation on where exactly that flag comes from, but it acts as you would expect.

All the images you find with your query have a cloud cover of zero, but also sentinel:valid_cloud_cover == False.

If you change your query to also check for valid cloud cover, you will get more sensible results.:

mysearch = catalog.search(collections=['sentinel-s2-l2a-cogs'],
                          bbox=[bbox['lonLower'],bbox['latLower'],bbox['lonHigher'],bbox['latHigher']],
                          query={"eo:cloud_cover":{"lt":10},
                                 "sentinel:valid_cloud_cover": {"eq": True}},
                          datetime='2020-01-01/2020-12-31', 
                          max_items=10)   

There's no returns if you search for a cloud cover less than 1% which isn't too surprising for a tropical location, but with less than 10% you get a few returns.

You might be able to set the parameter sortby to get the least cloudy location first, but the pystac documentation isn't too forthcoming with information on how to use this parameter.

Answer 2

As an alternative solution, you could use sentinelsat for your download. There is an option (-c) to specify the cloud percentage with Sentinel-2 images. You can use it in command line or through its Python API.

Note that this relies on the metadata, so you rely on the quality of the default cloud mask on the Copernicus hub, but this is usually enough for screening images but not perfect: you might need another cloud screening algorithm in difficult areas, which is then always very time consuming.