Exporting and Serving Model Outputs

The array that comes out of a model is not yet a product. It has speckle a human would never draw, it carries no legend, no nodata convention anyone else knows about, no record of which model produced it, and no internal layout that lets a tile server read one neighbourhood without reading the whole file. Turning it into something other systems can consume is the last stage of Raster Machine Learning & Model Inference, and it is the stage that determines whether the work gets used.

The target format is settled: a Cloud-Optimized GeoTIFF with a class band, a probability band, internal overviews built with mode resampling, and provenance in the tags. Everything in this topic is in service of producing that reliably.

---

Prerequisites

pip install "rasterio>=1.3.0" "rio-cogeo>=5.0" "scipy>=1.11" "shapely>=2.0" "geopandas>=0.14" "pystac>=1.10"

The COG mechanics assumed here are covered in writing and validating Cloud-Optimized GeoTIFFs, and the catalog side in querying STAC catalogs programmatically.

---

Step-by-Step Workflow

Step 1 — Clean the prediction before writing it

Raw pixel-based output is speckled: isolated pixels of one class inside a solid field of another. Some of that speckle is real — a farm building in a field is genuinely built-up — but most of it is model noise at the margin of a decision boundary, and it makes the product look unserious.

import numpy as np
from scipy import ndimage


def majority_filter(classes: np.ndarray, size: int = 3) -> np.ndarray:
    """Replace each pixel with the most common class in its neighbourhood."""
    out = np.zeros_like(classes)
    counts = np.zeros((int(classes.max()) + 1, *classes.shape), dtype="uint8")
    for cls in range(1, int(classes.max()) + 1):
        counts[cls] = ndimage.uniform_filter(
            (classes == cls).astype("float32"), size=size
        ) * size * size + 0.5
    out = counts.argmax(axis=0).astype("uint8")
    return np.where(classes == 0, 0, out)      # never invent data over nodata


def apply_minimum_mapping_unit(classes: np.ndarray, min_pixels: int = 10) -> np.ndarray:
    """Dissolve connected patches smaller than min_pixels into their surroundings."""
    out = classes.copy()
    for cls in np.unique(classes[classes > 0]):
        labelled, n = ndimage.label(classes == cls)
        if n == 0:
            continue
        sizes = ndimage.sum_labels(np.ones_like(labelled), labelled, range(1, n + 1))
        small = np.isin(labelled, np.flatnonzero(sizes < min_pixels) + 1)
        out[small] = 0                          # re-filled by the majority pass below
    return np.where(out == 0, majority_filter(classes), out)

The np.where(classes == 0, 0, out) guard is the important line. A majority filter over the whole array will happily fill nodata holes with whatever surrounds them, quietly extending the map beyond the imagery that supports it.

Step 2 — Write the class raster with a legend

A class raster without a colour table is a grey smear in every viewer that opens it. GDAL supports a palette on a uint8 band, and writing one costs three lines.

import rasterio

COLORS = {
    0: (0, 0, 0, 0),           # nodata, fully transparent
    1: (38, 115, 0, 255),      # forest
    2: (168, 204, 84, 255),    # crops
    3: (196, 40, 27, 255),     # built
    4: (28, 92, 168, 255),     # water
}


def write_class_cog(classes, profile, path: str) -> None:
    profile = profile | {
        "count": 1, "dtype": "uint8", "nodata": 0,
        "compress": "deflate", "predictor": 2,
        "tiled": True, "blockxsize": 512, "blockysize": 512,
    }
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(classes, 1)
        dst.write_colormap(1, COLORS)
        dst.set_band_description(1, "land_cover_class")
        dst.build_overviews([2, 4, 8, 16, 32], rasterio.enums.Resampling.mode)
        dst.update_tags(ns="rio_overview", resampling="mode")

Mode resampling for the overviews is not a stylistic choice. Average resampling on class codes produces pixels with values like 2.6, which round to a class that may be nowhere near either neighbour — a forest and a water pixel averaging to “built”. At low zoom the whole map becomes wrong in a way nobody notices until someone screenshots it.

Step 3 — Keep the probability, cheaply

The probability band is what lets a consumer re-threshold without re-running the model, and it is what any uncertainty discussion rests on. Storing it as float32 doubles the product size for precision nobody uses; storing the winning probability scaled to a byte costs one extra band and loses nothing that matters.

import numpy as np
import rasterio


def write_probability_band(prob: np.ndarray, profile, path: str) -> None:
    """Store the winning class probability as 1-255, reserving 0 for nodata."""
    scaled = np.clip(np.round(prob * 254) + 1, 1, 255).astype("uint8")
    scaled[~np.isfinite(prob)] = 0
    profile = profile | {"count": 1, "dtype": "uint8", "nodata": 0,
                         "compress": "deflate", "tiled": True}
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(scaled, 1)
        dst.set_band_description(1, "class_probability")
        dst.update_tags(1, scale="1/254", offset="-1/254",
                        note="value 0 is nodata; p = (v - 1) / 254")

Recording the decoding rule in the band tags is what makes the scaling survive contact with a consumer who did not read your documentation. The same convention applies to packed reflectance, as described in packing float rasters into int16 with scale and offset.

Step 4 — Vectorize when polygons are the deliverable

Many consumers want features, not pixels: a count of buildings, a field boundary layer, an area per administrative unit. rasterio.features.shapes converts a mask to polygons, and the work afterwards is all about cleaning.

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import shapes
from shapely.geometry import shape


def vectorize_class(path: str, target_class: int,
                    simplify_m: float = 10.0) -> gpd.GeoDataFrame:
    with rasterio.open(path) as src:
        arr = src.read(1)
        mask = (arr == target_class).astype("uint8")
        geoms = [
            {"geometry": shape(geom), "class_id": target_class}
            for geom, value in shapes(mask, mask=mask.astype(bool),
                                      transform=src.transform)
            if value == 1
        ]
        gdf = gpd.GeoDataFrame(geoms, crs=src.crs)

    gdf["geometry"] = gdf.geometry.simplify(simplify_m, preserve_topology=True)
    gdf["area_ha"] = gdf.geometry.area / 10_000
    return gdf[gdf.area_ha > 0.1].reset_index(drop=True)

Two details make the difference between a usable layer and a million-feature mess. Passing mask= to shapes restricts the walk to the pixels of interest instead of polygonising the background as one enormous multi-part shape. And simplifying with a tolerance tied to the pixel size — roughly one pixel — removes the staircase edges without moving any boundary further than the imagery can resolve.

Step 5 — Record what produced it

import pystac
import rasterio
from datetime import datetime, timezone


def build_stac_item(class_path: str, prob_path: str, *, item_id: str,
                    model: str, contract_sha: str,
                    source_ids: list[str]) -> pystac.Item:
    with rasterio.open(class_path) as src:
        bounds, crs = src.bounds, src.crs

    item = pystac.Item(
        id=item_id,
        geometry=None, bbox=list(bounds),
        datetime=datetime.now(timezone.utc),
        properties={
            "model:name": model,
            "model:contract_sha": contract_sha,
            "derived_from": source_ids,
            "proj:epsg": crs.to_epsg(),
        },
    )
    item.add_asset("class", pystac.Asset(href=class_path, media_type=pystac.MediaType.COG,
                                         roles=["data"], title="Predicted class"))
    item.add_asset("probability", pystac.Asset(href=prob_path, media_type=pystac.MediaType.COG,
                                               roles=["data"], title="Class probability"))
    return item

---

Serving the Product Once It Exists

A published prediction usually has two audiences with opposite needs. Analysts want the file, byte-accessible, so they can run their own statistics; everyone else wants to look at it on a map without downloading two gigabytes.

Both are served by the same COG, which is the main reason the format is worth the extra care at write time. An analyst reads windows over HTTP with the patterns in reading a COG over S3 without downloading; a viewer requests tiles from a dynamic tile server that reads the same overviews. Nothing needs to be duplicated, and there is exactly one file whose contents can be authoritative.

Two properties of the file make or break the tiling path. The internal block size must match what the tiler requests — 512 is the common choice and it is why every write in this topic sets it — and the overviews must go down far enough that a whole-country zoom level reads a few hundred kilobytes rather than the full resolution. A COG with no overviews serves its first tile perfectly and its zoomed-out view by reading everything.

For categorical products there is one more consideration: the colour table travels with the file, so a tile server that honours palettes renders the legend without configuration. That is a strong argument for writing the palette even when your own viewer applies its own styling — it makes the file self-describing for consumers you will never meet. The tiling side of this is covered in serving raster tiles with TiTiler.

---

Parameter Reference

---

Verification & Testing

Validate the COG structure rather than assuming it:

from rio_cogeo.cogeo import cog_validate

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

Then check the three properties that matter for a categorical product specifically. The set of values present must be exactly the declared classes plus zero — a stray 7 in a four-class map means the post-processing introduced a value. The overviews must contain only those same values, which is the direct test that mode resampling was used. And the nodata footprint of the prediction must match the nodata footprint of the input, because any difference means the cleanup either eroded valid data or extended the map past its evidence.

import numpy as np
import rasterio

with rasterio.open("class.tif") as src:
    full = np.unique(src.read(1))
    over = np.unique(src.read(1, out_shape=(src.height // 16, src.width // 16)))
assert set(full.tolist()) <= {0, 1, 2, 3, 4}, full
assert set(over.tolist()) <= set(full.tolist()), "overviews invented class values"

---

Troubleshooting

The colour table disappears after a gdal_translate

Palettes are dropped by several conversion paths, notably any that change the dtype or apply -expand rgb. Re-apply the colormap after conversion, or convert with rio cogeo create which preserves it.

Vectorizing produces one polygon with a million holes

The mask was not passed to shapes, so the background was polygonised as a single feature containing every hole. Pass mask= as the code above does.

The output file is enormous for a four-class map

predictor=2 is missing, or the file was written untiled so every block is stored independently. A cleaned four-class raster at 10,980² should compress to a few megabytes; hundreds of megabytes means the speckle survived and is defeating the compressor.

Overviews look like a different map at low zoom

Average or bilinear resampling was used. Rebuild with Resampling.mode and set the rio_overview tag so later tools do not undo it.

The STAC item is rejected by the catalog

Usually a missing datetime or a bbox in projected coordinates. STAC bboxes are in WGS84 degrees regardless of the asset CRS — transform the bounds before setting them, using the approach in transforming point coordinates with pyproj.

---

Frequently Asked Questions

Q: Should I ship the class raster or the polygons? Ship the raster as the authoritative product and the polygons as a convenience derivative. The raster is lossless with respect to what the model produced; every vectorization applies a simplification tolerance and a minimum area that throws information away, and different consumers want different tolerances.

Q: How do I store per-class probabilities without a huge file? Store only the winning class and its probability, scaled to a byte, which costs two bands instead of one per class. Full per-class probability stacks are worth keeping only when downstream consumers genuinely re-combine classes, and then they should be written as a separate asset.

Q: What minimum mapping unit should I apply? Whatever the product specification says, and if there is none, choose it from the sensor rather than the model: five to ten pixels at 10 m is a sensible floor for land cover, because a three-pixel patch is within the registration error of the imagery itself.

Q: Should the unfiltered prediction be published too? Keep it, publish it only if someone asks. The cleaned raster is the product; the raw one is provenance. Storing it costs little under compression and it is the only way to re-derive a different minimum mapping unit later without re-running inference.

Q: Can the same pipeline publish an update without breaking existing links? Yes, if the item id encodes the model version rather than only the area and date. Write a new item and new assets alongside the old ones, then move a stable alias — a collection-level “latest” link or a catalog search on the version property — rather than overwriting files in place. Consumers who pinned a version keep working, and nobody is silently served different numbers than they were last week.

---

Related

• Writing Prediction Rasters as COGs — profiles, palettes and overview policy in detail.
• Vectorizing Predicted Masks to Polygons — shapes, simplification and topology cleanup.
• Smoothing and Post-Processing Classification Rasters — majority filters, sieving and minimum mapping units.
• Attaching Model Metadata to STAC Items — the provenance fields worth standardising on.
• Writing and Validating Cloud-Optimized GeoTIFFs — the general COG authoring rules this topic applies.