Searching Sentinel-2 Scenes with pystac-client

One-sentence answer: search the catalogue with the field boundary as GeoJSON in EPSG:4326, a season-long datetime range and a generous scene cloud filter, then decide scene by scene using the cloud fraction inside the boundary — which is the only cloud number that describes your field.

Context

A season’s imagery pipeline starts with a list. Getting it wrong in either direction is expensive: too strict a filter and a wet fortnight leaves a three-week gap that the pipeline reports as “no data” rather than “we discarded it”; too loose and every downstream stage spends time on scenes whose field is under cloud. The search itself is cheap — it touches an index, not pixels — so the right shape is a wide search followed by a narrow, evidence-based rejection.

This is the entry point for satellite imagery APIs and STAC catalogues, and the accepted items feed the windowed reads in windowed COG reads over HTTP with rasterio.

What a season's search actually yields Four stages of a filter funnel over one temperate season: about thirty acquisitions intersect the field, twenty-two survive a sixty percent scene cloud pre-filter, eighteen remain after duplicate tile acquisitions are grouped, and eleven pass a ten percent cloud test measured inside the field boundary. 30 acquisitions intersect the field 150-day season, 5-day revisit 22 after scene cloud < 60% index-only, no bytes read 18 after tile deduplication field straddling two tiles 11 usable — field cloud < 10% ≈ one every 13 days The gaps are not evenly spaced: clear runs in midsummer, and three-week holes in wet spells — exactly when disease pressure peaks.

Prerequisites

Beyond the parent topic’s stack: pystac-client 0.8.*, a field boundary as GeoJSON in EPSG:4326, and a STAC endpoint for a Sentinel-2 Level-2A collection.

Step-by-step

1. Open the catalogue and check the collection. A renamed collection is the second most common cause of an empty result.

What each stage of the search actually transfers Horizontal bars comparing bytes transferred per stage for one field: a catalogue search over a whole season is under a megabyte because it touches only the index; a classification-band window is about a megabyte and decides acceptance; four reflectance windows are a few megabytes; and a single full band tile is over a hundred megabytes. Catalogue search, whole season ≈ 0.4 MB · index only Classification band, one field window ≈ 1.2 MB · decides accept or reject Four reflectance windows, one date ≈ 3.6 MB · only for accepted scenes One full band tile, one date ≈ 104 MB · what downloading costs Reading the classification band first is why a rejected scene costs about one megabyte instead of four hundred.

2. Search wide — boundary, season, scene cloud under 60%.

3. Group by acquisition datetime to catch a field split across two tiles.

4. Test cloud inside the boundary using the scene classification band, which is one 20 m band rather than four reflectance bands.

5. Report the revisit you actually achieved before building anything on top of it.

PYTHON 
from datetime import date, timedelta
from collections import defaultdict

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.mask import mask as rio_mask
from pystac_client import Client
from shapely.geometry import shape

CATALOG_URL = "https://earth-search.aws.element84.com/v1"
COLLECTION = "sentinel-2-l2a"
CLOUD_CLASSES = (3, 8, 9, 10)          # shadow, cloud medium, cloud high, cirrus


def usable_scenes(boundary_geojson: dict, start: date, end: date,
                  scene_cloud_max: int = 60, field_cloud_max: float = 0.10) -> list:
    """Every acquisition over a field that is actually clear enough to use."""
    client = Client.open(CATALOG_URL)
    collections = {c.id for c in client.get_collections()}
    assert COLLECTION in collections, (
        f"{COLLECTION} not in this catalogue; available: {sorted(collections)[:8]}")

    geom = shape(boundary_geojson)
    assert -180 <= geom.bounds[0] <= 180 and -90 <= geom.bounds[1] <= 90, (
        f"boundary bounds {geom.bounds} are not degrees — reproject to EPSG:4326 first")

    search = client.search(
        collections=[COLLECTION],
        intersects=boundary_geojson,
        datetime=f"{start.isoformat()}/{end.isoformat()}",
        query={"eo:cloud_cover": {"lt": scene_cloud_max}},
        limit=100,
    )
    items = sorted(search.items(), key=lambda i: i.datetime)
    print(f"{len(items)} item(s) after the scene cloud pre-filter")

    # A field near a tile edge appears twice at the same instant — keep the fuller footprint.
    by_time = defaultdict(list)
    for item in items:
        by_time[item.datetime.replace(second=0, microsecond=0)].append(item)
    deduped = [max(group, key=lambda i: shape(i.geometry).intersection(geom).area)
               for group in by_time.values()]
    print(f"{len(deduped)} acquisition(s) after tile deduplication")

    accepted = []
    for item in deduped:
        frac = field_cloud_fraction(item, boundary_geojson)
        if frac <= field_cloud_max:
            accepted.append(item)
        print(f"  {item.datetime:%Y-%m-%d}  scene {item.properties['eo:cloud_cover']:5.1f}%  "
              f"field {frac * 100:5.1f}%  {'accept' if frac <= field_cloud_max else 'reject'}")
    return accepted


def field_cloud_fraction(item, boundary_geojson: dict) -> float:
    """Share of field pixels flagged cloud or shadow, read from the 20 m classification band."""
    with rasterio.open(item.assets["scl"].href) as src:
        proj = gpd.GeoSeries([shape(boundary_geojson)], crs="EPSG:4326").to_crs(src.crs).iloc[0]
        arr, _ = rio_mask(src, [proj], crop=True, filled=True, nodata=0)
    scl = arr[0]
    inside = scl != 0
    if not inside.any():
        return 1.0                       # field outside the footprint — treat as unusable
    return float(np.isin(scl[inside], CLOUD_CLASSES).mean())

Inline verification — check the revisit you got, not the one the mission advertises:

PYTHON 
scenes = usable_scenes(boundary, date(2026, 4, 15), date(2026, 9, 15))
dates = [s.datetime.date() for s in scenes]
gaps = [(b - a).days for a, b in zip(dates, dates[1:])]
print(f"{len(scenes)} usable scenes; median gap {np.median(gaps):.0f} d, worst {max(gaps)} d")
assert len(scenes) >= 6, (
    f"only {len(scenes)} usable scenes in the season — relax the field cloud threshold or "
    "plan for drone flights over the gaps")
assert max(gaps) < 35, f"a {max(gaps)}-day gap will not support in-season decisions alone"

Settings worth being deliberate about

Setting Default here Why it matters
scene_cloud_max 60% A pre-filter over a 110 km tile. Tightening it to 20% saves a few classification-band reads and throws away scenes that are clear over the field, which is the wrong trade in a wet season
field_cloud_max 10% The decision that matters. At 20% cloud shadow starts entering the index and reads as crop stress; at 0% a single flagged headland pixel discards an otherwise perfect scene
Cloud class set 3, 8, 9, 10 Shadow, medium and high cloud probability, and thin cirrus. Dropping class 3 keeps scenes whose shadows will be read as stress; adding class 10 alone catches haze that the others miss
Deduplication key acquisition minute Groups the same overpass appearing in two overlapping tiles. Keying on the item identifier instead treats them as separate scenes and truncates the field
Search limit 100 Page size, not a result cap — items() follows pagination transparently
Season window planting to harvest Searching a calendar year doubles the request count and returns scenes over bare soil that no index will use

Cache the per-field cloud fraction the first time it is computed. It is the only part of the search that reads pixels, and a season over 400 fields is around two hours of I/O that should be paid once rather than on every pipeline run.

Gotchas and edge cases

  • An empty result is almost always a geometry problem. Degrees versus metres, or latitude before longitude. The bounds assertion above catches both in one line, which is worth more than any amount of staring at a search that “should” work.
Diagnosing an empty result in one check A decision diagram for an empty search. If the boundary's coordinates fall outside the range of degrees, the geometry is in the wrong coordinate reference system or axis order, which intersects nothing anywhere. If they are valid degrees, the collection identifier is the next thing to check. Search returned zero items boundary bounds within ±180, ±90? yes Check the collection identifier a renamed collection is the next cause; list the catalogue and compare no The geometry is not in EPSG:4326 UTM metres or latitude-first order intersect nothing, and STAC answers with an empty list An empty result is never an error — the search succeeded and found nothing, which is why the assertion has to be yours.

  • limit is a page size, not a cap. search.items() follows pagination on its own, so setting limit=10 does not return ten scenes — it returns everything, ten at a time. Use max_items if you genuinely want a cap.

  • The same field can appear in two tiles at the same instant. Sentinel-2 tiles overlap by several kilometres. Taking the first item silently truncates the field; the deduplication above keeps the footprint with the larger intersection, and a field split roughly evenly needs both windows mosaicked.

  • Cloud classes are conservative and imperfect. The classification band misses thin cirrus over bright soil and over-flags bright field edges. It is good enough to rank scenes and not good enough to be a mask on its own — that job belongs to cloud masking for agricultural imagery.

  • Reading the classification band for every candidate costs real time. At 22 candidates per field per season and 400 fields it is roughly two hours of I/O. Cache the fraction per item and field so a re-run of the pipeline never recomputes it.

This guide is part of Satellite Imagery APIs & STAC Catalogues for Field Monitoring — see there for asset selection, reflectance scaling and the manifest registration that follow.

Parent section Satellite Imagery APIs & STAC Catalogues for Field Monitoring Related Windowed COG Reads over HTTP with Rasterio