Converting GeoTIFFs to Cloud-Optimized GeoTIFF with rio-cogeo

One-sentence answer: rewrite with the GDAL COG driver at 512-pixel blocks, DEFLATE or ZSTD compression, predictor 2 for integers and 3 for floats, then validate the output and time a window read — because a badly laid-out file is perfectly readable and quietly ten times more expensive.

Context

Farm imagery arrives from photogrammetry software, from a colleague’s laptop, from a processing step that used gdal_translate with default options. All of it is readable, and none of it is necessarily laid out for remote access. The symptom of getting this wrong is not an error but a slow pipeline and a bill: reading a field window out of a striped mosaic pulls whole image rows, so a dashboard that shows twelve fields transfers hundreds of megabytes to display a few hundred kilobytes of pixels.

This guide implements the conversion described in cloud-optimized storage for field imagery, and its output is what windowed COG reads over HTTP with rasterio depends on.

What the conversion changes inside the file Two file layouts drawn as horizontal bars. The source file has image strips followed by a metadata directory at the very end, so a reader must seek to the end before it can read anything. The converted file starts with a header and tile index, followed by overview levels from coarsest to finest and then the full-resolution tiles. Source — striped, directory at the end full-width image strips — a window read pulls entire rows IFD reader must seek to the end first → an extra round trip on every open Converted — COG driver, 512 px blocks header ovr 1:32 ovr 1:8 ovr 1:2 full-resolution 512 × 512 tiles, addressable individually one small header request reveals every tile offset → 2–4 range requests for a field window

Prerequisites

Beyond the parent topic’s stack: rio-cogeo 5.3.* and GDAL 3.6+ with the COG driver. Source rasters must carry a CRS; a file without one is a bug to fix upstream, not a conversion input.

Step-by-step

1. Inspect the source to choose the predictor and confirm the nodata convention.

Two choices that follow from the data type alone A decision diagram driven by the band's data type. Floating-point index rasters take predictor 3 and average-resampled overviews; integer reflectance and classified rasters take predictor 2, and classified data needs nearest-neighbour overviews so it does not invent intermediate classes. Raster to convert dtype known is the band floating point? yes Predictor 3, average overviews floating-point differencing suits index rasters; averaging is right for continuous values no Predictor 2, nearest overviews horizontal differencing suits integer reflectance; nearest keeps classified rasters from inventing intermediate classes Both choices are silent when wrong: the file is valid either way, and only the compression ratio or an invented class label gives it away.

2. Convert with the COG driver.

3. Validate the output.

4. Time a window read on both files.

5. Put the validation in the publish path so a regression cannot ship.

PYTHON 
import time

import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rio_cogeo.cogeo import cog_translate, cog_validate, cog_info
from rio_cogeo.profiles import cog_profiles


def convert_to_cog(src_path: str, dst_path: str, *, blocksize: int = 512,
                   compress: str = "DEFLATE") -> dict:
    """Rewrite a raster as a cloud-optimized GeoTIFF, choosing the predictor from the dtype."""
    with rasterio.open(src_path) as src:
        assert src.crs is not None, f"{src_path} has no CRS — fix upstream, do not guess"
        dtype = src.dtypes[0]
        is_float = np.issubdtype(np.dtype(dtype), np.floating)
        if not is_float:
            assert src.nodata is not None, (
                f"{src_path} is {dtype} with no nodata value; declare one before conversion")

    profile = cog_profiles.get(compress.lower())
    profile.update({
        "blockxsize": blocksize,
        "blockysize": blocksize,
        "predictor": 3 if is_float else 2,
        "BIGTIFF": "IF_SAFER",
    })

    cog_translate(
        src_path, dst_path, profile,
        overview_level=5,
        overview_resampling="average" if is_float else "nearest",
        web_optimized=False,          # keep the native CRS; do not silently reproject to 3857
        in_memory=False,
        quiet=True,
    )
    return assert_cog(dst_path)


def assert_cog(path: str, min_overviews: int = 3) -> dict:
    """Refuse to publish anything that is not genuinely cloud-optimized."""
    valid, errors, warnings = cog_validate(path)
    assert valid, f"{path} is not a valid COG: {errors}"

    info = cog_info(path)
    overviews = len(info.IFD) - 1
    block = info.Profile.Blocksize
    assert overviews >= min_overviews, f"{path}: {overviews} overview level(s), expected ≥ {min_overviews}"
    assert block[0] >= 256 and block[1] >= 256, f"{path}: block size {block} is too small"
    if warnings:
        print(f"{path}: {len(warnings)} warning(s) — first: {warnings[0]}")
    return {"overviews": overviews, "blocksize": block, "compression": info.Profile.Compression}

Inline verification — the before-and-after read, which is the only measurement that reflects why any of this matters:

PYTHON 
def time_window_read(path: str, bounds) -> float:
    start = time.perf_counter()
    with rasterio.open(path) as src:
        win = from_bounds(*bounds, transform=src.transform)
        arr = src.read(1, window=win)
    assert arr.size > 0, f"{path}: empty window — bounds are outside the raster"
    return time.perf_counter() - start


before = time_window_read(SRC, FIELD_BOUNDS)
meta = convert_to_cog(SRC, DST)
after = time_window_read(DST, FIELD_BOUNDS)

print(f"{meta} | window read {before * 1000:.0f} ms → {after * 1000:.0f} ms")
assert after < before, "conversion did not improve the window read — was the source already tiled?"

Settings worth being deliberate about

Setting Default here Why it matters
Block size 512 px Small enough that a field window touches few blocks, large enough that per-request overhead stays negligible. Below 256 the request count dominates
Overview levels 5 Covers farm- and county-scale maps. Too few and a multi-field dashboard reads full-resolution pixels for a thumbnail
Predictor 2 integer, 3 float The wrong one silently costs a third to a half of the achievable compression ratio
Compression DEFLATE or ZSTD Lossless, 2–3.5× on reflectance and index data. Lossy codecs alter measurement values and therefore move zone boundaries
Overview resampling average for continuous, nearest for classes Averaging a zone raster produces class 2.4, which is not a class
web_optimized off On, it reprojects to Web Mercator — right for a tile server, wrong for anything a prescription derives from

Validation belongs in the publish path rather than in a notebook. A pipeline that writes with the wrong driver produces perfectly readable files, and the only symptom is a latency graph nobody is watching.

Gotchas and edge cases

  • web_optimized=True reprojects to EPSG:3857. It is the right choice for a tile server and the wrong one for analysis: resampling to Web Mercator distorts areas by a factor that grows with latitude and permanently alters the pixels. Keep the native CRS for anything a prescription will be derived from.
What each codec and predictor pairing actually saves Bars comparing compression ratios by data kind: classified zone rasters compress about twelvefold, integer reflectance about threefold, and float index rasters about 2.4 times with the correct predictor but only 1.5 times with the wrong one. Lossy JPEG reaches fifteenfold and is acceptable only for display-only browse imagery. Classified zones, uint8 ≈ 12× · DEFLATE, predictor 2 Reflectance, uint16 ≈ 3× · DEFLATE, predictor 2 Index raster, float32 ≈ 2.4× · DEFLATE, predictor 3 Index raster, float32, wrong predictor ≈ 1.5× · predictor 2 on float data Browse imagery, RGB ≈ 15× · JPEG, display only The fourth bar is the cost of predictor 2 on float32 — the file is valid, and a third of the achievable saving is simply gone.

  • Average resampling on a categorical raster invents classes. Overviews of a management-zone raster built with average produce zone 2.4. Use nearest for classified data — the convert function above branches on dtype, which covers the common case but not a categorical raster stored as float.

  • Lossy compression is not an option for measurement data. JPEG at quality 85 looks identical and perturbs every value; a zone boundary derived from it moves. Keep browse imagery as a separate derived product under its own prefix.

  • Converting in place through a temporary file needs disk headroom. A 12 GB drone mosaic needs the source, the temporary and the destination present simultaneously. in_memory=False above avoids holding it all in RAM, but the disk still has to fit two copies.

  • A validated COG can still be slow remotely. Validation checks the layout, not the reader’s configuration. If reads are slow after conversion, the cause is almost always GDAL_DISABLE_READDIR_ON_OPEN being unset — see windowed COG reads over HTTP with rasterio.

  • Run the validation in the pipeline, not in a notebook. A publish step that writes with the wrong driver produces files nothing complains about for months.

  • Conversion is a good moment to fix metadata, and a bad moment to change pixels. Setting a missing nodata value, correcting a band description or adding acquisition tags costs nothing here and saves guesswork later. Resampling, reprojecting or rescaling at the same time does not: the output then differs from the published original in ways no filename records, and two artefacts that should be identical no longer are. Keep the conversion a pure re-layout, and make any pixel change a separate, versioned derived product.

This guide is part of Cloud-Optimized Storage for Field Imagery — see there for key design, compression choices and lifecycle rules.

Parent section Cloud-Optimized Storage for Field Imagery Related Storing Multi-Temporal Index Stacks as Zarr