Improved searching for S2 imagery from Element 84 STAC sources hosted in AWS

[iferencik]

Context

The data source is described at:
https://github.com/Element84/earth-search?tab=readme-ov-file#sentinel-2-l1c-and-l2a-sentinel-2-l1c-and-sentinel-2-l2a

Within rapida, the goal is to quickly search for imagery that intersects a specific project polygon and then find the most suitable imagery (little clouds inn the specified date-time range).

While this sounds simple the practice has shown this to be a complex task that can be approached in several ways.

Approach

• Find candidate S2 tiles

  • generate a list of MGRS 100K tiles that intersect the geom
  • group by 5 and do stac search for each tile asynchronously to identify candidates
    • extend (by one week in every side of the range) temporal interval in case there are not enough candidates to fully cover each MGRS tile
  • order in space
  • spatial optimization + nodata coverage
    • should be skipped by using async approach by 5 tiles
    • nodata could be moved earlier (stac search by tile)
  • select the candidates base on cc, nodata and time, highest weigh has cloud cover, then time the nodata
    • some tiles can have nodata pixel, in this case more tiles are necessary to join them

• Identify the best candidate for each tile considering time window, cloud cover and nodata coverage

  • Sentinel2 MGRS Item class to
    • accept more the one stac item/json S2MGRSItem (s2id=None, grid_cand_ids=[], )
    • create a mosaic aligned perfectly to MGRS

[Thuhaa]

Approach

1. With the project polygon supplied
2. Get all the tiles that intersect with the bbox. All the tiles from earth search.
3. From these tiles, quickly search for all unique tiles that intersect and from the unique, get the best ones in terms of cloud cover, time range and general differences

Details.

# the step is a distance in degrees to move by to get random points within the bbox
step = 0.5
minx, miny, maxx, maxy = aoi.bounds
    while x <= maxx:
        y = miny
        while y <= maxy:
            mgrs_code = m.toMGRS(y, x, MGRSPrecision=0)
            candidate_tiles.add(mgrs_code)
            y += step
        x += step

This code will give us the codes of the MGRS grid that we have and that can be used to filter in pystac client using an approach like follows

search = client.search(
        collections=["sentinel-2-l1c"],
        query={
            "mgrs:utm_zone": {"eq": 32},
            "eo:cloud_cover": {"lt": 1},
            "mgrs:latitude_band": {"eq": "V"},
            "mgrs:grid_square": {"eq": "KP"},
        },
        max_items=100,
        datetime="2022-01-01/2022-03-28"
    )

[Thuhaa]

1. Start with all 100km tiles in a bbox and time range
2. Group by unique MGRS tile (grid:code in the Sentinel-2 STAC metadata)
3. Within each tile group, rank candidates by: Lowest cloud cover (eo:cloud_cover)
4. The s2:datatake_id shows scenes captured in the same pass. This is likely to show the least differences as you move from a tile to another as the time range is short - prefer items with similar s2:datatake_id if cloud cover percentage is not over the threshold, otherwise, look for the ones close to the selected datetime.

After some research, I found out that while it is simple enough to get non duplicated sentinel tiles, at first, there's no guarantee that the footprint of the tile covers all the area. Example from these 2 items with EXACTLY similar MGRS grid code taken 3 days apart:

1. https://earth-search.aws.element84.com/v1/collections/sentinel-2-l1c/items/S2B_21LXE_20250626_0_L1C
2. https://earth-search.aws.element84.com/v1/collections/sentinel-2-l1c/items/S2B_21LXE_20250629_0_L1C

Using the tileinfo.json has a property dataCoveragePercentage which shows what percentage of the tile has data.