Attaching Model Metadata to STAC Items

Transform the footprint into WGS84, then record who made the prediction and from what:

import pystac
from rasterio.warp import transform_bounds

bbox = list(transform_bounds(src.crs, "EPSG:4326", *src.bounds))   # degrees, always

item = pystac.Item(id="landcover_36NYF_2026", bbox=bbox, geometry=bbox_to_polygon(bbox),
                   datetime=processed_at, properties={
                       "model:name": "landcover", "model:version": "3.1.0",
                       "model:contract_sha": contract_sha,
                       "proj:epsg": src.crs.to_epsg(),
                   })

A prediction without provenance is a picture. This page belongs to exporting and serving model outputs in Raster Machine Learning & Model Inference.

What Someone Needs to Trust a Prediction

Six months after a run, the questions that arrive are always the same: which model made this, what did it read, when, and how good was it. Each maps onto a field, and each is unanswerable if it was not recorded at write time.

Four questions, four groups of fields Identity fields say which model and version produced the output and hash its input contract. Input fields list the source scene identifiers and the date range they cover. Quality fields carry an accuracy summary and a link to the assessment. Access fields hold the asset links and the projection information a consumer needs to read the file. What the item records identity model:name, model:version model:contract_sha "which model made this" inputs derived_from item ids start and end datetime "what did it read" quality overall accuracy, method link to the full report "how good is it" access — assets, media types, roles, proj:epsg, proj:transform the bbox and geometry are WGS84; the projection fields say what the asset is actually in Every field here is cheap to write and impossible to reconstruct later.

Environment & Setup

Package Version pin Used for
pystac >=1.10 Item construction and validation
rasterio >=1.3.0 Bounds, CRS and transform from the asset
shapely >=2.0 The footprint geometry
jsonschema >=4.0 Offline validation of the item against the spec 
pip install "pystac[validation]>=1.10" "rasterio>=1.3.0" "shapely>=2.0"

Complete Working Example

import hashlib
import json
from datetime import datetime, timezone
from pathlib import Path

import pystac
import rasterio
from rasterio.warp import transform_bounds
from shapely.geometry import box, mapping


def contract_hash(contract_path: str) -> str:
    return hashlib.sha256(Path(contract_path).read_bytes()).hexdigest()[:16]


def build_prediction_item(class_path: str, prob_path: str, *, item_id: str,
                          model_name: str, model_version: str,
                          contract_path: str, source_ids: list[str],
                          accuracy: dict | None = None) -> pystac.Item:
    with rasterio.open(class_path) as src:
        # STAC bbox and geometry are ALWAYS WGS84, whatever the asset CRS is
        bbox = list(transform_bounds(src.crs, "EPSG:4326", *src.bounds, densify_pts=21))
        epsg = src.crs.to_epsg()
        transform = list(src.transform)[:6]
        shape_hw = [src.height, src.width]

    item = pystac.Item(
        id=item_id,
        geometry=mapping(box(*bbox)),
        bbox=bbox,
        datetime=datetime.now(timezone.utc),
        properties={
            "model:name": model_name,
            "model:version": model_version,
            "model:contract_sha": contract_hash(contract_path),
            "derived_from": source_ids,
            "proj:epsg": epsg,
            "proj:transform": transform,
            "proj:shape": shape_hw,
            "processing:level": "L4",
        },
    )
    if accuracy:
        item.properties["accuracy:overall"] = accuracy["overall"]
        item.properties["accuracy:method"] = accuracy["method"]
        item.properties["accuracy:assessed"] = accuracy["date"]

    item.add_asset("class", pystac.Asset(
        href=class_path, media_type=pystac.MediaType.COG,
        roles=["data"], title="Predicted class",
        extra_fields={"classification:classes": [
            {"value": 1, "name": "forest"}, {"value": 2, "name": "cropland"},
            {"value": 3, "name": "built"}, {"value": 4, "name": "water"},
        ]}))
    item.add_asset("probability", pystac.Asset(
        href=prob_path, media_type=pystac.MediaType.COG,
        roles=["data", "quality"], title="Class probability, scaled 1-255"))

    for sid in source_ids:
        item.add_link(pystac.Link(rel="derived_from", target=sid,
                                  media_type=pystac.MediaType.JSON))
    return item

densify_pts=21 in the bounds transform is a small detail with real consequences. A projected rectangle is not a rectangle in degrees; transforming only the four corners cuts the curved edges and produces a bbox that excludes slivers of the actual footprint, so a catalog search along that edge misses the item.

Versioning and Updates

New versions alongside, never in place Three items for the same area and period exist, one per model version, each with its own assets. A latest alias points at the newest. Consumers who pinned version two keep reading exactly what they validated against, while new consumers follow the alias and get version three. One area, three versions, one alias …_v1 superseded, still readable …_v2 pinned by two consumers …_v3 current latest alias Overwriting v2's assets in place would silently change numbers someone already published.

Encoding the model version in the item id rather than only in the properties is what makes this work. A consumer who needs stability pins the versioned id; one who wants the current best follows the alias or searches on model:version. Nothing is ever overwritten, so a result that was reproducible last year stays reproducible.

The cost is storage, and it is usually modest: a cleaned class raster compresses to a couple of megabytes, so keeping three versions of a national product costs less than one scene of the imagery it was derived from. Where storage genuinely matters, expire old versions on a documented schedule rather than by overwriting, so consumers get a deprecation rather than a surprise. The catalog-side mechanics of finding these items again are covered in querying STAC catalogs programmatically.

Verification

item.validate()                       # raises on a spec violation

assert -180 <= item.bbox[0] < item.bbox[2] <= 180, "bbox not in degrees"
assert -90 <= item.bbox[1] < item.bbox[3] <= 90
assert item.properties["model:contract_sha"], "no contract hash recorded"
assert item.properties["derived_from"], "no source items recorded"
assert all(a.href for a in item.assets.values()), "asset without an href"
print(json.dumps(item.to_dict()["properties"], indent=2))
Why the bbox transform needs densified edges A rectangular footprint in a projected CRS becomes a slightly curved quadrilateral in degrees. Transforming only the four corners produces a bounding box that cuts the bulging edges, so parts of the real footprint fall outside the declared bbox. Densifying the edges before taking the bounds includes them. Four corners versus a densified edge corners only misses the bulges densified contains the footprint The gap is small in degrees and large enough to lose an item from an edge-of-area search.

item.validate() catches the structural errors; the assertions catch the semantic ones the schema cannot see. The bbox range check in particular catches the single most common mistake — projected metres left in the bbox — which validates fine against the schema and places the item somewhere off the coast of West Africa, where a surprising number of orphaned geospatial records end up.

Common Errors

STACValidationError about a missing datetime

An item needs either datetime or both start_datetime and end_datetime. For a composite spanning a season, set the range and leave datetime as None.

The bbox is in the asset’s projected CRS. Transform it with transform_bounds and a densified edge, as in the example.

The same item id is generated twice for different areas

The id encodes only the model and the date, not the tile. Include the tile or area identifier as well as the model version, so an id is unique across everything the pipeline will ever produce; a collision silently replaces one product with another when the items are written into a catalog.

Consumers read the wrong band meanings

The class codes live only in your head. Put them in the asset’s classification:classes field and in the raster’s own tags, as described in reading and writing GDAL tags and band descriptions.

Frequently Asked Questions

Q: What is the minimum provenance a prediction item needs? Model name and version, a hash of the input contract, the identifiers of every source scene, and the processing datetime. With those four a run can be reproduced; without any one of them it cannot.

Q: Why must the bbox be in WGS84? Because the specification says so, and every catalog search relies on it. The asset can be in any CRS — record that separately in the projection fields — but a bbox in projected metres will either be rejected or silently place the item in the Gulf of Guinea.

Q: Should accuracy figures live in the item? A summary belongs there — overall accuracy, the assessment date and a link to the full report. It is what lets someone searching the catalog decide whether the product is fit for their purpose without opening it.

Q: Does any of this require a running STAC API? No. A static catalog — a tree of JSON files on object storage with a root catalog linking to collections and items — supports the same search tooling for a read-only archive and costs nothing to host. Add an API only when you need server-side filtering over a large, frequently updated holding; until then the static form is easier to version, easier to back up and impossible to take down by accident.