Vectorizing Predicted Masks to Polygons

Pass a mask so only the class you want is walked, and build geometries with the raster’s own transform:

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

with rasterio.open("class.tif") as src:
    arr = src.read(1)
    mask = arr == 3
    geoms = [shape(g) for g, v in shapes(arr, mask=mask, transform=src.transform) if v == 3]

gdf = gpd.GeoDataFrame({"class_id": 3}, geometry=geoms, crs=src.crs)

Omitting the mask is the difference between a clean layer and one enormous multi-part polygon full of holes. This page belongs to exporting and serving model outputs in Raster Machine Learning & Model Inference.

Why the Mask Argument Matters So Much

shapes walks every connected region of equal value, including the background. On a scene where one class covers 90% of the pixels, that background region is a single feature with tens of thousands of interior rings, and it is both useless and slow to write.

With and without a mask Without a mask the walk returns the three target patches plus one background polygon that surrounds them and carries three interior rings. With a mask it returns only the three patches. On a real scene the background feature can carry tens of thousands of rings and dominate both the file size and the write time. no mask — 4 features mask=arr==3 — 3 features the grey background is feature 4 with one interior ring per patch only the target patches are returned no rings, no background geometry

The second consequence is subtler: without a mask, the value returned by shapes includes the background code, so any downstream filter has to know about it. With a mask, every geometry that comes back is one you asked for.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 features.shapes and the transform
shapely >=2.0 Geometry construction, simplification, validity
geopandas >=0.14 Attributes and writing GeoPackage or Parquet
numpy >=1.23 Mask construction 
pip install "rasterio>=1.3.0" "shapely>=2.0" "geopandas>=0.14" "numpy>=1.23"

Complete Working Example

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

CLASS_NAMES = {1: "forest", 2: "cropland", 3: "built", 4: "water"}


def vectorize_prediction(path: str, *, classes: list[int] | None = None,
                         min_area_ha: float = 0.1,
                         simplify_px: float = 1.0) -> gpd.GeoDataFrame:
    """Polygonise selected classes from a prediction raster."""
    records = []
    with rasterio.open(path) as src:
        arr = src.read(1)
        px_area_ha = abs(src.transform.a * src.transform.e) / 10_000
        tol = simplify_px * abs(src.transform.a)
        crs = src.crs

        for code in (classes or sorted(CLASS_NAMES)):
            mask = arr == code
            if not mask.any():
                continue
            for geom, value in shapes(arr, mask=mask, transform=src.transform):
                if int(value) != code:
                    continue
                poly = shape(geom)
                records.append({
                    "class_id": code,
                    "class_name": CLASS_NAMES.get(code, str(code)),
                    # pixel count BEFORE simplification keeps area consistent with the map
                    "pixels": int(round(poly.area / abs(src.transform.a * src.transform.e))),
                    "geometry": poly.simplify(tol, preserve_topology=True),
                })

    gdf = gpd.GeoDataFrame(records, crs=crs)
    gdf["area_ha"] = gdf["pixels"] * px_area_ha
    gdf = gdf[gdf["area_ha"] >= min_area_ha].reset_index(drop=True)
    gdf["geometry"] = gdf.geometry.make_valid()
    return gdf


if __name__ == "__main__":
    layer = vectorize_prediction("class.tif", classes=[3, 4])
    print(layer.groupby("class_name").agg(n=("area_ha", "size"),
                                          total_ha=("area_ha", "sum")))
    layer.to_file("prediction_polygons.gpkg", layer="predictions", driver="GPKG")

Recording the pixel count before simplification is the detail that keeps the vector layer honest. Simplification moves vertices, so polygon area drifts by a fraction of a per cent per feature, and a total computed from simplified geometry will not match the raster the layer was derived from. Storing the count means both numbers are available and the discrepancy is visible rather than mysterious.

Choosing the Simplification Tolerance

Tolerance against fidelity At zero tolerance the boundary is a pixel staircase with hundreds of vertices. At one pixel it becomes a smooth outline with a tenth of the vertices and almost no area change. At four pixels corners are cut, the area changes by several per cent and small features are lost entirely. tolerance 0 tolerance 1 px tolerance 4 px 412 vertices area exact 46 vertices area within 0.3% 11 vertices area off by 4%; corners lost One pixel is the tolerance that removes the artefact without removing the shape.

The rule of thumb is that the tolerance should equal the pixel size, because the staircase is exactly a pixel-scale artefact and nothing smaller than a pixel was ever real. Going further is a cartographic decision rather than a data one, and it should be made at display time rather than baked into the stored layer.

preserve_topology=True costs a little speed and prevents the failure that makes a simplified layer unusable: self-intersections created when a narrow neck collapses. Even with it, run make_valid afterwards, because simplification of a polygon with holes can still produce a ring that touches its own exterior.

For very large layers, simplification is also the main lever on file size. A national built-up layer at zero tolerance can run to hundreds of millions of vertices; at one pixel it is an order of magnitude smaller and renders interactively in any desktop GIS.

Verification

import numpy as np
import rasterio

with rasterio.open("class.tif") as src:
    arr = src.read(1)
    px_area_ha = abs(src.transform.a * src.transform.e) / 10_000

for code, group in layer.groupby("class_id"):
    raster_ha = (arr == code).sum() * px_area_ha
    vector_ha = group["area_ha"].sum()
    print(f"class {code}: raster {raster_ha:,.0f} ha, vector {vector_ha:,.0f} ha")
    assert abs(vector_ha - raster_ha) / raster_ha < 0.05, "area drifted more than 5%"

assert layer.geometry.is_valid.all(), "invalid geometries in the output"
assert not layer.geometry.is_empty.any()
Raster area against vector area per class For three classes the raster and vector totals agree closely, with the vector slightly lower because the minimum mapping unit removed small patches. A fourth class shows a much larger gap, which indicates a problem in the mask or an over-aggressive tolerance rather than the expected filtering loss. Vector totals should trail raster totals slightly, not sharply forest -2% cropland -3% built -38% — investigate raster pixel area vector polygon area

The area comparison should differ only by whatever the minimum mapping unit removed. A larger gap means either the tolerance is too aggressive or a class was partly dropped by the mask — both are easy to fix and impossible to notice once the layer is in someone else’s hands. The counterpart check on the raster side is described in computing zonal statistics with rasterstats.

Common Errors

TopologyException when writing the layer

Simplification produced a self-intersecting ring. Call make_valid() after simplifying and before writing, and keep preserve_topology=True.

The polygon count is in the millions

Speckle survived into the vectorization. Apply the majority filter and minimum mapping unit first, as in smoothing and post-processing classification rasters; a cleaned raster typically vectorizes to a hundredth as many features.

Polygons are offset from the imagery

The transform passed to shapes was not the transform of the array being walked — usually a window array with the scene transform. Use src.window_transform(win) for windowed work.

Frequently Asked Questions

Q: Why does vectorizing produce one giant polygon? Because no mask was passed, so the background was walked as a single feature with a hole for every object. Pass mask= to shapes with a boolean array of the pixels you actually want, and the background disappears from the output entirely.

Q: How much should the polygons be simplified? About one pixel. That removes the staircase edges without moving any boundary further than the imagery can resolve. Larger tolerances start cutting corners off real shapes, and the area attribute stops matching the raster it came from.

Q: Should area be computed from the raster or the polygon? From the raster if the number must be consistent with the map, because simplification changes polygon area by a fraction of a per cent per feature. Carry a pixel count attribute alongside the geometry so both are available and the difference is explicit.

Q: What format should the polygon layer ship in? GeoPackage for general use and GeoParquet where the consumer is a data pipeline rather than a desktop GIS. Both carry the CRS properly, neither truncates field names, and both handle layers of millions of features — which shapefiles, still the default in some organisations, do not.