Spatial Data Formats for ML Annotation: COG, GeoParquet, and STAC

A team ingests a few hundred gigabytes of aerial imagery as plain striped GeoTIFFs and stores its vector labels as Shapefiles on a shared bucket. It works for the first sprint. Then the dataset grows to tens of thousands of scenes, and three failures land in the same week. The annotation tool takes twelve seconds to open a single 2 GB scene because every tile read pulls the whole file across the network. A batch of newly exported labels silently loses attribute values because Shapefile .dbf columns truncate field names at ten characters and collapse annotation_confidence and annotation_class into the same key. And when an auditor asks which sensor and acquisition date produced a given label, nobody can answer, because the provenance lived only in a spreadsheet that stopped being updated months ago.

None of these are annotation-quality problems. They are format problems. This guide shows how three modern spatial formats — Cloud-Optimized GeoTIFF for imagery, GeoParquet for vector labels, and SpatioTemporal Asset Catalog (STAC) for the index — replace that fragile stack, and how to produce and validate each one with runnable Python. It sits inside the broader Geospatial Annotation Fundamentals & Architecture foundation, which frames why explicit data contracts matter at every stage of a spatial ML pipeline.

# Prerequisites & Format Contracts

Before converting anything, pin an environment so that COG structure, GeoParquet metadata, and STAC validation behave identically on every machine and in CI. The versions below are the ones the code in this guide is written against.

pip install \
  rasterio==1.3.10 \
  rio-cogeo==5.3.0 \
  geopandas==0.14.4 \
  pyarrow==16.1.0 \
  pystac==1.10.1 \
  shapely==2.0.6

System dependencies

• GDAL 3.6+ with the bindings that match your rasterio wheel.
• A working PROJ data directory — verify with python -c "import pyproj; print(pyproj.datadir.get_data_dir())".
• Enough scratch disk for one uncompressed copy of the largest scene you convert.

Format contracts to agree on up front

Every asset in the pipeline should declare the same projection. The first time a batch mixes projections is the first time labels land in the wrong place, so document a canonical coordinate reference system — usually a local UTM zone such as EPSG:32633 — and require imagery, labels, and catalog to agree on it. Decide a tile block size (512 pixels is a safe default), an overview resampling method, and a STAC identifier scheme that ties each label file back to the scene it annotates.

# Core Format Workflow

The three formats are not competitors; they occupy three roles. COG makes imagery streamable, GeoParquet makes labels queryable, and STAC makes the whole set discoverable. Build them in that order.

# Produce a Cloud-Optimized GeoTIFF for Imagery

A plain GeoTIFF stores pixels in horizontal strips with the metadata directory anywhere in the file. To read a 256×256 window a client may have to walk most of the file. A COG fixes this by writing internally tiled pixels, precomputed overviews (downsampled zoom levels), and a leading image-file-directory, so a client can read the file header, learn the byte ranges it needs, and fetch only those with HTTP range requests. That is the single change that turns a twelve-second scene open into a sub-second one.

from pathlib import Path

from rio_cogeo.cogeo import cog_translate
from rio_cogeo.profiles import cog_profiles


def build_cog(src_path: str, dst_path: str, blocksize: int = 512) -> Path:
    """Translate a source raster into a Cloud-Optimized GeoTIFF.

    Uses DEFLATE compression, internal tiling at ``blocksize``, and
    power-of-two overviews so annotation clients can stream tiles.
    """
    profile = cog_profiles.get("deflate")
    profile.update(blockxsize=blocksize, blockysize=blocksize)

    cog_translate(
        src_path,
        dst_path,
        profile,
        overview_level=5,
        overview_resampling="average",
        web_optimized=False,
        in_memory=False,
    )
    return Path(dst_path)


if __name__ == "__main__":
    build_cog("data/raw/scene_0421.tif", "data/cog/scene_0421.tif")

overview_resampling="average" suits continuous imagery; for a categorical raster such as a rasterized label mask use "nearest" so class values are never blended. Keep web_optimized=False unless you are serving a web map, because web optimization reprojects to Web Mercator and you generally want annotation imagery to stay in its native projection.

# Validate the COG Structure

Producing a COG is not the same as producing a valid one — a wrong block size, missing overviews, or a trailing IFD all silently degrade to full-file reads. Validate structure explicitly rather than trusting the writer.

rio cogeo validate data/cog/scene_0421.tif

Wrap the same check in Python so it can run as a gate over a whole directory:

from rio_cogeo.cogeo import cog_validate


def assert_valid_cog(path: str) -> None:
    """Raise if ``path`` is not a spec-compliant COG."""
    is_valid, errors, warnings = cog_validate(path)
    if warnings:
        print(f"{path}: warnings -> {warnings}")
    if not is_valid:
        raise ValueError(f"{path} is not a valid COG: {errors}")
    print(f"{path}: valid COG")


if __name__ == "__main__":
    assert_valid_cog("data/cog/scene_0421.tif")

cog_validate returns errors for hard spec violations (no internal tiling, overviews out of order) and warnings for softer issues (block size larger than the image, too few overview levels). Treat errors as build-breaking and log warnings for review.

# Store Vector Labels as GeoParquet

Shapefiles are where annotation attributes go to die: ten-character field names, a 2 GB size ceiling, no native typing beyond a handful of .dbf types, and a CRS held in a fragile external .prj. GeoParquet replaces all of that with a columnar, compressed, strongly typed format that embeds the projection as PROJJSON inside the file. Columns are read independently, so filtering a million-feature label set by class touches only the class column.

import geopandas as gpd
from shapely.geometry import Polygon


def write_label_geoparquet(dst_path: str, target_epsg: int = 32633) -> None:
    """Write example annotation polygons to CRS-aware GeoParquet."""
    features: list[dict[str, object]] = [
        {
            "annotation_id": "a-0001",
            "annotation_class": "building",
            "annotation_confidence": 0.94,
            "geometry": Polygon([(0, 0), (0, 10), (10, 10), (10, 0)]),
        },
        {
            "annotation_id": "a-0002",
            "annotation_class": "solar_panel",
            "annotation_confidence": 0.81,
            "geometry": Polygon([(20, 20), (20, 26), (28, 26), (28, 20)]),
        },
    ]
    gdf = gpd.GeoDataFrame(features, geometry="geometry", crs=f"EPSG:{target_epsg}")
    # Full attribute names survive; no .dbf truncation.
    gdf.to_parquet(dst_path, index=False)


def read_and_check(dst_path: str, expected_epsg: int = 32633) -> None:
    gdf = gpd.read_parquet(dst_path)
    actual = gdf.crs.to_epsg()
    assert actual == expected_epsg, f"CRS drift: expected {expected_epsg}, got {actual}"
    print(f"Loaded {len(gdf)} labels, CRS EPSG:{actual}, columns={list(gdf.columns)}")


if __name__ == "__main__":
    write_label_geoparquet("data/labels/scene_0421.parquet")
    read_and_check("data/labels/scene_0421.parquet")

The full attribute names annotation_class and annotation_confidence round-trip intact — the exact failure the Shapefile stack could not survive. If your labels start life as Shapefiles, the dedicated walkthrough on converting Shapefiles to GeoParquet for annotation covers CRS preservation and field-name recovery across the roundtrip.

# Catalog Assets with STAC

With imagery in COGs and labels in GeoParquet, you still need to know which label file annotates which scene, in which projection, captured when, by which sensor. STAC answers that. A STAC Item is a GeoJSON Feature describing one spatiotemporal asset group; a STAC Collection groups related Items and holds shared metadata. Point the Item’s assets at the COG and the GeoParquet, and the catalog becomes the single source of truth for provenance.

from datetime import datetime, timezone

import pystac
from shapely.geometry import box, mapping


def build_stac_item(
    item_id: str,
    cog_href: str,
    label_href: str,
    bbox: tuple[float, float, float, float],
    acquired: datetime,
    epsg: int = 32633,
) -> pystac.Item:
    """Create a STAC Item linking a COG scene to its GeoParquet labels."""
    geometry = mapping(box(*bbox))
    item = pystac.Item(
        id=item_id,
        geometry=geometry,
        bbox=list(bbox),
        datetime=acquired,
        properties={"proj:epsg": epsg},
    )
    item.add_asset(
        "imagery",
        pystac.Asset(href=cog_href, media_type=pystac.MediaType.COG, roles=["data"]),
    )
    item.add_asset(
        "labels",
        pystac.Asset(
            href=label_href,
            media_type="application/x-parquet",
            roles=["labels"],
        ),
    )
    return item


def build_collection(items: list[pystac.Item]) -> pystac.Collection:
    spatial = pystac.SpatialExtent([[-180.0, -90.0, 180.0, 90.0]])
    temporal = pystac.TemporalExtent([[datetime(2024, 1, 1, tzinfo=timezone.utc), None]])
    collection = pystac.Collection(
        id="annotation-scenes",
        description="COG imagery with GeoParquet annotation labels.",
        extent=pystac.Extent(spatial, temporal),
    )
    for item in items:
        collection.add_item(item)
    return collection


if __name__ == "__main__":
    item = build_stac_item(
        item_id="scene_0421",
        cog_href="data/cog/scene_0421.tif",
        label_href="data/labels/scene_0421.parquet",
        bbox=(300000.0, 5800000.0, 301000.0, 5801000.0),
        acquired=datetime(2024, 6, 1, 10, 30, tzinfo=timezone.utc),
    )
    collection = build_collection([item])
    collection.normalize_hrefs("data/catalog")
    collection.save(catalog_type=pystac.CatalogType.SELF_CONTAINED)

The proj:epsg property comes from the STAC Projection extension and records the asset’s native CRS, so a consumer never has to open the COG to learn how to reproject the labels. normalize_hrefs followed by save writes a self-contained catalog whose links resolve relative to the catalog root.

# Format Role & Validation Reference

# Edge Cases & Gotchas

# COG overviews and block size mismatches

The most common invalid COG has correct pixels but wrong structure. If the block size exceeds the image dimensions, the file has a single internal tile and behaves like a striped GeoTIFF under range reads. If overviews are missing, a zoomed-out annotation view forces a full-resolution read of the entire scene. Pick a block size (512) smaller than your smallest scene, request enough overview levels to reach a browsable thumbnail, and let cog_validate warnings flag both conditions before the file ships.

# GeoParquet CRS metadata and GeoArrow encoding

GeoParquet keeps CRS as PROJJSON in file-level geo metadata, not per row. If you build a GeoDataFrame without setting crs, the written file carries a null projection and every downstream reader has to guess — the same silent failure a missing .prj causes. Always set crs before writing. Be aware that GeoParquet 1.1 introduced GeoArrow-encoded geometry columns alongside the original well-known-binary encoding; pin your geopandas and pyarrow versions together so writer and reader agree on which encoding they emit and expect.

# STAC required fields and null datetime

A STAC Item that omits any of type, stac_version, id, geometry, bbox, properties.datetime, assets, or links will fail validation and be rejected by most catalog servers. The trap is datetime: for a scene captured over a time window rather than an instant, datetime may be null, but only if you supply start_datetime and end_datetime in properties instead. Validate items rather than assuming the constructor enforced this.

# Mixing formats with inconsistent projections

The formats compose cleanly only when their projections agree. A COG in one UTM zone with labels in a neighboring zone will render but place every annotation tens of thousands of meters off. Assert that each Item’s proj:epsg matches the CRS recorded in its GeoParquet asset, and reproject explicitly where they differ rather than trusting a downstream tool to reconcile them silently.

# Integration & Automation Hooks

# Serving COGs to Label Studio

When you connect Label Studio integrated with geospatial workflows to a bucket of COGs, the tiling backend issues range requests against each scene, so a valid COG is the difference between an annotator waiting on a full download and one panning smoothly. Gate the ingest step with assert_valid_cog so no striped GeoTIFF ever reaches the labeling queue, and expose the STAC Collection so tasks can be created directly from catalog Items rather than a hand-maintained file list.

# Versioning formatted assets with DVC

Once imagery is in COGs and labels in GeoParquet, track them with DVC for geospatial training data. Content-addressed storage plays well with these formats: a COG’s bytes are stable across re-exports when you fix the compression and overview settings, and GeoParquet’s deterministic column layout gives reproducible hashes. Commit the STAC catalog alongside the data so every dataset version resolves to the exact imagery and labels its Items reference.

# Validation & Testing

Structure validation and geometry validation are separate concerns — a file can be a perfect COG full of nonsense, or valid GeoParquet holding self-intersecting polygons. Run both.

# Catalog and item validation

import pystac


def validate_catalog(catalog_path: str) -> None:
    """Validate every Item and Collection under a STAC catalog root."""
    catalog = pystac.Catalog.from_file(catalog_path)
    n = catalog.validate_all()
    print(f"Validated {n} STAC objects successfully")


if __name__ == "__main__":
    validate_catalog("data/catalog/collection.json")

# Geometry and CRS integrity check

import geopandas as gpd


def check_label_integrity(path: str, expected_epsg: int = 32633) -> None:
    """Assert labels are valid, non-empty, and in the expected CRS."""
    gdf = gpd.read_parquet(path)

    assert gdf.crs is not None, f"{path}: missing CRS metadata"
    assert gdf.crs.to_epsg() == expected_epsg, (
        f"{path}: CRS mismatch, expected EPSG:{expected_epsg}"
    )
    assert gdf.geometry.is_valid.all(), f"{path}: contains invalid geometries"
    assert not gdf.geometry.is_empty.any(), f"{path}: contains empty geometries"
    print(f"{path}: {len(gdf)} labels valid in EPSG:{gdf.crs.to_epsg()}")


if __name__ == "__main__":
    check_label_integrity("data/labels/scene_0421.parquet")

Wire the COG validator, the catalog validator, and the geometry check into one pre-training gate. If any asset fails, the batch never enters the training queue — the cheapest place to catch a format defect is before it becomes a corrupted label a model learns from. Reprojection edge cases, such as reconciling a Shapefile whose .prj disagrees with its coordinates, are covered in depth by the coordinate reference systems guide.

# Frequently Asked Questions

# Related

• Converting Shapefiles to GeoParquet for Annotation — a step-by-step migration off legacy Shapefiles, preserving CRS and recovering truncated field names
• Coordinate Reference Systems in Annotation Pipelines — how to normalize the projections these formats embed so imagery and labels stay aligned
• Vector vs Raster Annotation Workflows — when labels belong in GeoParquet vectors versus rasterized COG masks
• Implementing DVC for Geospatial Training Data — version-controlling COGs, GeoParquet, and the STAC catalog together

This guide is part of the broader Geospatial Annotation Fundamentals & Architecture foundation that underpins every production ML pipeline working with spatial data.