Satellite Imagery APIs & STAC Catalogues for Field Monitoring

Drone flights give you centimetres when someone flies; satellites give you a usable field observation every week or two, for free, over every field you manage, without anyone leaving the office. The output of this pipeline is a per-field, per-date stack of surface reflectance windows — the input to band math and raster algebra and, once several dates accumulate, to temporal aggregation of vegetation indices. This topic sits under farm data platform engineering and assumes the manifest and object-key conventions described there.

The discipline that makes this work at farm scale is refusing to download tiles. A Sentinel-2 granule covers 110 km by 110 km; a field covers a few hundred metres. Reading the field window directly out of the remote asset turns a day-long backfill into a coffee break, and it is the difference between a platform you can re-run and one you run once and hope.

Prerequisites

• Python 3.11+, pystac-client 0.8., rasterio 1.3., shapely 2.0., numpy 1.26., geopandas 1.0.*
• A STAC API endpoint for the archive you are using, and credentials if it requires them
• Field boundaries in EPSG:4326 for search, and a metre-based CRS per field for area work — see understanding CRS in precision agriculture
• GDAL HTTP environment variables configured as set out in the section overview — without GDAL_DISABLE_READDIR_ON_OPEN, every open lists the bucket prefix
• Familiarity with band identifiers for the sensor: Sentinel-2 B04 red, B05 red edge, B08 NIR at 10 m, B8A NIR at 20 m — as distinguished in parsing Sentinel-2 vs drone multispectral bands

1. Concept: A Catalogue Is an Index, Not an Archive

A SpatioTemporal Asset Catalog item is a JSON document describing one scene: its footprint as a polygon, its acquisition datetime, properties such as estimated cloud cover and processing baseline, and an assets map of band name to HTTP URL. The pixels are somewhere else. This separation is the whole point — you query a small, fast index to decide which of several million scenes are worth touching, and only then pay for bytes.

A search takes three filters that matter for field monitoring: an intersects geometry (the field boundary, or a bounding box around a farm), a datetime range, and a property query on cloud cover. The result is a list of items ordered by acquisition time, each of which you can accept or reject before reading a pixel.

2. Step-by-Step Implementation

Step 1 — Search by footprint and season

from datetime import date
from pystac_client import Client

CATALOG_URL = "https://earth-search.aws.element84.com/v1"
COLLECTION = "sentinel-2-l2a"


def search_scenes(boundary_geojson: dict, start: date, end: date, max_cloud: int = 60) -> list:
    """Return catalogue items intersecting a field boundary, newest last."""
    client = Client.open(CATALOG_URL)
    search = client.search(
        collections=[COLLECTION],
        intersects=boundary_geojson,
        datetime=f"{start.isoformat()}/{end.isoformat()}",
        query={"eo:cloud_cover": {"lt": max_cloud}},
        limit=100,
    )
    items = sorted(search.items(), key=lambda i: i.datetime)
    assert all(i.geometry is not None for i in items), "catalogue item without a footprint"
    return items

The limit parameter is a page size, not a result cap — search.items() follows pagination transparently. A whole season over one field typically returns 60–75 items before cloud filtering and 8–15 usable ones after.

Step 2 — Decide on cloud over the field, not over the tile

The scene classification layer is a 20 m categorical band; classes 3 (cloud shadow), 8 (cloud medium probability), 9 (cloud high probability) and 10 (thin cirrus) are the ones that ruin an index. Read only the field window of that band and compute the masked fraction:

import numpy as np
import rasterio
from rasterio.mask import mask as rio_mask
from shapely.geometry import shape

CLOUD_CLASSES = (3, 8, 9, 10)


def cloud_fraction_over_field(item, boundary_geojson: dict) -> float:
    """Share of field pixels flagged as cloud or shadow in the scene classification band."""
    href = item.assets["scl"].href
    with rasterio.open(href) as src:
        geom = shape(boundary_geojson)
        proj = (
            __import__("geopandas")
            .GeoSeries([geom], crs="EPSG:4326")
            .to_crs(src.crs)
            .iloc[0]
        )
        arr, _ = rio_mask(src, [proj], crop=True, filled=True, nodata=0)
    scl = arr[0]
    valid = scl != 0
    assert valid.any(), "field window is entirely outside the scene footprint"
    return float(np.isin(scl[valid], CLOUD_CLASSES).mean())

Reading the classification band first is a deliberate economy: it is one 20 m band, so the rejection decision costs a fraction of what reading four reflectance bands would. On a season’s search that ordering roughly halves total bytes transferred.

Step 3 — Read the reflectance windows

def read_field_window(item, band: str, boundary_geojson: dict):
    """Return (array, transform, crs) for one band clipped to the field boundary."""
    import geopandas as gpd

    href = item.assets[band].href
    with rasterio.open(href) as src:
        proj = gpd.GeoSeries([shape(boundary_geojson)], crs="EPSG:4326").to_crs(src.crs).iloc[0]
        arr, transform = rio_mask(src, [proj], crop=True, filled=True, nodata=0)
        crs = src.crs
    assert arr.ndim == 3 and arr.shape[0] == 1, f"unexpected shape {arr.shape} for {band}"
    return arr[0], transform, crs

rasterio.mask with crop=True issues range requests for exactly the internal tiles the window touches. Watch the resolution mix: B04 and B08 are 10 m while B05 and SCL are 20 m, so arrays from different bands do not align. Resample once, explicitly, to the resolution the index needs — never let a broadcast do it implicitly, which is the alignment failure discussed in fixing nodata and NaN propagation in band math.

Step 4 — Apply the processing-baseline offset

def to_reflectance(dn: np.ndarray, item) -> np.ndarray:
    """Convert digital numbers to surface reflectance, honouring the baseline offset."""
    offset = 0.0
    baseline = str(item.properties.get("s2:processing_baseline", "00.00"))
    if float(baseline.replace("N", "") or 0) >= 4.0:
        offset = float(item.properties.get("boa_add_offset", -1000.0))
    refl = (dn.astype("float32") + offset) / 10_000.0
    refl[dn == 0] = np.nan                       # 0 is nodata, not zero reflectance
    finite = refl[np.isfinite(refl)]
    assert finite.size == 0 or (-0.2 < np.nanmedian(finite) < 1.4), (
        f"implausible reflectance median {np.nanmedian(finite):.3f} — offset applied twice?")
    return refl

A doubled offset is the classic symptom: reflectance shifts by 0.1, vegetation NDVI drifts down by two to four hundredths, and the time series develops a step exactly at the baseline change. The assertion above catches it on the first scene rather than in a zone map three months later.

Step 5 — Register the scene

Each accepted scene becomes one manifest row keyed on collection, item identifier and the asset checksum, written in the transaction that publishes the derived window. The pattern is set out in the section overview and implemented in writing idempotent tasks for imagery ingestion.

3. Key Parameters and Tuning

4. Edge Cases and Failure Modes

A field straddling two tiles. Sentinel-2 tiles overlap, so a field near a tile edge appears in two items with the same datetime. Reading either alone truncates the field. Detect the duplicate by acquisition time and either mosaic the two windows or keep only the item whose footprint fully contains the boundary — and assert containment rather than assuming it.

Nodata is zero, and zero is a valid measurement. Sentinel-2 L2A uses 0 as nodata, which collides with legitimately dark pixels only in pathological cases but collides constantly with the arrays produced by clipping to a polygon, where everything outside the field is filled with the nodata value. Mask before you compute anything, never after.

Sun angle and seasonal bias. Reflectance from a low winter sun is not comparable with midsummer reflectance even after atmospheric correction, and the effect is strong enough at high latitudes to fake a green-up. When comparing dates across a season, compare like with like — the same sensor, similar view geometry — and treat any index step that coincides with a satellite change as suspect until proven otherwise.

Baseline changes mid-series. Scenes reprocessed under a new baseline are re-published with the same acquisition date and a different item identifier. Manifest identity keyed on item identifier alone will happily ingest both; key on collection, tile, datetime and baseline, and prefer the newest baseline when duplicates appear.

The catalogue endpoint moves or rate-limits. Public STAC endpoints are best-effort. Wrap searches in the same retry-with-jitter policy used for management-system APIs, cache item metadata locally once fetched, and never make a user-facing request depend on a live catalogue call.

Empty search results are not always empty. A search that returns nothing usually means the geometry was passed in the wrong axis order or the wrong CRS. STAC expects GeoJSON in EPSG:4326 with longitude first; a boundary in UTM metres silently intersects nothing. This is the axis-order and CRS confusion catalogued in debugging CRS and projection errors.

5. Verification and Output Validation

def validate_scene_window(refl: dict[str, np.ndarray], field_ha: float, pixel_m: float = 10.0) -> None:
    """Sanity-check a set of band windows before they are published."""
    shapes = {b: a.shape for b, a in refl.items()}
    assert len(set(shapes.values())) == 1, f"bands are not aligned: {shapes}"

    any_band = next(iter(refl.values()))
    covered_ha = np.isfinite(any_band).sum() * (pixel_m ** 2) / 10_000
    assert covered_ha > 0.5 * field_ha, (
        f"only {covered_ha:.1f} ha of valid pixels for a {field_ha:.1f} ha field — "
        "clipped to the wrong footprint or mostly cloud")

    red, nir = refl["B04"], refl["B08"]
    with np.errstate(invalid="ignore", divide="ignore"):
        ndvi = (nir - red) / (nir + red)
    median = float(np.nanmedian(ndvi))
    assert -1.0 <= median <= 1.0, f"NDVI median {median} outside physical range"
    assert median > 0.05, (
        f"NDVI median {median:.3f} suggests bare soil or a band mix-up — "
        "confirm B08 is NIR and B04 is red for this collection")

The last assertion is worth keeping even though it will occasionally fire legitimately, on a genuinely bare field in April. A swapped red and NIR assignment produces a near-perfect mirror image of a correct index — every value negated — and nothing else in the pipeline notices, because the arrays have the right shape, the right CRS and a plausible range.

Feed the validator a deliberately swapped pair in the test suite and assert it raises. As everywhere else on this site, a gate that has never rejected anything proves nothing.

6. Integration with the Broader Pipeline

Accepted windows go straight into index calculation — calculating NDVI and NDRE with rasterio works on exactly these arrays — and their masks come from the same classification band that automating cloud removal in Sentinel-2 time series uses. Once a season accumulates, the stack is stored as described in cloud-optimized storage for field imagery and reduced to per-field features for management zone classification.

Satellite and drone imagery are complements, not substitutes: satellites establish the temporal baseline cheaply, and a drone flight resolves what a 10 m pixel cannot when the baseline shows something worth flying. Choosing which index to compute at which growth stage is covered in vegetation index selection for crop stages.

Frequently Asked Questions

How many usable scenes should I expect per season? Over a temperate maize or soybean season of roughly 150 days, a five-day nominal revisit yields about 30 acquisitions per field, of which typically 8–15 pass a 10% field-cloud threshold. The distribution is uneven — long clear runs in July, and gaps of three weeks in wet spells — which is why gap-filling and temporal smoothing matter more than raw revisit.

Should I store the raw windows or only the computed index? Store the reflectance windows. They are small once clipped to a field, and every future index — a new red-edge formulation, a soil-adjusted variant with a different L factor — can be recomputed from them. Storing only NDVI locks you out of exactly the improvements this site’s index pages describe.

Is a commercial imagery provider worth it? For daily revisit at 3 m during a critical window, sometimes. But run the free archive first: most operations discover their limiting factor is not revisit but the absence of a pipeline that reliably turns available scenes into decisions, and paying for more scenes does not fix that.

This topic is part of Farm Data Platform Engineering: APIs, Storage & Orchestration — see there for the manifest, storage and scheduling context.

Related

• Searching Sentinel-2 Scenes with pystac-client — a runnable season search with cloud filtering and result inspection
• Windowed COG Reads over HTTP with Rasterio — the GDAL settings and window arithmetic behind a 200 ms field read
• Cloud Masking for Agricultural Imagery — what to do with the pixels the classification band flags
• Cloud-Optimized Storage for Field Imagery — where the accepted windows live and how they are keyed
• Temporal Aggregation of Vegetation Indices — turning an irregular series of accepted dates into a comparable seasonal signal