Retrieving corresponding cloud mask for Sentinel-2

Question

I am processing 2A -level Sentinel2 products (directly from Copernicus using the sentinelsat package in Python). Some of them contain clouds, and I want to mask those clouds. Since the cloud masks in the SCL (Scene Classification Layer) product that come with the 2A products from Copernicus are not great, I want to use masks generated by s2cloudless.

To apply the s2cloudless package, I need the corresponding 1C products because the algorithm only works on the 1C product. When I download the product metadata for my bounding box using sentinelsat's SentinelAPI 'query' function, I get about 6000 results (of which exactly half is 2A and the other half 1C). I do not understand however which 2A product corresponds to which 1C product (so essentially; which 2A product was produced from which 1C product). I need this so I know to which 2a product the cloud mask I produce using the 1C product will belong. The 2A products do have a 'level1cpdiidentifier', though these do not correspond to the 1C product identifiers and there is no description I could find that explains how I can use this identifier. Can anyone explain to me how I can best connect the 2A products to their 1C counterparts?

I've also tried directly downloading the pre-calculated s2cloudless masks from Earth Engine and sentinelhub for the L2A products I already have (I would like to do this without interrupting my current workflow). This did not seem to work without getting a paid subscription.

Answer 1

You are analyzing S2 L2A products and you want to use a S2 L1C output. I'd do it in two steps: 1) download S2L2A with sentinelsat. 2) Download S2L1C cloud mask from google earth eninge.

First step (downloading s2 L2A, producttype is the key for getting only L2A):

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

api = SentinelAPI('user', 'pass', 'https://apihub.copernicus.eu/apihub')

footprint = geojson_to_wkt(read_geojson('some_area.geojson'))

products = api.query(footprint,
                     date=(date(yyyy, mm, dd),date(yyyy, mm, dd)),
                     platformname='Sentinel-2',
                     producttype='S2MSI2A',
                     )

api.download_all(products)

Second step (downloading cloud mask from gee):

import ee
# ! pip install geetools if you don't have it
from geetools import batch 

ee.Authenticate()
ee.Initialize()

sentinel = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY")

# coordinates of the polygon to extract
yourarea = ee.Geometry.Polygon(
    [[[x1, y1],
      [x2, y2],
      [x3, y3],
      [x4, y4],
      [x1, y1]]])

# date range
SrtDate = 'yyyy-mm-dd'
EndDate = 'yyyy-mm-dd'

collection = (sentinel
              .filterBounds(yourarea)
              .filterDate(SrtDate,EndDate))

exported_images = batch.Export.imagecollection.toDrive(
    collection=collection,
    folder='S2_Clouds',
    region=yourarea,
    scale=10,
    dataType='int',
    crs= 'EPSG:32XXX', # the CRS of your zone
    maxPixels=10000000000000
)

Then you can merge both products

Answer 2

The names for L1A and L2A products are very similar. Look for the L1C product which has the same aquisition datetime and tile, it should be sufficient

Here is the naming convention of S2 products: https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-2-msi/naming-convention

For example, one L2A derived from the L1C:

S2B_MSIL1C_20210517T103619_N7990_R008_T30QVE_20210929T075738.SAFE
S2B_MSIL2A_20210517T103619_N7990_R008_T30QVE_20211004T113819.SAFE

Only the processing level changes in the name.

Answer 3

Both answers below are very helpful, and helped me understand the process better. I deducted the following viable methods based on aldo_tapia's suggestions. The following works well for me, but if someone has better methods please comment on my code or add an answer:

satimg_cloud = imageCollection.filter(
              ee.Filter.stringEndsWith('system:index', tileid)) \
              .filter(ee.Filter.stringStartsWith('system:index', 
                                               datatakestarttime))

The 'tileid' and 'datatakestarttime' are metadata you can find in your ESA image; the API returns that information (with those exact names) and are unique to your ESA satellite image. You can get the image in this collection using .first(), and the probabilities by selecting the cloud probability band specifically. Alternatively, you can also filter the collection beforehand on cloud probabilities only (for example using ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')).

You could also do the following after filtering the image collection per the first code snippet above, though this will include an additional server call with getInfo():

ee_img = ee.Image("COPERNICUS/S2_CLOUD_PROBABILITY/" + str(satimg_cloud.aggregate_array('system:index').getInfo()[0]))