🧭 Concept Data sources & downloads Raster & rasterio Finding & Fetching Data · step 12 of 14 Updated

STAC Explained: How Satellite Imagery Catalogues Work

Problem statement

You need a cloud-free Sentinel-2 scene over one city, from last July. Traditionally that meant creating an account on a portal, learning its search form, working out which of six product levels you wanted, queuing an order, waiting, and downloading a 1.2 GB ZIP — to use one 10 m band over an area 1/200th the size of the scene.

STAC — the SpatioTemporal Asset Catalog — replaces all of that with a JSON API and a URL you can read directly:

from pystac_client import Client

catalog = Client.open("https://earth-search.aws.element84.com/v1")
search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[-2.30, 53.44, -2.18, 53.52],
    datetime="2026-06-01/2026-08-01",
    query={"eo:cloud_cover": {"lt": 20}},
)
print(search.matched(), "scenes")
5 scenes

No account, no order queue, no download. The problem STAC solves is not storage — it is search. Imagery archives are petabytes; the hard part was always finding the twelve scenes that matter.

Quick answer

STAC is four nested things, and knowing which is which explains every API call:

Concept What it is Analogy
Catalog the entry point, listing collections the library
Collection one dataset, with shared metadata one periodical
Item one scene at one time and place one issue
Asset one downloadable file on an item one page — a band, a thumbnail, metadata

An Item is a GeoJSON Feature with extra required fields. That is the whole trick: the geometry is the scene footprint, the properties are the metadata, and the assets dict points at the actual pixels.

item = next(search.items())
print(item.id, "|", item.datetime.date(), "| cloud", round(item.properties["eo:cloud_cover"], 2))
print("assets:", list(item.assets)[:8])
S2B_30UWE_20260712_0_L2A | 2026-07-12 | cloud 0.06
assets: ['aot', 'blue', 'cloud', 'coastal', 'granule_metadata', 'green', 'nir', 'nir08']
The four levels of a STAC catalogue — catalog, collection, item and asset — with an example at each level.
Four levels. A search returns Items; the pixels live one level further down, in Assets.

Step-by-step solution

1. Find the catalogue and list its collections 

catalog = Client.open("https://earth-search.aws.element84.com/v1")
print(catalog.title)
print([c.id for c in catalog.get_collections()])
Earth Search by Element 84
['sentinel-2-pre-c1-l2a', 'cop-dem-glo-30', 'naip', 'cop-dem-glo-90',
 'landsat-c2-l2', 'sentinel-2-l2a', 'sentinel-2-l1c', 'sentinel-2-c1-l2a',
 'sentinel-1-grd']

Nine collections, and the names encode the product level. sentinel-2-l1c is top-of-atmosphere reflectance; sentinel-2-l2a is surface reflectance, atmospherically corrected. Unless you have a reason to do your own correction, you want L2A.

Note cop-dem-glo-30 in the same catalogue — a global 30 m elevation model, searchable by exactly the same API. STAC is not only for optical imagery.

2. Search on the three things every collection supports 

search = catalog.search(
    collections=["sentinel-2-l2a"],   # what
    bbox=[-2.30, 53.44, -2.18, 53.52],  # where — always WGS84 degrees
    datetime="2026-06-01/2026-08-01",   # when — RFC 3339, `..` for open-ended
)
print(search.matched())
23

bbox is always in WGS84 degrees, regardless of what CRS the imagery is stored in. This is fixed by the specification, and it is one of the few things in geospatial data access that is never ambiguous.

datetime accepts a single instant, a closed range (start/end), or a half-open one (2026-06-01/..).

3. Filter on collection-specific properties

Beyond the big three, every collection publishes its own properties, and you can query them:

search = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[-2.30, 53.44, -2.18, 53.52],
    datetime="2026-06-01/2026-08-01",
    query={"eo:cloud_cover": {"lt": 20}},
)
print(search.matched(), "of 23 scenes under 20% cloud")
5 of 23 scenes under 20% cloud

Eighteen scenes eliminated without downloading a byte. eo:cloud_cover comes from the EO extension — STAC's core is deliberately small, and domain fields arrive through named extensions (eo: for optical, sar: for radar, proj: for projection, view: for geometry).

4. Read the assets, not the item

An Item is metadata. The pixels are in item.assets, one entry per band or file:

red = item.assets["red"]
print(red.media_type)
print(red.href[:80])
image/tiff; application=geotiff; profile=cloud-optimized
https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/30/U/WE/…

profile=cloud-optimized is the part that changes how you work. A Cloud-Optimized GeoTIFF is internally tiled with overviews, and HTTP range requests let a client read one tile without fetching the file. That is why you never need to download the scene.

import geopandas as gpd
from shapely.geometry import shape, box

footprint = shape(item.geometry)
study = box(-2.30, 53.44, -2.18, 53.52)
print(f"study area covered: {study.intersection(footprint).area / study.area:.1%}")
study area covered: 100.0%

bbox search returns scenes whose footprint intersects your box — a scene clipping one corner matches as readily as one covering everything. For anything mosaicked, check coverage per scene and combine, as in merging and mosaicking rasters.

A search funnel narrowing 23 scenes to 5 by cloud cover, then to 1 by footprint coverage, with no pixels downloaded.
Every narrowing here happens on metadata. The first byte of imagery is read after the funnel, not before.

Code examples

Example 1 — choosing the best scene by an explicit rule 

from pystac_client import Client
from shapely.geometry import shape, box

STUDY = (-2.30, 53.44, -2.18, 53.52)


def best_scene(collection, bbox, datetime, *, max_cloud=20, min_coverage=0.95):
    catalog = Client.open("https://earth-search.aws.element84.com/v1")
    search = catalog.search(
        collections=[collection], bbox=list(bbox), datetime=datetime,
        query={"eo:cloud_cover": {"lt": max_cloud}},
    )
    area = box(*bbox)

    scored = []
    for item in search.items():
        coverage = area.intersection(shape(item.geometry)).area / area.area
        if coverage < min_coverage:
            continue
        scored.append((item.properties["eo:cloud_cover"], coverage, item))

    if not scored:
        raise LookupError(f"no scene under {max_cloud}% cloud covering {min_coverage:.0%}")

    scored.sort(key=lambda t: t[0])
    for cloud, coverage, item in scored:
        print(f"  {item.id:28} cloud {cloud:5.2f}%  coverage {coverage:5.1%}")
    return scored[0][2]


chosen = best_scene("sentinel-2-l2a", STUDY, "2026-06-01/2026-08-01")
print("chosen:", chosen.id, chosen.datetime.date())
  S2B_30UWE_20260712_0_L2A     cloud  0.06%  coverage 100.0%
  S2A_30UWE_20260627_0_L2A     cloud  4.31%  coverage 100.0%
  S2B_30UWE_20260622_0_L2A     cloud 11.87%  coverage 100.0%
chosen: S2B_30UWE_20260712_0_L2A 2026-07-12

The selection rule is code, not a click in a portal — so it is reproducible, reviewable and rerunnable next season.

A caveat on eo:cloud_cover: it is a scene-wide percentage. A scene at 0.06% cloud can still have its only cloud sitting on your study area. For serious work, read the scene's cloud mask over your extent instead of trusting the summary number.

Example 2 — reading only your study area from a 10,980-pixel scene 

import os
import time

os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"   # don't list the bucket
os.environ["AWS_NO_SIGN_REQUEST"] = "YES"                  # public data, no credentials

import pyproj
import rasterio
from rasterio.windows import from_bounds

start = time.perf_counter()
with rasterio.open(chosen.assets["red"].href) as src:
    print(f"scene: {src.width} x {src.height} px, {src.crs}, {src.dtypes[0]}")
    print(f"internal tiles: {src.block_shapes[0]}, overviews: {src.overviews(1)}")

    to_scene = pyproj.Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)
    left, bottom = to_scene.transform(STUDY[0], STUDY[1])
    right, top = to_scene.transform(STUDY[2], STUDY[3])

    window = from_bounds(left, bottom, right, top, src.transform)
    red = src.read(1, window=window)

print(f"read {red.shape} in {time.perf_counter() - start:.1f}s "
      f"({red.size / (src.width * src.height):.2%} of the scene)")
scene: 10980 x 10980 px, EPSG:32630, uint16
internal tiles: (1024, 1024), overviews: [2, 4, 8, 16]
read (898, 787) in 5.7s (0.59% of the scene)

Six seconds and half a percent of the scene, over HTTP, with no file on disk. The GDAL_DISABLE_READDIR_ON_OPEN setting matters more than it looks: without it, GDAL lists the whole S3 prefix looking for sidecar files, which can take longer than the read itself.

Example 3 — a time series from the metadata alone 

import pandas as pd

catalog = Client.open("https://earth-search.aws.element84.com/v1")
items = list(catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=list(STUDY),
    datetime="2026-01-01/2026-08-01",
).items())

frame = pd.DataFrame([{
    "date": item.datetime.date(),
    "cloud": round(item.properties["eo:cloud_cover"], 1),
    "platform": item.properties["platform"],
    "id": item.id,
} for item in items]).sort_values("date")

print(f"{len(frame)} scenes; {(frame['cloud'] < 20).sum()} usable")
print(frame.groupby(frame["date"].astype("datetime64[ns]").dt.month)
           .agg(scenes=("id", "size"), usable=("cloud", lambda s: (s < 20).sum()))
           .to_string())
44 scenes; 9 usable
       scenes  usable
date
1           6       0
2           5       1
3           6       1
4           5       0
5           6       2
6           5       2
7           6       3
8           5       0

Forty-four scenes, nine of them usable, all of that established from metadata in a couple of seconds. That monthly table is the honest answer to "can I do a monthly time series here?" — and for January and April, in this climate, the answer is no.

Explanation

Why STAC is just GeoJSON

A STAC Item is a valid GeoJSON Feature. Its geometry is the scene footprint and its properties hold the metadata, with a handful of required fields (datetime, id) and an assets object alongside.

The consequences are practical. Any GeoJSON tool can read an Item. A search response is a FeatureCollection, so it drops straight into GeoPandas:

gdf = gpd.GeoDataFrame.from_features(
    [item.to_dict() for item in items], crs="EPSG:4326"
)
print(gdf[["datetime", "eo:cloud_cover"]].head(3).to_string(index=False))

Footprints become a layer you can plot, intersect and clip like any other. That reuse of an existing format is why STAC spread so quickly.

Why the API is a profile of OGC API Features

STAC's search API is a constrained version of the same standard covered in spatial web services explained — /collections, /collections/{id}/items, limit, bbox, paging by next link. It adds datetime and query on top.

So the truncation discipline from GeoJSON downloaded from an API is empty or truncated applies here too. search.matched() is the STAC name for numberMatched, and search.items() handles paging for you — which is precisely why you should use the client rather than hand-rolled requests.

Why COGs make the download step disappear

A Cloud-Optimized GeoTIFF is an ordinary GeoTIFF with two constraints: the data is stored in internal tiles rather than strips, and reduced-resolution overviews are included. The header sits at the front.

A client reads the header with one small range request, works out which tiles cover the window, and fetches only those. For the scene above that is 0.59% of the file. And because overviews are present, a client that only needs a preview can read the 16× overview instead — a few hundred kilobytes for a whole scene.

This is why STAC catalogues do not offer a "download" button. The file already behaves like a service.

A cloud-optimized GeoTIFF read over HTTP: a header request, then range requests for only the tiles overlapping the study area.
Three small range requests instead of a 1.2 GB download. The file is unchanged; the layout is what makes it possible.

The commonest STAC mistake is not a bad query — it is searching the wrong collection. For Sentinel-2:

  • L1C — top-of-atmosphere reflectance. What the sensor saw, including haze.
  • L2A — surface reflectance, atmospherically corrected, with a scene classification band.

Compute an index like NDVI from L1C and the values are not comparable between dates, because the atmosphere differed. That is a silent error: the arithmetic works, the map looks plausible, and the time series is measuring weather.

Unless you are doing your own atmospheric correction, use L2A.

Edge cases or notes

  • bbox is always WGS84 degrees, even when the imagery is in UTM. Convert your study area before searching, not after.
  • eo:cloud_cover is scene-wide. Use the scene classification or cloud mask asset for a per-pixel answer over your extent.
  • Asset keys are not standardised across collections. Sentinel-2 on Earth Search uses red/nir; other catalogues use B04/B08. Inspect item.assets rather than assuming.
  • Some catalogues require signed URLs. Microsoft's Planetary Computer needs planetary_computer.sign(item) before the asset hrefs work; Earth Search does not.
  • A scene footprint is not a coverage guarantee. Sentinel-2 granules include large nodata areas along orbit edges — a footprint can contain your area while the pixels there are empty.
  • datetime is RFC 3339 and timezone-aware. A bare date is interpreted as UTC midnight, which can shift results by one scene near the boundary.
  • Static catalogues exist too. Not every STAC has a search API; some are just a tree of JSON files on object storage, which you crawl rather than query.

FAQ

What is the difference between a STAC Item and an Asset?

An Item is one scene's metadata — a GeoJSON Feature with a footprint and properties. Assets are the files hanging off it: individual bands, thumbnails, metadata documents.

Do I need an account?

For Earth Search and most public catalogues, no. Some — including Microsoft's Planetary Computer — need a token or a URL-signing step before assets are readable.

Why is bbox in degrees when the imagery is in UTM?

The specification fixes it that way, so one query works across catalogues whose data is in many different CRSs. Convert your study area to WGS84 before searching.

Can I search for elevation data, not just imagery?

Yes. cop-dem-glo-30 and cop-dem-glo-90 are STAC collections on the same catalogue and search identically.

Why should I use L2A rather than L1C?

L2A is atmospherically corrected surface reflectance. Indices computed from L1C are not comparable between dates because the atmosphere differed — a silent error in any time series.

Do I have to download the whole scene?

No, and that is the point. Assets are Cloud-Optimized GeoTIFFs, so a windowed read fetches only the tiles covering your area.

How do I know a scene actually covers my study area?

Intersect item.geometry with your area and compute the fraction. A bbox search returns anything that touches your box, including scenes clipping one corner.