Searching STAC Catalogs for Sentinel-2 with pystac-client

A STAC search is cheap, which is exactly why it is easy to get wrong: a query that returns 4,000 items when you wanted 40 costs nothing to run and quietly wastes the rest of the afternoon in the stacking step. This guide covers building queries that return the scenes you actually meant for a forest study area — the geometry, the property filters, the paging, and the asset signing — as the entry point to Cloud-Native Raster Workflows for Forestry, part of Ecological GIS Data Foundations in Python.

When to use each search parameter

The bbox versus intersects distinction matters more in forestry than in most fields, because ownerships are often long and thin. A riparian corridor 200 km long and 400 m wide has a bounding box covering an entire region, so a bbox search returns every scene over that region while intersects returns only the handful the corridor actually crosses.

Minimal reproducible example

One connection, one search, one list of items. Passing the geometry rather than the bounds is the single most valuable change from the naive version.

import json

import geopandas as gpd
import pystac_client
from shapely.ops import unary_union

CATALOG = "https://earth-search.aws.element84.com/v1"


def find_scenes(stands_path: str, start: str, end: str,
                max_cloud: int = 40, simplify_m: float = 500.0):
    """STAC items covering a stand layer, filtered by date and scene cloud cover."""
    stands = gpd.read_file(stands_path)
    # STAC geometries are WGS84; simplify first so the request body stays small.
    geom = unary_union(stands.to_crs(4326).geometry.values)
    geom = geom.simplify(simplify_m / 111_000.0).buffer(0)

    client = pystac_client.Client.open(CATALOG)
    search = client.search(
        collections=["sentinel-2-l2a"],
        intersects=json.loads(gpd.GeoSeries([geom], crs=4326).to_json())["features"][0]["geometry"],
        datetime=f"{start}/{end}",
        query={"eo:cloud_cover": {"lt": max_cloud}},
        sortby=[{"field": "properties.datetime", "direction": "asc"}],
    )
    return list(search.items())


items = find_scenes("stands.gpkg", "2026-04-01", "2026-09-30")
print(f"{len(items)} scenes")

search.items() pages transparently, so the list comprehension may issue several requests. That is usually what you want, but on a wide query it can return tens of thousands of items — set max_items while exploring, and remove it only once the filters are tight.

Simplifying the geometry before sending it is not cosmetic. Stand boundaries digitised from survey data routinely carry tens of thousands of vertices, and some STAC servers reject a request body above a size limit while others accept it and take a very long time.

Filter on properties, not after the fact

Every property in an item can be filtered server-side. The two that matter most for forestry are scene cloud cover and the MGRS tile id — the latter because pinning the tile keeps a time series on one grid and avoids the seam between adjacent tiles’ differing processing.

def find_scenes_one_tile(tile: str, start: str, end: str, max_cloud: int = 30):
    """Constrain a time series to a single MGRS tile for grid stability."""
    client = pystac_client.Client.open(CATALOG)
    search = client.search(
        collections=["sentinel-2-l2a"],
        datetime=f"{start}/{end}",
        query={
            "eo:cloud_cover": {"lt": max_cloud},
            "grid:code": {"eq": f"MGRS-{tile}"},
        },
    )
    return list(search.items())

Newer servers prefer CQL2, which expresses the same filter with more capable syntax and supports combinations the simple query extension cannot:

def find_scenes_cql2(bbox, start: str, end: str):
    """The same filter in CQL2, which most current STAC APIs prefer."""
    client = pystac_client.Client.open(CATALOG)
    search = client.search(
        collections=["sentinel-2-l2a"],
        bbox=bbox,
        datetime=f"{start}/{end}",
        filter_lang="cql2-json",
        filter={
            "op": "and",
            "args": [
                {"op": "<", "args": [{"property": "eo:cloud_cover"}, 30]},
                {"op": ">", "args": [{"property": "view:sun_elevation"}, 25]},
            ],
        },
    )
    return list(search.items())

The sun-elevation filter is worth knowing about at high latitudes. Winter scenes at low sun angles carry deep terrain shadow that no cloud mask flags, and in mountainous forest they will corrupt an index time series far more than a partly cloudy summer scene would.

Parameter reference

Sign the assets before reading

Some catalogues hand out URLs that need a signature or a session. The Microsoft Planetary Computer signs per item; requester-pays buckets need credentials in the GDAL environment. Getting this wrong produces a 403 at read time, long after the search looked fine.

import planetary_computer


def signed_hrefs(items, assets=("B04", "B08", "SCL")) -> list[dict]:
    """Return signed asset URLs per item for the Planetary Computer catalogue."""
    out = []
    for item in items:
        signed = planetary_computer.sign(item)
        out.append({
            "id": item.id,
            "datetime": item.properties["datetime"],
            "cloud": item.properties.get("eo:cloud_cover"),
            **{a: signed.assets[a].href for a in assets if a in signed.assets},
        })
    return out

For the Element 84 Earth Search catalogue the hrefs are already public and no signing is needed — which is a reason to prefer it for reproducible published work, since a signed URL expires and a provenance record containing one is worthless a day later. Record item ids, never signed hrefs.

Expected output and verification

Three checks turn a search from “it returned something” into a result you can build on:

from collections import Counter


def audit_search(items) -> dict:
    """Sanity-check a STAC result set before spending time stacking it."""
    dates = sorted(i.properties["datetime"][:10] for i in items)
    clouds = sorted(c for c in (i.properties.get("eo:cloud_cover") for i in items)
                    if c is not None)
    tiles = Counter(i.properties.get("grid:code") or i.properties.get("s2:mgrs_tile")
                    for i in items)
    return {
        "n_items": len(items),
        "date_range": (dates[0], dates[-1]) if dates else None,
        "months_covered": len({d[:7] for d in dates}),
        "median_cloud": clouds[len(clouds) // 2] if clouds else None,
        "tiles": dict(tiles),
    }


report = audit_search(items)
assert report["n_items"] > 0, "no scenes — check the geometry CRS and the collection id"
assert len(report["tiles"]) <= 4, f"scenes span {len(report['tiles'])} tiles; pin one for a time series"
assert report["months_covered"] >= 3, "sparse coverage — relax the cloud threshold"

The tile assertion catches the failure that hurts most later. A study area straddling an MGRS boundary returns scenes from two tiles at the same date, and stacking them naively produces two entries per date with a hard edge where they meet.

That figure is the argument for a permissive scene-level threshold. Scene cloud cover describes the whole 110 km tile, most of which is not your forest, so rejecting a scene at 40% throws away imagery that may be entirely clear over the stands. Filter loosely here and do the real work per pixel, as described in Masking Clouds with the Sentinel-2 Scene Classification Layer.

Common pitfalls

• Sending a projected geometry. STAC expects WGS84. A geometry in UTM is accepted as degrees, lands off the coast of Africa, and returns nothing.
• A cloud threshold that is too strict. Scene-level cloud cover is a whole-tile statistic. Filtering at 10% discards scenes that are perfectly clear over a small ownership.
• Forgetting collections. Some servers then search everything, mixing Landsat and Sentinel items whose band names and scaling differ.
• Storing signed URLs in a provenance record. They expire. Store item ids and the catalogue URL, and re-sign when reading.
• Assuming sortby works. Support varies. If ordering matters, sort the returned list in Python rather than trusting the server.

Frequently Asked Questions

Related

• Cloud-Native Raster Workflows for Forestry — the wider pattern this search begins
• Reading Cloud Optimized GeoTIFFs with rasterio Windows — what to do with the asset URLs a search returns
• Chunked Raster Processing with xarray and dask — turning an item list into a lazy stack
• Calculating NDVI from Sentinel-2 with rasterio — the analysis these scenes usually feed

Up: Cloud-Native Raster Workflows · Ecological GIS Data Foundations in Python