Cloud-Optimized Storage for Field Imagery

Storage layout is the difference between a platform whose imagery costs a few dollars a month to query and one that reads a hundred megabytes to answer a question about half a hectare. The output of this topic is a concrete layout: internally tiled GeoTIFFs with overviews for per-date imagery, Zarr stores chunked on time for index stacks, object keys designed around both the queries and the expiry rules, and validation that proves the writer actually produced what was asked for. It sits under farm data platform engineering, and it is the storage half of the metadata table described in PostGIS schema design for farm data.

Prerequisites

• Python 3.11+, rasterio 1.3., rio-cogeo 5.3., xarray 2024., zarr 2.18., numcodecs 0.13.*
• GDAL 3.6+ with the COG driver, and the HTTP environment settings from the section overview
• An object store with range-request support and lifecycle rules
• Imagery to store: orthomosaics from orthomosaic stitching workflows, windows from satellite imagery APIs, or index rasters from band math and raster algebra
• All rasters carrying an explicit CRS — a file without one is not storable, it is a bug, as fixing rasterio CRS mismatch errors explains

1. Concept: The Layout Is the API

A remote reader can only be fast if the file tells it where to look and lets it fetch a small piece. Three properties provide that.

Internal tiling. A striped GeoTIFF stores pixels row by row across the full image width, so reading a 500 × 500 m window pulls every row it touches from edge to edge. A tiled GeoTIFF stores 512 × 512 pixel blocks, and a window read fetches only the blocks it overlaps — for a field inside a large mosaic that is often a hundredfold less data.

Overviews. Reduced-resolution copies embedded in the same file let a map at farm scale read a 1:64 pyramid level instead of decimating full-resolution pixels. Without them, a dashboard showing twelve fields reads twelve full-resolution mosaics.

Header at the front. The tile offsets must be readable in one small request near the start of the object. GDAL then knows exactly which byte ranges to ask for. A file whose directory sits at the end costs an extra round trip per open, which at four hundred fields times thirty dates is a meaningful fraction of the runtime.

2. Step-by-Step Implementation

Step 1 — Write a genuine cloud-optimized GeoTIFF

import rasterio
from rasterio.shutil import copy as rio_copy

COG_PROFILE = {
    "driver": "COG",
    "blocksize": 512,
    "compress": "DEFLATE",
    "predictor": 2,          # 2 for integer bands; 3 for float32 — see step 2
    "overview_resampling": "average",
    "overview_levels": 5,
    "num_threads": "ALL_CPUS",
    "BIGTIFF": "IF_SAFER",
}


def write_cog(src_path: str, dst_path: str, *, float_data: bool = False) -> None:
    """Rewrite a raster as a cloud-optimized GeoTIFF with overviews."""
    profile = dict(COG_PROFILE, predictor=3 if float_data else 2)
    with rasterio.open(src_path) as src:
        assert src.crs is not None, f"{src_path} has no CRS — refusing to publish"
        assert src.nodata is not None or float_data, (
            f"{src_path} has no nodata value; integer rasters must declare one")
        rio_copy(src, dst_path, **profile)

The COG driver builds overviews and orders the file correctly on its own; the older recipe of writing a tiled TIFF and then calling gdaladdo still works but produces a file whose directory placement depends on GDAL version. Prefer the driver.

Step 2 — Choose compression by what the pixels mean

Lossy compression on measurement data is the mistake worth naming twice. A JPEG-compressed NDVI raster looks identical on screen and has had every value perturbed by up to a few hundredths — enough to move a zone boundary, and therefore enough to change an application rate. Keep browse imagery as a separate, clearly named derived product if you need one.

Step 3 — Design keys around queries and expiry

An object key is not a filename; it is the only index the object store itself provides, and lifecycle rules act on prefixes. A layout that works:

published/s2/{tile}/{acquired_date}/{band}.tif
published/flights/{field_id}/{flight_ts}/ortho.tif
derived/ndvi/{code_version}/{field_id}/{acquired_date}.tif
derived/zones/{code_version}/{field_season_id}.tif
staging/{run_id}/{whatever}
browse/{field_id}/{acquired_date}.jpg

Three properties fall out of it. Lifecycle rules are expressible: expire staging/ after seven days, transition published/s2/ to infrequent access after ninety, delete browse/ after a season. Derived products are versioned by algorithm, so recomputing writes beside the old rather than over it. And the field-scoped prefixes mean “everything for this field” is one prefix listing rather than a database round trip — useful precisely when the database is the thing that is down.

Avoid dates at the front of a key (2026/06/12/field_a/...). It reads naturally and it means every object written in a week shares a prefix, which is exactly the hot-spotting pattern that trips per-prefix request limits during a backfill.

Step 4 — Zarr for the time dimension

A season of index rasters for one field is a small three-dimensional array. Written as thirty COGs it answers “show me 12 June” instantly and “compute the seasonal integral per pixel” slowly. Written as one Zarr store chunked on time it does the reverse.

import xarray as xr
import numpy as np
from numcodecs import Blosc

def write_index_stack(da: xr.DataArray, store: str) -> None:
    """Persist a (time, y, x) index stack chunked for per-pixel time series."""
    assert da.dims == ("time", "y", "x"), f"unexpected dims {da.dims}"
    assert da.dtype == np.float32, f"store indices as float32, not {da.dtype}"

    chunked = da.chunk({"time": da.sizes["time"], "y": 256, "x": 256})
    encoding = {da.name: {
        "compressor": Blosc(cname="zstd", clevel=5, shuffle=Blosc.SHUFFLE),
        "chunks": (da.sizes["time"], 256, 256),
        "_FillValue": np.float32("nan"),
    }}
    chunked.to_dataset().to_zarr(store, mode="w", encoding=encoding, consolidated=True)

The whole time axis in one chunk is the point: a per-pixel seasonal statistic touches one chunk per 256 × 256 spatial block instead of one object per date. consolidated=True writes a single metadata object so opening the store is one request rather than one per array — on object storage that alone is usually a tenfold difference in open time.

Step 5 — Validate the layout, do not assume it

from rio_cogeo.cogeo import cog_validate, cog_info

def assert_cog(path: str, min_overviews: int = 3) -> None:
    valid, errors, warnings = cog_validate(path)
    assert valid, f"{path} is not a valid COG: {errors}"
    info = cog_info(path)
    ovr = len(info.IFD) - 1
    assert ovr >= min_overviews, f"{path} has {ovr} overview levels, expected ≥ {min_overviews}"
    assert info.Profile.Blocksize[0] >= 256, f"{path} block size {info.Profile.Blocksize} is too small"
    if warnings:
        print(f"{path}: {len(warnings)} COG warning(s) — {warnings[0]}")

Run this in the publication step, not in a notebook. A pipeline that writes with the wrong driver produces perfectly readable files that are quietly slow, and nothing downstream will ever complain — the symptom is a bill and a latency graph, not an error.

3. Key Parameters and Tuning

4. Edge Cases and Failure Modes

A COG that is not one. Files pass through resampling steps, gdal_translate invocations and other tools that silently rewrite them as striped. Validate on publication, and validate again if a file has been touched by anything outside your pipeline.

Nodata versus NaN. Integer rasters need an explicit nodata value; float rasters can use NaN, but only if every consumer honours it. Mixing the two — a uint16 reflectance raster with nodata 0 feeding a float32 index raster with NaN — is the propagation problem in fixing nodata and NaN propagation in band math. Decide per data type and record the choice in the metadata table.

Eventual consistency on listing. An object written a second ago may not appear in a prefix listing. Never use a listing to decide whether work is done — that is what the manifest is for.

Requester-pays and cross-region reads. Reading a public archive from a different region works and costs egress on every range request. At backfill scale that is real money. Co-locate compute with the archive, or copy the windows you need once and read locally thereafter.

Small-object explosion. One COG per field per band per date is 400 × 4 × 30 = 48,000 objects a season. That is fine for storage and awkward for anything that lists. It is the reason index stacks belong in Zarr and the reason browse imagery gets its own prefix with an aggressive expiry.

Overwriting published objects. A rerun that writes the same key with different bytes silently changes what an old prescription was based on. Published keys are write-once; if the bytes must change, the identity was wrong.

5. Verification and Output Validation

import time
import rasterio
from rasterio.windows import from_bounds

def assert_window_read_is_cheap(url: str, bounds, max_seconds: float = 1.5) -> None:
    """A field window on a well-formed COG should be fast even over HTTP."""
    start = time.perf_counter()
    with rasterio.open(url) as src:
        window = from_bounds(*bounds, transform=src.transform)
        arr = src.read(1, window=window)
    elapsed = time.perf_counter() - start
    assert arr.size > 0, "window read returned no pixels — bounds are outside the raster"
    assert elapsed < max_seconds, (
        f"window read took {elapsed:.2f}s — file is probably striped, or "
        "GDAL_DISABLE_READDIR_ON_OPEN is unset")

This is the check that catches the failure nothing else does. A striped file, a missing GDAL setting and a cross-region read all produce correct pixels and a slow pipeline; asserting the latency of a representative read turns an invisible regression into a failing test. Run it against one published object per kind on every deploy, and give it a deliberately striped file in the test suite to prove it fires.

6. Integration with the Broader Pipeline

Everything that produces pixels writes here: windows from satellite imagery APIs and STAC catalogues, orthomosaics from batch and async orthomosaic processing, index rasters from calculating NDVI and NDRE with rasterio, and classified zones from threshold mapping for crop health. Everything that consumes pixels reads from here, with the raster_artifact row in PostGIS schema design for farm data as the index, and the writes are sequenced by orchestrating seasonal pipelines with Airflow.

Format choice for the vector half of the platform is a separate decision, covered in GeoPackage vs GeoParquet vs Shapefile for farm data.

Frequently Asked Questions

How much does a season of imagery actually cost to store? For 400 fields averaging 45 ha, a season of four Sentinel-2 bands clipped to field windows across 30 dates is roughly 20–40 GB after compression — a few dollars a month. Full-resolution drone orthomosaics are the expensive item: a single 45 ha five-band mosaic at 5 cm is 8–15 GB before compression, so flight archives are where lifecycle policy earns its keep.

Should I keep raw drone imagery or only the mosaic? Keep the mosaic and the capture log always. Keep raw frames only while you might reprocess — a season is a reasonable default — because raw frames are an order of magnitude larger than the mosaic they produce, and reprocessing them a year later usually means a different photogrammetry version anyway.

Is Zarr worth it for a single farm? Below about ten fields, no — thirty COGs and a loop is simpler and fast enough. The crossover comes when per-pixel seasonal statistics across many fields become routine, at which point the difference between one contiguous chunk and thirty remote opens is the difference between seconds and minutes per field.

This topic is part of Farm Data Platform Engineering: APIs, Storage & Orchestration — see there for the manifest and metadata patterns that index these objects.

Related

• Converting GeoTIFFs to Cloud-Optimized GeoTIFF with rio-cogeo — a runnable conversion with validation and before-and-after timing
• Storing Multi-Temporal Index Stacks as Zarr — chunking for per-pixel time series and appending a new date safely
• Satellite Imagery APIs & STAC Catalogues — the windowed reads this layout is designed to make cheap
• PostGIS Schema Design for Farm Data — the metadata rows that point at these objects
• GeoPackage vs GeoParquet vs Shapefile for Farm Data — the same decision for the vector half of the platform