Calculating IoU Thresholds for Geospatial Object Detection

Reliable IoU evaluation for aerial and satellite object detection requires three things that standard computer vision toolkits omit: reprojection from angular coordinates to a local metric coordinate reference system before any area computation, topology validation with shapely 2.0+ make_valid(), and threshold values calibrated to ground sample distance (GSD) and mission type. Reproject all geometries to a UTM or state-plane CRS using pyproj 3.6+, then apply adaptive cutoffs in the 0.35–0.75 range rather than a fixed 0.50. Skipping any of these steps introduces projection-induced bias that silently corrupts precision/recall curves.

# Why Unprojected IoU Breaks Aerial Pipelines

EPSG:4326 (WGS84) stores coordinates as decimal degrees; a degree of longitude spans roughly 111 km at the equator but shrinks near zero at high latitudes, so intersection area computed in degree-squared units is geographically meaningless and produces artificially suppressed IoU scores that shift as your dataset spans different latitudes. Varying GSD compounds the problem: an annotation tolerance of ±2 pixels at 10 cm/pixel is ±20 cm on the ground, but at 50 cm/pixel that same pixel tolerance is ±100 cm — a fixed 0.50 threshold that passes valid detections at high resolution will reject them at coarser resolution purely because boundary pixelation widens the mismatch, not because the model degraded.

# Step-by-Step Implementation

Install the required packages once:

pip install shapely==2.0.6 pyproj==3.6.1 numpy==1.26.4

# Step 1 — Reproject All Geometries to a Metric CRS

Choose a UTM zone that covers your tile’s centroid. For a dataset at longitude 13.4°E (central Europe), EPSG:32633 (UTM Zone 33N) is appropriate. Pass always_xy=True to force (longitude, latitude) input order regardless of CRS authority axis definitions:

from pyproj import Transformer

def make_transformer(source_crs: str = "EPSG:4326",
                     target_crs: str = "EPSG:32633") -> Transformer:
    """Return a Transformer that always expects (lon, lat) / (x, y) input order."""
    return Transformer.from_crs(source_crs, target_crs, always_xy=True)

For global datasets spanning multiple UTM zones, derive the zone automatically from each annotation’s centroid longitude:

def utm_epsg_from_lon_lat(lon: float, lat: float) -> str:
    zone = int((lon + 180) / 6) + 1
    hemisphere = "326" if lat >= 0 else "327"
    return f"EPSG:{hemisphere}{zone:02d}"

# Step 2 — Validate Topology Before Computing Area

shapely’s make_valid() repairs self-intersecting rings and unclosed exteriors. Call it on every geometry — prediction and ground-truth alike — before any set operation:

from shapely.geometry import Polygon, box
from shapely.validation import make_valid
from typing import Union

def to_valid_polygon(coords: Union[list, tuple]) -> Polygon:
    """
    Accept either [minx, miny, maxx, maxy] (axis-aligned box)
    or a list of (x, y) ring coordinates (arbitrary polygon).
    Always returns a topologically valid shapely Polygon.
    """
    if len(coords) == 4 and not isinstance(coords[0], (list, tuple)):
        geom = box(*coords)
    else:
        geom = Polygon(coords)
    return make_valid(geom)

# Step 3 — Compute IoU in Metric Space

After reprojection and validation, intersection and union areas are in square metres. Assign confidence scores alongside the IoU value when building evaluation logs — they let you weight borderline matches rather than applying a hard binary cut:

from shapely.ops import transform
from shapely.geometry import Polygon

def geospatial_iou(
    pred: Polygon,
    gt: Polygon,
    transformer: "Transformer",
) -> float:
    """
    Compute IoU between two valid shapely Polygons after projecting
    both to the metric CRS defined by `transformer`.
    Returns a float in [0.0, 1.0].
    """
    project = lambda x, y: transformer.transform(x, y)
    pred_m = transform(project, pred)
    gt_m   = transform(project, gt)

    intersection_area = pred_m.intersection(gt_m).area
    union_area        = pred_m.union(gt_m).area

    if union_area == 0.0:
        return 0.0
    return intersection_area / union_area

# Step 4 — Apply a GSD-Calibrated Threshold

Wrap steps 1–3 in a single evaluation function that accepts explicit CRS and threshold parameters. Consult the reference table in the next section to choose iou_threshold:

def evaluate_detection(
    pred_coords: list,
    gt_coords: list,
    source_crs: str = "EPSG:4326",
    target_crs: str = "EPSG:32633",
    iou_threshold: float = 0.50,
) -> tuple[float, bool]:
    """
    Full projection-aware IoU evaluation.

    Args:
        pred_coords: [minx,miny,maxx,maxy] or polygon ring [(lon,lat), …]
        gt_coords:   same format as pred_coords
        source_crs:  CRS of the input coordinates (default WGS84)
        target_crs:  Local metric CRS for area computation (default UTM 33N)
        iou_threshold: Match cutoff, calibrated to GSD and mission type

    Returns:
        (iou_score, is_match)
    """
    t = make_transformer(source_crs, target_crs)
    pred_poly = to_valid_polygon(pred_coords)
    gt_poly   = to_valid_polygon(gt_coords)
    iou       = geospatial_iou(pred_poly, gt_poly, t)
    return iou, iou >= iou_threshold

# Step 5 — Size-Stratified Batch Evaluation

Aggregate IoU scores mask scale-dependent failure modes. Bin predictions by projected area — the same approach used in polygon vs. bounding-box annotation quality assessment — and compute per-bin match rates after transformation:

import numpy as np
from dataclasses import dataclass

@dataclass
class DetectionResult:
    iou: float
    projected_area_m2: float   # ground-truth area in square metres
    is_match: bool

def stratified_map(
    results: list[DetectionResult],
    bins: dict[str, tuple[float, float]] | None = None,
) -> dict[str, float]:
    """
    Compute per-size-bin match rate (proxy mAP) from a list of DetectionResults.
    Default bins: small < 100 m², medium 100–10 000 m², large > 10 000 m².
    """
    if bins is None:
        bins = {"small": (0, 100), "medium": (100, 10_000), "large": (10_000, float("inf"))}

    summary: dict[str, float] = {}
    for name, (lo, hi) in bins.items():
        subset = [r for r in results if lo <= r.projected_area_m2 < hi]
        if subset:
            summary[name] = float(np.mean([r.is_match for r in subset]))
        else:
            summary[name] = float("nan")
    return summary

# Threshold and CRS Reference

IoU threshold by mission type:

GSD scaling rule: lower the threshold by approximately 0.05 per 10 cm/pixel increase in GSD above 20 cm/pixel. At 50 cm/pixel, sub-pixel annotation disagreement between labellers spans 25–50 cm on the ground — more than enough to drop IoU below a fixed 0.50 for a correctly localised detection.

EPSG quick reference:

• EPSG:4326 — WGS84, angular degrees, input format only
• EPSG:32633 — UTM Zone 33N, metric (central Europe / Africa)
• EPSG:32737 — UTM Zone 37S, metric (East Africa / Madagascar)
• EPSG:3857 — Web Mercator, metric fallback for multi-zone datasets (distortion acceptable for tile-level area ratios under ~50 km wide)

Ensure the CRS and GSD metadata needed for threshold selection travels with every dataset export — preserving metadata across dataset versions covers the mechanics of attaching that context to each versioned snapshot.

# Common Errors and Fixes

TopologicalError: The operation 'GEOSIntersection_r' produced a null geometry Root cause: self-intersecting polygon (bowtie ring) passed to .intersection(). Fix: call make_valid(geom) on both operands before the intersection call.

IoU is always 0.0 for visually overlapping boxes Root cause: input coordinates in (lat, lon) order passed to a transformer expecting (lon, lat). Fix: add always_xy=True to Transformer.from_crs().

IoU scores drop sharply for tiles above 55°N Root cause: area computed in degree-squared units — EPSG:4326 was never reprojected. Fix: ensure to_valid_polygon receives metric coordinates after transform(project, geom) has been applied.

ShapelyDeprecationWarning: The array interface is deprecated / wrong area returned Root cause: shapely 1.x geometry object passed to a shapely 2.x function. Fix: upgrade to shapely==2.0.6 and re-create all geometry objects from raw coordinates rather than unpickling from 1.x.

# Related

• Coordinate Reference Systems in Annotation Pipelines — parent page covering CRS contracts, datum management, and reprojection patterns across an entire annotation pipeline
• Confidence Scoring for Geospatial Labels — assign per-annotation uncertainty values that complement IoU during model evaluation and active-learning triage
• Best Practices for Polygon vs Bounding Box Annotation — annotation geometry choices that determine how tight IoU scores can realistically be
• Preserving Metadata Across Dataset Versions — ensure the CRS and GSD metadata needed for threshold selection travel with every dataset export

This page covers one specialised calculation within Coordinate Reference Systems in Annotation Pipelines, which is itself part of Geospatial Annotation Fundamentals & Architecture.