#Cloud-Optimized GeoTIFF for Raster Workloads

Serving raster imagery from object storage exposes a hard constraint that desktop GIS never had to face: a client should never download a 40 GB scene to render a 256-pixel map tile. Cloud-Optimized GeoTIFF (COG) solves this by reorganising an ordinary GeoTIFF so that its internal structure — tiles, overviews, and the directory that indexes them — supports partial reads. A browser, a tile server, or a Python analytics job can fetch a single tile at a single zoom level with one bounded HTTP range request, reading kilobytes instead of gigabytes. This is the raster counterpart to the columnar and streaming vector formats covered across Geospatial Storage Fundamentals & Format Comparison.

The format is deceptively simple: it is still a TIFF, readable by any GDAL-based tool with zero special handling. What makes it “cloud optimized” is the discipline of layout — square internal tiling instead of scanline strips, a full pyramid of reduced-resolution overviews, and a header placed at the front of the file so a reader learns where every tile lives before it fetches any pixels. Get that layout wrong and you have a technically valid GeoTIFF that behaves catastrophically over the network. This guide walks GIS data engineers through the structure, the GDAL COG driver and rasterio APIs, the codec trade-offs, and the boundary where COG is the right tool versus where a chunked array store like Zarr takes over.

#Prerequisites

• Python 3.10+ with rasterio>=1.3, GDAL>=3.4 (the standalone COG driver landed in GDAL 3.1), and rio-cogeo>=5.0 for validation
• A source raster: a multi-band scene (Sentinel-2, Landsat, orthophotography) or a single-band continuous grid (DEM, land-surface temperature) with a defined CRS and nodata value
• Object storage access: an S3-, GCS-, or Azure-compatible endpoint reachable through GDAL’s /vsi virtual filesystem layer for range-read testing
• Familiarity with raster fundamentals: bands, data types (uint8, int16, float32), resampling methods, and the difference between tiled and stripped TIFF organisation
• A tile-serving or analytics target: TiTiler, a WMTS endpoint, or a rasterio/xarray batch job whose access pattern you can profile

#Architectural Foundations

A GeoTIFF is a TIFF with georeferencing tags. Internally, every TIFF is described by one or more Image File Directories (IFDs) — small tables of tags that record width, height, data type, and, crucially, the byte offsets where the actual pixel data lives. In a classic stripped GeoTIFF, pixels are stored in horizontal strips and the IFD often trails at the end of the file. A reader must therefore parse the whole file to discover the layout, which over HTTP means downloading everything.

Cloud-Optimized GeoTIFF re-imposes three rules on that structure:

• Internal tiling. Pixels are stored as square tiles (commonly 256×256 or 512×512) rather than full-width strips. Each tile is an independently addressable, independently compressed block. Fetching the tile that covers a map viewport touches only that block’s byte range.
• Overviews (pyramids). The file carries a stack of decimated copies of the raster — half resolution, quarter resolution, and so on — each stored as its own tiled IFD. A zoomed-out request reads a small overview tile instead of downsampling full-resolution pixels on the fly. The mechanics of choosing levels and resampling are covered in depth in Overview Pyramids and Internal Tiling for COG.
• Front-loaded directory. All IFDs and the tile offset/byte-count arrays are written at the very start of the file. A client reads this header once — 8 to 64 KB for most scenes — learns the byte offset and length of every tile, and thereafter issues exactly one range request per tile it needs.

The diagram below shows how these three rules turn a viewport request into a single bounded read against object storage.

Because each tile is compressed independently, codec choice is a per-tile decode cost paid on every fetch — which is why raster compression tuning differs so sharply from the whole-column compression used in Parquet columnar storage for GIS. A tile server may decode thousands of tiles per second, so a codec that shaves 5% off file size but doubles decode latency is a poor trade. That tension between ratio and decode speed is the subject of ZSTD vs DEFLATE Compression for Cloud-Optimized GeoTIFF.

#Step-by-Step Workflow

#1. Assess the Source Raster

Before writing anything, characterise the input. Band count, data type, and entropy drive every downstream decision. Eight-bit RGB orthophotography behaves very differently from a 32-bit float elevation model: the former tolerates lossy WEBP for visualisation tiles, the latter demands a lossless codec and usually a horizontal predictor.

# rasterio>=1.3, GDAL>=3.4
import rasterio

with rasterio.open("scene.tif") as src:
    print("bands:", src.count)
    print("dtype:", src.dtypes[0])
    print("size:", src.width, "x", src.height)
    print("crs:", src.crs)
    print("nodata:", src.nodata)
    print("blockxsize:", src.profile.get("blockxsize"))
    # A stripped source reports blockysize of 1 — a strong signal it is not COG-ready

If blockxsize is absent or blockysize is 1, the source is stored in scanline strips and must be re-tiled. Note the data type: continuous data (int16, float32) benefits from predictor 2 or 3, whereas categorical or palette rasters do not.

#2. Choose Internal Tile Geometry

The internal tile (GDAL calls it the block) is the atomic unit of every read. A 256×256 block matches the web mapping tile convention and keeps per-request payloads small — ideal for interactive panning. A 512×512 block halves the number of tiles and the size of the offset directory, which helps analytics jobs that read large contiguous windows. The full trade-off, including how block size interacts with overview levels, is worked through in the overview and internal tiling reference.

# Rule of thumb: match block size to the dominant read pattern
block_size = 256 if serving_web_map_tiles else 512

#3. Generate Overviews

Overviews are what make zoom-out cheap. Without them, rendering a continental view means reading and downsampling every full-resolution tile in frame. Build power-of-two levels down to roughly the size of a single tile, and pick a resampling method that suits the data — average for continuous imagery, nearest for categorical classifications, mode for thematic land-cover.

# Overview factors halve resolution each step until the top level is ~1 tile
factors = [2, 4, 8, 16, 32]
resampling = "average"   # use "nearest" for categorical/palette rasters

#4. Select a Compression Codec

Every tile is compressed independently, so the codec sets both file size and per-tile decode cost. DEFLATE and LZW are the universally compatible lossless options; ZSTD is lossless and faster with better ratios where the reader supports it; WEBP offers dramatic size reductions for 8-bit visualisation tiles at the cost of being lossy (or near-lossless). Continuous data almost always wants a predictor. The selection logic below is expanded — with benchmark numbers — in the dedicated ZSTD vs DEFLATE comparison.

# Predictor 2 = horizontal differencing (integer continuous data)
# Predictor 3 = floating-point predictor (float32/float64 grids)
codec_profile = {
    "compress": "ZSTD",
    "zstd_level": 9,
    "predictor": 2,      # 3 for float rasters, 1 (off) for RGB imagery
}

#5. Validate Range-Read Behaviour

A file is only cloud optimized if it behaves that way over the network. Use rio cogeo validate to confirm the header layout, then profile actual range requests against object storage. The production section below wraps writing and validation into one reusable function.

# rio-cogeo>=5.0
from rio_cogeo.cogeo import cog_validate

is_valid, errors, warnings = cog_validate("scene_cog.tif")
print("valid COG:", is_valid, errors, warnings)

#Production-Ready Implementation

The function below converts an arbitrary source raster into a validated Cloud-Optimized GeoTIFF using the GDAL COG driver through rasterio. It picks a predictor from the data type, builds overviews, streams tiles rather than loading the whole array, and fails loudly if validation does not pass.

# rasterio>=1.3, rio-cogeo>=5.0, GDAL>=3.4
from __future__ import annotations

import logging
from pathlib import Path

import rasterio
from rasterio.enums import Resampling
from rasterio.shutil import copy as rio_copy
from rio_cogeo.cogeo import cog_validate

logger = logging.getLogger(__name__)

_FLOAT_DTYPES = {"float32", "float64"}
_CONTINUOUS_INT = {"int16", "uint16", "int32", "uint32"}


def build_cog(
    src_path: str | Path,
    dst_path: str | Path,
    *,
    compress: str = "ZSTD",
    zstd_level: int = 9,
    block_size: int = 512,
    overview_resampling: str = "average",
) -> None:
    """Convert a raster to a validated Cloud-Optimized GeoTIFF.

    Args:
        src_path: Source raster readable by GDAL/rasterio.
        dst_path: Destination COG path (.tif). Overwritten if present.
        compress: Tile codec — ZSTD, DEFLATE, LZW, or WEBP.
        zstd_level: ZSTD encoder level (ignored for non-ZSTD codecs).
        block_size: Internal tile edge in pixels (256 or 512 recommended).
        overview_resampling: 'average' for continuous, 'nearest' for categorical.

    Raises:
        ValueError: If the produced file fails COG validation.
    """
    src_path, dst_path = Path(src_path), Path(dst_path)

    with rasterio.open(src_path) as src:
        dtype = src.dtypes[0]
        if dtype in _FLOAT_DTYPES:
            predictor = 3           # floating-point predictor
        elif dtype in _CONTINUOUS_INT:
            predictor = 2           # horizontal differencing
        else:
            predictor = 1           # off — palette/8-bit RGB

        profile = src.profile.copy()
        profile.update(
            driver="COG",
            compress=compress,
            blocksize=block_size,
            overview_resampling=overview_resampling,
            BIGTIFF="IF_SAFER",
            NUM_THREADS="ALL_CPUS",
        )
        if compress.upper() == "ZSTD":
            profile["level"] = zstd_level
            profile["predictor"] = predictor
        elif compress.upper() in {"DEFLATE", "LZW"}:
            profile["predictor"] = predictor

        # The COG driver builds overviews and front-loads the directory itself.
        try:
            rio_copy(src, dst_path, **profile)
        except rasterio.errors.RasterioError as exc:
            logger.error("COG write failed for %s: %s", src_path, exc)
            raise

    is_valid, errors, warnings = cog_validate(str(dst_path))
    for w in warnings:
        logger.warning("COG warning (%s): %s", dst_path.name, w)
    if not is_valid:
        raise ValueError(f"{dst_path} is not a valid COG: {errors}")

    logger.info(
        "Wrote valid COG %s (codec=%s, block=%d, predictor=%d)",
        dst_path, compress, block_size, predictor,
    )


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    build_cog("scene.tif", "scene_cog.tif", compress="ZSTD", block_size=512)

#Codec and Layout Reference

Representative figures for a 3-band, 8-bit orthophoto scene and a single-band float32 DEM, written with the COG driver on a 4-core cloud VM. Ratios and decode speeds vary with imagery entropy and hardware — benchmark on your own data.

*WEBP supports a near-lossless mode; true lossless WEBP is limited to 8-bit data and gives lower ratios than the figures above.

For the ratio-versus-decode-speed trade-off between the lossless codecs in detail, see ZSTD vs DEFLATE for COG, which also covers how these choices echo the vector-versus-raster ZSTD level guidance from the compression topic area.

#When COG Beats Zarr and NetCDF

COG and chunked array stores like Zarr solve overlapping problems, and teams routinely over-reach with one where the other fits. The dividing line is dimensionality and access pattern:

• Reach for COG when the data is 2D or few-band imagery served as map tiles or read as spatial windows: satellite scenes, orthophotography, DEMs, single-date classifications. One self-describing file carries the pyramid, any GDAL client reads it with no runtime, and every browser mapping stack and tile server speaks it natively.
• Reach for Zarr or NetCDF when the data is a high-dimensional cube — stacked time steps, atmospheric levels, or dozens of spectral bands — where chunking must span more than the two spatial axes. A COG can hold many bands, but it has no first-class time or depth axis, so temporal analytics degrade into managing thousands of separate files.

For a mosaic of thousands of scenes, a common production pattern keeps each scene as a COG and layers a STAC catalogue plus a dynamic tiler on top, rather than rechunking everything into one array store. The single-file portability that makes COG awkward for cubes is exactly what makes it ideal for tile serving, much as GeoParquet and FlatGeobuf trade off analytical scanning against streaming range reads on the vector side.

#Failure Modes and Gotchas

#Frequently Asked Questions

When does Cloud-Optimized GeoTIFF beat Zarr or NetCDF for raster serving?

Cloud-Optimized GeoTIFF wins for 2D or few-band imagery served as map tiles, because it packs a georeferenced pyramid and internal tiles into one file that any GDAL client, tile server, or browser mapping stack reads directly over HTTP range requests. Zarr and NetCDF win for high-dimensional array cubes — time, depth, or spectral axes — where chunking across more than three dimensions matters more than single-file portability.

Do I have to use the GDAL COG driver, or can I write a valid COG with GTiff?

Both produce valid files, but the COG driver is safer. It enforces tiling, builds overviews, and front-loads the IFD and tile-offset directory so the header sits in the first few kilobytes. Writing with the GTiff driver requires you to set tiled, blockxsize, and blockysize manually, build overviews explicitly, and run a rewrite pass to move the directory to the front. The COG driver folds all of that into one call.

How large is the header a client must read before the first tile?

For a well-formed Cloud-Optimized GeoTIFF the leading header — image file directory, geokeys, and the tile offset and byte-count arrays — is typically 8 to 64 KB. GDAL fetches it in one or two range requests, caches it, and thereafter issues exactly one bounded range request per tile. A file whose directory sits at the end forces a full-file scan and defeats the format.

#Related

• Overview Pyramids and Internal Tiling for COG
• ZSTD vs DEFLATE Compression for Cloud-Optimized GeoTIFF
• Comparing GeoParquet vs FlatGeobuf Performance
• Understanding Parquet Columnar Storage for GIS
• ZSTD Compression Levels for Geospatial Data

← Back to Geospatial Storage Fundamentals & Format Comparison