Smoothing and Post-Processing Classification Rasters

Filter, sieve, then restore the footprint — in that order:

import numpy as np
from scipy import ndimage

nodata = classes == 0
smoothed = majority_filter(classes, size=3)
cleaned = sieve_small_patches(smoothed, min_pixels=10)
cleaned[nodata] = 0                    # never extend the map past its evidence

Every one of those three lines prevents a defect that ships otherwise. This page belongs to exporting and serving model outputs in Raster Machine Learning & Model Inference.

Speckle Is a Property of the Model, Not the Landscape

A pixel-based classifier makes each decision independently, so nothing in it prefers spatially coherent output. Near a decision boundary — a field that is marginally more soil than crop this week — small variations in reflectance flip individual pixels, and the result is a scatter of isolated values inside otherwise uniform regions.

The important judgement is which single pixels are noise and which are real. A farm building inside a field is a genuine one-pixel feature at 10 m; a single “water” pixel in the middle of a wheat field is not. Because a filter cannot tell them apart, the minimum mapping unit should be a documented product decision rather than a tuning parameter.

Environment & Setup

pip install "numpy>=1.23" "scipy>=1.11" "rasterio>=1.3.0"

Complete Working Example

import numpy as np
import rasterio
from rasterio.features import sieve
from scipy import ndimage


def majority_filter(classes: np.ndarray, size: int = 3) -> np.ndarray:
    """Most common class in each neighbourhood, computed per class in one pass each."""
    codes = np.unique(classes[classes > 0])
    if codes.size == 0:
        return classes
    votes = np.zeros((codes.max() + 1, *classes.shape), dtype="float32")
    for code in codes:
        votes[code] = ndimage.uniform_filter(
            (classes == code).astype("float32"), size=size, mode="nearest")
    out = votes.argmax(axis=0).astype("uint8")
    return np.where(classes == 0, 0, out)          # nodata is never filled


def clean_prediction(in_path: str, out_path: str, *,
                     filter_size: int = 3, min_pixels: int = 10) -> dict[str, int]:
    with rasterio.open(in_path) as src:
        classes = src.read(1)
        profile = src.profile
        colormap = src.colormap(1) if src.colormap(1) else None

    nodata = classes == 0
    smoothed = majority_filter(classes, size=filter_size)
    # rasterio's sieve dissolves connected regions below min_pixels into neighbours
    cleaned = sieve(smoothed, size=min_pixels, connectivity=8)
    cleaned[nodata] = 0

    changed = int((cleaned != classes).sum())
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(cleaned, 1)
        if colormap:
            dst.write_colormap(1, colormap)
        dst.update_tags(post_processing=f"majority{filter_size};sieve{min_pixels}")

    return {"changed_pixels": changed,
            "changed_fraction_pct": round(100 * changed / classes.size, 2)}

rasterio.features.sieve is worth preferring over a hand-rolled connected-component pass: it is implemented in GDAL, handles the neighbour-assignment rule properly, and is an order of magnitude faster on a full scene. Writing the operations into the file’s tags means a consumer can tell what was done without asking.

Probability-Aware Smoothing

Where the probability band from inference is available, a better cleanup is possible: smooth the probabilities and re-take the argmax, so confident pixels resist being overruled by their neighbours and uncertain ones give way easily.

import numpy as np
from scipy import ndimage


def smooth_probabilities(probs: np.ndarray, sigma: float = 1.0) -> np.ndarray:
    """Gaussian-smooth each class probability plane, then re-decide."""
    smoothed = np.stack([ndimage.gaussian_filter(p, sigma=sigma) for p in probs])
    return smoothed.argmax(axis=0).astype("uint8") + 1

A sigma of about one pixel is a good starting point and corresponds roughly to a 3×3 majority filter in strength, while behaving far better at boundaries: because the underlying surface is continuous, an edge moves smoothly rather than jumping a whole pixel when the vote tips.

The one caveat is that Gaussian smoothing of probabilities does not preserve their sum exactly at the edges of the valid region, so re-normalise if the probabilities themselves are being published rather than only the argmax.

Verification

import numpy as np

assert set(np.unique(cleaned).tolist()) <= set(np.unique(classes).tolist()), \
    "cleanup introduced a class code that was not predicted"
assert ((cleaned == 0) == (classes == 0)).all(), "nodata footprint changed"

frac = (cleaned != classes).mean()
print(f"{frac:.2%} of pixels changed")
assert frac < 0.15, "cleanup is rewriting too much of the map"

The changed-fraction bound is the guard that catches an over-aggressive filter. A healthy 3×3 majority pass on a moderately speckled prediction changes two to six per cent of pixels; anything above ten per cent means either the model output is extremely noisy — worth fixing upstream rather than downstream — or the window is far too large. The footprint assertion catches the other classic failure, where a filter quietly fills nodata holes with whatever surrounds them and the published map covers ground the imagery never saw. That is worse than a wrong class, because a consumer has no way to tell the invented area from the observed area once the file has left your hands.

Common Errors

Nodata holes are filled after filtering

The np.where(classes == 0, 0, out) guard is missing. A filter has no concept of “no observation” and will happily interpolate across it.

Small real features disappear

The minimum mapping unit is larger than the features you care about. It is a product decision: state it, and if single-pixel detections matter, do not sieve at all.

The file grows after cleaning

Impossible for the raster, so the file was rewritten without predictor=2 or without tiling. Check the profile, as described in writing prediction rasters as COGs.

Frequently Asked Questions

Q: Does smoothing improve accuracy or just appearance? Usually both, by a small amount. Isolated misclassified pixels are more often wrong than right, so removing them raises overall accuracy slightly while removing far more of the visual noise. It lowers accuracy only where the real landscape genuinely contains single-pixel features.

Q: Majority filter or probability smoothing? Probability smoothing is better where probabilities exist, because it respects how confident the model was rather than treating every pixel’s vote equally. A majority filter is the only option when all you have is the class raster.

Q: What order should the operations run in? Smooth first, sieve second, restore nodata last. Sieving before smoothing wastes work on patches the filter would have removed anyway, and restoring nodata at any point other than the end lets a later step spill across the boundary.

Q: Should the unsmoothed raster be kept? Keep it as provenance even if it is never published. It is the only way to re-derive a different minimum mapping unit later without re-running inference, and under compression it costs little next to the imagery it came from.

Related

• Exporting and Serving Model Outputs — the publishing workflow this cleanup precedes.
• Writing Prediction Rasters as COGs — writing the cleaned result properly.
• Vectorizing Predicted Masks to Polygons — the step that benefits most from a cleaned raster.
• Thresholding Change Maps and Removing Noise — the same filters applied to change products.