Writing Prediction Rasters as COGs

Derive the output profile from the input so the grids match by construction, then write tiled uint8 with a palette and mode overviews:

import rasterio
from rasterio.enums import Resampling

profile = src.profile | {"count": 1, "dtype": "uint8", "nodata": 0,
                         "compress": "deflate", "predictor": 2,
                         "tiled": True, "blockxsize": 512, "blockysize": 512}
profile.pop("photometric", None)

with rasterio.open("class.tif", "w", **profile) as dst:
    dst.write(classes, 1)
    dst.write_colormap(1, COLORS)
    dst.build_overviews([2, 4, 8, 16, 32], Resampling.mode)

Every one of those settings has a specific failure it prevents. This page belongs to exporting and serving model outputs in Raster Machine Learning & Model Inference.

A Class Raster Is Not a Reflectance Raster

Most COG advice is written for continuous imagery, and several of its defaults are actively wrong for categorical output. Averaging is forbidden, lossy compression is forbidden, and the dtype should be the smallest that holds the class codes.

Profile settings that differ for categorical output For reflectance a float or int16 dtype, average overviews and lossy compression are all reasonable. For a class raster the dtype should be uint8, overviews must use mode resampling, compression must be lossless, and a colour table is worth adding. Only tiling and block size are the same for both. Same format, different rules setting reflectance class raster dtype int16 or float32 uint8 overview resampling average mode — never average compression ZSTD, or lossy for preview DEFLATE + predictor 2 colour table not applicable yes, with 0 transparent Tiling and a 512 pixel block size are the two settings both cases share.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Profile handling, colormaps, overviews
rio-cogeo >=5.0 COG validation and single-pass conversion
numpy >=1.23 The class array 
pip install "rasterio>=1.3.0" "rio-cogeo>=5.0" "numpy>=1.23"

Complete Working Example

import numpy as np
import rasterio
from rasterio.enums import Resampling

CLASSES = {1: "forest", 2: "cropland", 3: "built", 4: "water"}
COLORS = {
    0: (0, 0, 0, 0),
    1: (38, 115, 0, 255),
    2: (168, 204, 84, 255),
    3: (196, 40, 27, 255),
    4: (28, 92, 168, 255),
}


def write_prediction_cog(classes: np.ndarray, reference_path: str, out_path: str,
                         *, model_name: str, source_ids: list[str]) -> None:
    """Write a categorical prediction aligned to a reference raster's grid."""
    if classes.dtype != np.uint8:
        raise TypeError(f"class array must be uint8, got {classes.dtype}")
    unexpected = set(np.unique(classes).tolist()) - ({0} | set(CLASSES))
    if unexpected:
        raise ValueError(f"unexpected class codes: {sorted(unexpected)}")

    with rasterio.open(reference_path) as ref:
        profile = ref.profile | {
            "count": 1, "dtype": "uint8", "nodata": 0, "driver": "GTiff",
            "compress": "deflate", "predictor": 2, "zlevel": 9,
            "tiled": True, "blockxsize": 512, "blockysize": 512,
            "BIGTIFF": "IF_SAFER",
        }
        profile.pop("photometric", None)     # inherited RGB hint breaks 1-band output

    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(classes, 1)
        dst.write_colormap(1, COLORS)
        dst.set_band_description(1, "predicted_class")
        dst.update_tags(
            model_name=model_name,
            classes=";".join(f"{k}={v}" for k, v in CLASSES.items()),
            derived_from=",".join(source_ids),
        )
        dst.build_overviews([2, 4, 8, 16, 32], Resampling.mode)
        dst.update_tags(ns="rio_overview", resampling="mode")

Validating the class codes before writing is cheap insurance. A stray value — 255 from an ignore region, or 5 from a model trained with one more class than the palette knows about — produces a file that renders as transparent holes in some viewers and as a random colour in others.

Compression and Overview Economics

What predictor 2 does to a class raster An uncompressed four-class raster of 10,980 squared pixels is 120 megabytes. Plain DEFLATE brings it to 14 megabytes, DEFLATE with horizontal differencing to 5, and the same file after speckle removal to under 2. Cleaning the output compresses it far more than any compression setting does. One tile, four ways of storing it uncompressed 120 MB deflate 14 MB deflate + predictor 2 5 MB cleaned, then the same 1.8 MB — speckle was most of the file Overviews add about a third on top of whichever bar you land on.

The last bar is the interesting one. Isolated misclassified pixels break every run in the raster, so an unfiltered prediction defeats the compressor almost as thoroughly as random noise would. Applying the majority filter and minimum mapping unit from the parent topic typically shrinks the file by a further factor of two or three — the cleanup pays for its own storage.

Overviews cost roughly a third of the base image for a factor-of-two pyramid, which is the standard geometric series. Skipping them saves that third and makes every zoomed-out view read the full resolution, so it is a false economy for anything that will be looked at. The compression trade-offs for continuous data are different and are set out in choosing COG compression: ZSTD vs DEFLATE.

Verification

import numpy as np
import rasterio
from rio_cogeo.cogeo import cog_validate

valid, errors, warnings = cog_validate("class.tif")
assert valid, errors

with rasterio.open("class.tif") as src:
    assert src.dtypes[0] == "uint8" and src.nodata == 0
    assert src.colormap(1)[0][3] == 0, "class 0 must be transparent"
    assert src.overviews(1), "no overviews built"
    full = set(np.unique(src.read(1)).tolist())
    coarse = set(np.unique(src.read(1, out_shape=(src.height // 32,
                                                  src.width // 32))).tolist())
    assert coarse <= full, "overviews contain values absent from the full band"
What averaging does to class codes A two by two block containing classes one, one, four and four averages to two point five, which rounds to class two — a class that is present in neither input pixel. Mode resampling picks one of the values actually present, so the overview stays a valid class map at every level. One 2x2 block, two resampling rules 1 1 4 4 forest and water 2 average = 2.5 -> cropland a class that was never there 1 mode = 1 -> forest a value that really occurred Repeat this over five pyramid levels and a zoomed-out map can show a composition the data never had.

The overview check is the one that catches average resampling, and it is worth running in continuous integration rather than by hand. An averaged overview introduces class codes that exist nowhere in the data, and the symptom — a map that changes composition as you zoom out — is easy to miss and embarrassing to ship. The general validation discipline is covered in validating COG structure in CI.

Common Errors

ValueError: photometric=RGB requires 3 bands

The profile was inherited from a three-band source. Pop photometric before opening the destination, as the example does.

The file is valid but enormous

predictor is missing, or the prediction still contains speckle. Check both; the second usually dominates. A quick diagnostic is to compare the compressed size against the size of the same array after a majority filter — if cleaning shrinks it by more than half, speckle rather than the compressor is the problem.

Consumers report that class 0 renders as black rather than transparent

The colour table entry for 0 has a non-zero alpha, or the viewer ignores alpha and honours only the nodata tag. Set both: alpha zero in the palette and nodata=0 in the profile, so a consumer that reads either one gets the right answer.

The colour table vanishes after processing

Several GDAL conversion paths drop palettes, particularly anything that changes dtype. Re-apply the colormap after conversion, or convert with rio cogeo create which preserves it.

Frequently Asked Questions

Q: Why DEFLATE with predictor 2 for a class raster? Horizontal differencing turns runs of identical class codes into runs of zeros, which DEFLATE compresses extremely well. A cleaned four-class raster typically drops from tens of megabytes to a few, and the decode cost is negligible.

Q: Does a prediction raster need a colour table? It makes the file self-describing, so any viewer or tile server renders the legend without configuration. It costs a few hundred bytes and saves every consumer a conversation about what class 3 means.

Q: Should I write the COG directly or convert afterwards? Write a tiled GeoTIFF during inference and convert once at the end. Maintaining strict COG byte order while blocks arrive in an arbitrary sequence forces rewrites; a single pass afterwards is far cheaper.

Q: Can several prediction bands live in one file? Yes, and for a class band plus a probability band it is usually the better arrangement: one file, one set of overviews, one thing to move. Keep the bands in a fixed documented order with descriptions set, give each its own tags describing any scaling, and remember that a colour table applies to a single band, so the probability band needs a rescale range rather than a palette. Split them into separate files only when consumers routinely want one without the other, since then the extra byte ranges are pure waste.