Automating Batch Pre-Labeling with SAM and QGIS

A headless Python pipeline runs Meta’s Segment Anything Model (SAM) on georeferenced raster tiles, transforms binary segmentation masks back to map coordinates using each tile’s affine transform, and exports spatially aligned GeoJSON that QGIS loads directly for human review. The pipeline removes manual digitizing as the rate-limiting step: a GPU pass over one hundred 1024×1024 tiles completes in minutes, producing polygon annotations that annotators correct rather than draw from scratch.

Why Unprojected Inference Breaks the Pipeline

When the source raster is in a geographic coordinate reference system such as EPSG:4326, one pixel does not represent a consistent ground distance across the image. SAM’s receptive field spans different physical areas in the north versus south of a scene, causing fragmented or oversized masks at scale. The second failure is quieter: converting pixel-space contours back to geographic coordinates using a geographic CRS affine transform produces geometries that look correct in a GIS viewer but carry distorted area and perimeter values, which collapses IoU threshold calculations at model evaluation time. Both failures are silent — the pipeline finishes, the GeoJSON loads, but the annotations do not match the ground truth.

SAM Batch Pre-Labeling Pipeline Four pipeline stages shown left to right: Raster Tiling and CRS Normalization, SAM Batch Inference, Mask to Vector Conversion, and QGIS Review and Export. Raster Tiling CRS → metric 1024×1024 px 15–20% overlap SAM Inference vit_h / vit_l IoU ≥ 0.85 stability ≥ 0.92 Mask → Vector affine transform topology clean GeoJSON export QGIS Review confidence triage topology check COCO / YOLO out Step 1 Step 2 Step 3 Step 4
Four deterministic pipeline stages: tile and reproject, run SAM, convert masks to georeferenced polygons, then validate in QGIS.

Step-by-Step Implementation

Step 1 — Tile the Raster and Normalize the CRS

Use rasterio and rio-cogeo to reproject to a local UTM zone and split into overlapping tiles. Reprojection must happen before tiling so every tile inherits a consistent metric CRS.

python 
# requirements: rasterio==1.3.10, rio-cogeo==3.6.0, numpy==1.26.4
import math
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.windows import Window
import os

TARGET_CRS = "EPSG:32634"          # UTM zone 34N — adjust to your AOI
TILE_SIZE = 1024                   # pixels
OVERLAP = int(TILE_SIZE * 0.175)   # ~18% overlap


def reproject_to_metric(src_path: str, dst_path: str) -> None:
    """Reproject an arbitrary raster to a metric CRS before tiling."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, TARGET_CRS, src.width, src.height, *src.bounds
        )
        profile = src.profile.copy()
        profile.update(crs=TARGET_CRS, transform=transform,
                       width=width, height=height)
        with rasterio.open(dst_path, "w", **profile) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=TARGET_CRS,
                    resampling=Resampling.bilinear,
                )


def tile_raster(src_path: str, out_dir: str) -> list[str]:
    """Slice a metric raster into overlapping tiles; return tile paths."""
    os.makedirs(out_dir, exist_ok=True)
    tile_paths: list[str] = []
    step = TILE_SIZE - OVERLAP
    with rasterio.open(src_path) as src:
        cols = math.ceil((src.width - OVERLAP) / step)
        rows = math.ceil((src.height - OVERLAP) / step)
        for r in range(rows):
            for c in range(cols):
                col_off = c * step
                row_off = r * step
                w = min(TILE_SIZE, src.width - col_off)
                h = min(TILE_SIZE, src.height - row_off)
                window = Window(col_off, row_off, w, h)
                transform = src.window_transform(window)
                profile = src.profile.copy()
                profile.update(width=w, height=h, transform=transform)
                tile_path = os.path.join(out_dir, f"tile_r{r:04d}_c{c:04d}.tif")
                with rasterio.open(tile_path, "w", **profile) as dst:
                    dst.write(src.read(window=window))
                tile_paths.append(tile_path)
    return tile_paths

Step 2 — Batch SAM Inference

Pre-load the checkpoint once and iterate tiles in a torch.no_grad() context. Mixed-precision cuts VRAM usage by roughly 40% without affecting mask quality.

python 
# requirements: segment-anything==1.0, torch==2.3.0, opencv-python==4.9.0.80
import torch
import numpy as np
import rasterio
from segment_anything import sam_model_registry, SamAutomaticMaskGenerator

CHECKPOINT = "checkpoints/sam_vit_h_4b8939.pth"
MODEL_TYPE = "vit_h"
DEVICE = "cuda" if torch.cuda.is_available() else "cpu"

sam = sam_model_registry[MODEL_TYPE](checkpoint=CHECKPOINT)
sam.to(device=DEVICE)

mask_generator = SamAutomaticMaskGenerator(
    model=sam,
    points_per_side=32,
    pred_iou_thresh=0.85,
    stability_score_thresh=0.92,
    crop_n_layers=1,
    min_mask_region_area=500,   # pixels; filters noise fragments
)


def run_sam_on_tiles(tile_paths: list[str]) -> list[tuple[str, list[dict]]]:
    """Return (tile_path, masks) pairs for every input tile."""
    results: list[tuple[str, list[dict]]] = []
    with torch.no_grad():
        with torch.autocast(device_type=DEVICE, dtype=torch.float16):
            for tile_path in tile_paths:
                with rasterio.open(tile_path) as src:
                    img = src.read().transpose(1, 2, 0)  # (H, W, C)
                    # SAM expects uint8 RGB; handle multi-band imagery
                    if img.shape[2] > 3:
                        img = img[:, :, :3]
                    img = img.astype(np.uint8)
                masks = mask_generator.generate(img)
                results.append((tile_path, masks))
    return results

Step 3 — Convert Pixel Masks to Georeferenced Polygons

Apply each tile’s affine transform to map pixel contour coordinates to real-world map coordinates. Assign per-polygon confidence scores so annotators can triage masks by certainty in QGIS.

python 
# requirements: shapely==2.0.4, geopandas==0.14.4, opencv-python==4.9.0.80
import cv2
from shapely.geometry import Polygon, mapping
from shapely.ops import transform as shapely_transform
import geopandas as gpd
from rasterio.transform import Affine

MIN_AREA_PX = 500


def mask_to_georef_polygons(
    mask_binary: np.ndarray,
    affine_tfm: Affine,
    predicted_iou: float,
    source_tile: str,
) -> list[dict]:
    """Convert one SAM binary mask to a list of GeoJSON-ready feature dicts."""
    contours, _ = cv2.findContours(
        mask_binary.astype(np.uint8),
        cv2.RETR_EXTERNAL,
        cv2.CHAIN_APPROX_SIMPLE,
    )
    features: list[dict] = []
    for cnt in contours:
        if cv2.contourArea(cnt) < MIN_AREA_PX:
            continue
        pts = cnt.squeeze()
        if pts.ndim < 2 or len(pts) < 3:
            continue
        # Build pixel-space polygon, then apply raster affine transform
        pixel_ring = [tuple(pt) for pt in pts]
        pixel_poly = Polygon(pixel_ring)
        if not pixel_poly.is_valid:
            pixel_poly = pixel_poly.buffer(0)

        def _affine(xs, ys):
            geo_xs, geo_ys = zip(*(affine_tfm * (col, row) for col, row in zip(xs, ys)))
            return list(geo_xs), list(geo_ys)

        geo_poly = shapely_transform(_affine, pixel_poly)
        geo_poly = geo_poly.simplify(0.5, preserve_topology=True)
        features.append({
            "geometry": mapping(geo_poly),
            "properties": {
                "source_tile": source_tile,
                "confidence": round(float(predicted_iou), 4),
                "area_m2": round(geo_poly.area, 2),
            },
        })
    return features


def build_geodataframe(
    sam_results: list[tuple[str, list[dict]]],
    output_crs: str = TARGET_CRS,
) -> gpd.GeoDataFrame:
    """Aggregate all tile results into a single GeoDataFrame."""
    all_features: list[dict] = []
    for tile_path, masks in sam_results:
        with rasterio.open(tile_path) as src:
            tfm = src.transform
        for m in masks:
            all_features.extend(
                mask_to_georef_polygons(
                    m["segmentation"], tfm, m["predicted_iou"], tile_path
                )
            )
    gdf = gpd.GeoDataFrame.from_features(all_features, crs=output_crs)
    # Remove duplicates from overlapping tiles: keep the highest-confidence copy
    gdf = gdf.sort_values("confidence", ascending=False)
    gdf = gdf.drop_duplicates(subset=["geometry"])
    return gdf.reset_index(drop=True)

Step 4 — Export GeoJSON and Load into QGIS

Write the GeoDataFrame to a GeoJSON file and configure QGIS for annotation review.

python 
OUTPUT_GEOJSON = "prelabels/sam_batch_output.geojson"

def export_prelabels(gdf: gpd.GeoDataFrame) -> None:
    gdf.to_file(OUTPUT_GEOJSON, driver="GeoJSON")
    print(f"Exported {len(gdf)} pre-label polygons → {OUTPUT_GEOJSON}")
    high_conf = (gdf["confidence"] >= 0.90).sum()
    print(f"  High-confidence (≥0.90): {high_conf} ({100*high_conf/len(gdf):.1f}%)")

In QGIS, open Layer > Add Layer > Add Vector Layer and select the GeoJSON. Apply a Graduated renderer on the confidence field — green for ≥ 0.90, amber for 0.80–0.89, red for < 0.80. Enable Snapping (Settings > Snapping Options, tolerance 5–10 map units) so editors can snap pre-label edges to existing survey boundaries. This review cycle integrates with the broader QGIS plugin ecosystem for annotation teams.

Because every polygon already carries a metric geotransform, the reviewed layer exports cleanly to COCO/YOLO formats with the geotransform preserved as metadata — no second reprojection round-trip, and pixel-space bounding boxes recovered by inverting the same affine used in Step 3.

Key Parameters Reference

Parameter Recommended value Effect
points_per_side 32 Grid density for automatic prompt points; raise to 64 for dense urban scenes
pred_iou_thresh 0.85 Minimum predicted IoU to keep a mask; lower increases recall, raises noise
stability_score_thresh 0.92 Mask consistency across threshold perturbations; lower if SAM misses soft boundaries
min_mask_region_area 500 px Pixel noise filter; scale with GSD (higher GSD → larger threshold)
Tile overlap 15–20% Prevents edge-cut artifacts; overlap > 25% creates excessive duplicates
simplify tolerance 0.5 m Vertex reduction; keep < 1 m to preserve parcel or building edges
MIN_AREA_M2 post-filter 25–100 m² Remove sub-parcel noise after metric CRS conversion

Common Errors and Fixes

TypeError: image must be uint8

SAM requires three-channel uint8 RGB input, but multispectral rasters read via rasterio are often uint16 and carry more than three bands. The mask generator raises this on the first tile. Fix: img = img[:, :, :3].astype(np.uint8) — and apply a per-band stretch first if the source is true uint16 radiance rather than a display-ready composite.

Masks look correct but coordinates land far off the basemap

The raster’s CRS was geographic (EPSG:4326) at inference time, so each tile’s affine transform maps pixel columns and rows to decimal-degree deltas instead of metres. Polygons render with the right shape but the wrong scale and position. Fix: run reproject_to_metric() before tiling so every tile inherits a uniform metric transform.

Topology Checker flags hundreds of overlapping polygons

Tile overlap produces the same object in two adjacent tiles without deduplication. The drop_duplicates(subset=["geometry"]) step removes exact duplicates, but near-duplicates with slightly different vertex counts survive. Fix: dissolve near-duplicates with gpd.overlay(gdf, gdf, how="union") followed by a groupby("source_tile").geometry.union_all() pass, or snap to a shared grid before deduplicating.

CUDA out of memory with vit_h

The vit_h checkpoint needs roughly 7 GB of VRAM per batch, which overflows a 16 GB GPU once the automatic mask generator allocates its crop pyramid. Fix: switch to vit_l (~4 GB), enable torch.autocast as shown in Step 2, or drop crop_n_layers to 0 and process tiles strictly sequentially.

This page is part of the QGIS Plugin Ecosystem for Annotation Teams section within Labeling Workflows & Toolchain Integration.

Related in QGIS Plugin Ecosystem for Annotation Teams

Label Studio vs QGIS for Geospatial Annotation
Choose the tool by the geometry, not the brand. Use Label Studio for high-throughput distributed labeling of tiled image chips with model-assisted pre-labeling, where a large annotator pool clears bounding boxes and coarse masks at speed. Use QGIS for precision topology editing, vertex snapping, and CRS-aware cadastral