Spatial Data Archival · Cold Storage

Indexing Archived Imagery Footprints for Fast Discovery

Discovery is the part of an archive that has to be fast even when the data behind it is not. A user searching a decade of imagery expects results in under a second, and delivering that against several hundred thousand scenes means an index built for the query shape — bounded in space, bounded in time, filtered on a handful of properties. This walkthrough is for the engineer building that index, and it covers the choice most archives get wrong: indexing bounding boxes when the footprints are not rectangles.

Bounding Boxes Lie About Imagery

A satellite scene is a rotated quadrilateral; an aerial run is a long thin strip; a mosaic tile with nodata edges covers less than its extent suggests. Indexing the bounding box alone returns scenes whose actual coverage misses the query entirely.

Footprint coverage of the bounding box, by scene shape Three scene shapes with the share of their bounding box they cover: a north-up tile at 100 percent, a rotated scene at 71 percent, and a diagonal flight strip at 23 percent, where a corner query matches the box but not the imagery. mosaic tile covers 100% of its bbox rotated satellite scene covers 71% aerial flight strip covers 23% amber query matches the box, not the imagery

Step-by-Step Procedure

Step 1 — Store the real footprint, not just the extent

The STAC item’s geometry member is the footprint and its bbox is the extent. Populate both: the bbox drives coarse filtering, the geometry decides the answer.

from osgeo import gdal
from shapely.geometry import Polygon, mapping
import numpy as np

def footprint(path: str, simplify_m: float = 5.0) -> dict:
    """Real coverage polygon from the valid-data mask, not the raster extent."""
    ds = gdal.Open(path)
    mask_ds = gdal.Translate("/vsimem/mask.tif", ds, bandList=[1],
                             outputType=gdal.GDT_Byte, scaleParams=[[0, 1, 255, 255]])
    poly_ds = gdal.Polygonize(mask_ds.GetRasterBand(1).GetMaskBand(), None,
                              *_ogr_layer("/vsimem/fp.shp"), 0)
    geom = _largest_ring("/vsimem/fp.shp").simplify(simplify_m)
    return mapping(geom)

Simplify the footprint before storing it. An unsimplified valid-data mask can hold tens of thousands of vertices per scene, which inflates the index without changing any answer at query resolution.

Step 2 — Build a two-stage index

Coarse filtering on a spatial key is fast and cheap; exact intersection is correct and slower. Running them in that order gives both.

-- Stage 1: an H3 cell list per scene, for index-only pruning
CREATE TABLE scene_cells AS
SELECT item_id, h3_cell
FROM   scenes, UNNEST(h3_polygon_to_cells(footprint, 6)) AS t(h3_cell);

CREATE INDEX ON scene_cells (h3_cell);
CREATE INDEX ON scenes USING GIST (footprint);      -- stage 2: exact geometry

-- The query uses both: cells prune, geometry decides
SELECT s.item_id, s.datetime, s.gsd, s.storage_tier
FROM   scenes s
WHERE  s.item_id IN (SELECT item_id FROM scene_cells
                     WHERE h3_cell = ANY (h3_polygon_to_cells(:query_geom, 6)))
  AND  ST_Intersects(s.footprint, :query_geom)
  AND  s.datetime BETWEEN :t0 AND :t1
ORDER BY s.datetime DESC
LIMIT 200;
Query latency and correctness across three index designs Bounding-box-only indexing is fast but returns false positives; exact geometry alone is correct but slow; a two-stage cell-then-geometry index is both fast and correct. index designp95 latencyfalse positivesverdict bbox only one r-tree on the extent 340 ms31% fast and wrong exact geometry only GIST on the footprint 1,850 ms0% correct and slow H3 cells → geometry prune, then decide 62 ms0% both 412,000 scenes · typical query is a 20 km box over a 5-year window.

Step 3 — Keep the index rebuildable from the objects

The index is derived data and should be reconstructible in an afternoon from the archive itself, because that is what makes it safe to lose — and what makes a recovery in another region possible.

python -m archive.index rebuild \
  --source-bucket spatial-archive --prefix archive/imagery/ \
  --h3-resolution 6 --simplify-m 5 --workers 64
# 412,037 scenes indexed in 38 min (footprints read from item documents, not from rasters)

Validation & Verification 

python -m archive.index verify --sample 500
# sampled queries              500
# results matching brute force 500
# false positives                0
# false negatives                0
# p95 latency                  62 ms

Expected output is exact agreement with a brute-force scan on the sample. False negatives are the serious failure — a scene that exists and is not returned — and they usually mean the cell coverage was computed at a coarser resolution than the query uses.

Troubleshooting

Symptom Root cause Fix
Results include scenes with no overlap Only the bbox was tested Add the exact geometry test as stage two
Missing scenes near cell edges Footprint cells computed with containment, not intersection Use polygon-to-cells with intersection semantics
Index build very slow Footprints recomputed from rasters each time Read footprints from the item documents; compute once at ingest
Index far larger than expected Unsimplified footprints with thousands of vertices Simplify to the query resolution before storing
Queries fast but stale Index not rebuilt after new ingest Append on write; rebuild fully only after a schema change

Operational Execution Checklist

Keeping the Index Honest as the Archive Grows

A footprint index degrades in two ways as an archive grows: it falls behind new ingest, and its simplification tolerance becomes wrong as resolutions change. Both are detectable from the index itself, so neither needs to be discovered by a user.

Two degradations, and how each is detected Two ways a footprint index degrades, the check that detects each, and the remedy. degradation detected by remedy behind on ingest index count vs catalogue count append on write, not nightly tolerance too coarse footprint area error per collection re-simplify that collection cells at the wrong resolution cells per scene distribution rebuild the cell table stale after a re-projection sampled footprint vs file rebuild from the items All four checks read the index and the catalogue only, so they run in seconds and can be part of the nightly job.

The tolerance row is the one that appears years in, when a collection of high-resolution frames arrives into an index simplified for satellite scenes. Measuring the area error per collection rather than for the index as a whole is what makes it visible.

Frequently Asked Questions

Should the index live in a database or in files?

Both, with different roles. A database gives interactive latency and is the right home for the serving index; a set of static footprint files in the archive is the durable record from which the database can be rebuilt. Archives that keep only the database eventually discover that its schema, its host and its backups are all outside the preservation boundary.

What H3 resolution suits a footprint index?

One where a typical scene covers a handful of cells rather than one or hundreds — resolution 6 for satellite scenes and mosaic tiles at the scales most archives hold, resolution 7 or 8 for high-resolution aerial frames. The test is the cell count per scene: single digits to low tens prunes well without inflating the index.

How are time-bounded queries handled?

As an ordinary index on the datetime column, applied after the spatial prune. Spatial selectivity is usually much higher than temporal selectivity in an archive — a query covers a small area and a wide time range — so pruning spatially first and filtering on time second is the ordering that does the least work.

Should the index store footprints in geographic or projected coordinates?

In one reference system for the whole index, almost always geographic, with per-scene projected details left in the item. Mixing systems inside the index means every query has to be transformed per candidate, which defeats the purpose. The conversion cost is paid once at index time rather than repeatedly at query time.

How are very large footprints handled?

Coarse-resolution cells and a size threshold. A continental mosaic covers thousands of cells at the index resolution, which bloats the cell table for one scene. Storing such scenes at a coarser cell resolution, with a flag, keeps the table small and costs a slightly less selective prune for the few scenes that are inherently unselective anyway.

Does the index need to support time-travel queries?

If the archive publishes generations, yes — and it is cheap. Adding a validity interval to each index row lets a query ask what the archive held as of a past date, which is exactly the question an audit asks. Without it, an index describes only the present and cannot support a reproducibility claim about an earlier analysis.

Should the index include collections that are not yet published?

Include them and mark them, rather than excluding them. An index that silently omits unpublished material makes it impossible to answer questions about the archive’s total holding, and the marking costs a boolean. Filtering at query time is the right place to enforce visibility, not at index time.

What is the cost of keeping the index continuously up to date?

An append per ingested scene, which is microseconds, plus the footprint computation that already happened during validation. The expensive version is the nightly rebuild, which many archives adopt by default and rarely need — appending on write is both cheaper and fresher, and the full rebuild becomes a recovery mechanism rather than a routine.

How should the index handle scenes with no valid data at all?

Index them with an empty footprint and a flag rather than omitting them. A scene that is entirely nodata is a fact about the archive worth being able to find, and silently excluding such scenes makes the index and the catalogue disagree on the holding’s size.

Up one level: Metadata Cataloging & Discovery.

Related pages

Automating STAC Catalog Generation for Archived Imagery Extracting ISO 19115 Metadata from Legacy GeoTIFFs Keeping STAC Catalogs in Sync with Lifecycle Transitions