#!/usr/bin/env python3 """Normalize governed Belgian building references for detector training. The script never asserts semantic parity between providers. It emits one canonical training class while retaining provider-native identifiers, properties and an explicit accept/reject decision for every source feature. """ from __future__ import annotations import argparse import hashlib import json from collections import Counter from datetime import UTC, datetime from pathlib import Path from typing import Any import rasterio from pyproj import Transformer from shapely.geometry import mapping, shape from shapely.ops import transform as shapely_transform, unary_union from shapely.validation import make_valid SUPPORTED_SOURCES = {"grb", "spw_picc", "urbis"} EXCLUDED_TOKENS = { "canopy": ("canopy", "afdak", "auvent"), "ruin": ("ruin", "ruine", "ruïne"), "underground": ("underground", "ondergronds", "souterrain"), } def _feature_id(feature: dict[str, Any], index: int) -> str: properties = feature.get("properties") or {} for key in ("source_feature_id", "OBJECTID", "INSPIRE_ID", "id", "gml_id"): value = properties.get(key) if value not in (None, ""): return str(value) if feature.get("id") not in (None, ""): return str(feature["id"]) return f"row-{index}" def _source_class(properties: dict[str, Any]) -> str | None: for key in ("source_class", "TYPE", "type", "OBJTYPE", "nature", "NATURE", "class"): value = properties.get(key) if value not in (None, ""): return str(value) return None def _semantic_exclusion(properties: dict[str, Any]) -> str | None: haystack = " ".join(str(value).lower() for value in properties.values() if value is not None) for reason, tokens in EXCLUDED_TOKENS.items(): if any(token in haystack for token in tokens): return f"excluded_{reason}" return None def _source_creation_at(source_name: str, properties: dict[str, Any]) -> datetime | None: raw: Any = None if source_name == "grb": raw = properties.get("BEGINDATUM") elif source_name == "spw_picc": raw = properties.get("DATE_CREAT") if raw in (None, ""): return None try: if isinstance(raw, (int, float)) or str(raw).isdigit(): value = float(raw) if value > 10_000_000_000: value /= 1000 return datetime.fromtimestamp(value, tz=UTC) parsed = datetime.fromisoformat(str(raw).replace("Z", "+00:00")) return parsed if parsed.tzinfo else parsed.replace(tzinfo=UTC) except (OSError, OverflowError, ValueError): return None def _polygonal(geometry: Any) -> Any | None: if geometry.geom_type in {"Polygon", "MultiPolygon"}: return geometry if geometry.geom_type == "GeometryCollection": polygons = [part for part in geometry.geoms if part.geom_type in {"Polygon", "MultiPolygon"}] if not polygons: return None from shapely.ops import unary_union return unary_union(polygons) return None def merge_touching_roof_instances(features: list[dict[str, Any]], source_name: str) -> list[dict[str, Any]]: """Dissolve only compact touching groups into imagery-visible roof instances.""" if not features: return [] source_geometries = [(feature, shape(feature["geometry"])) for feature in features] dissolved = unary_union([geometry for _feature, geometry in source_geometries]) components = list(dissolved.geoms) if dissolved.geom_type == "MultiPolygon" else [dissolved] merged: list[dict[str, Any]] = [] for component in components: contributors = [ feature for feature, geometry in source_geometries if geometry.intersects(component) ] source_ids = sorted(str(feature["properties"]["source_feature_id"]) for feature in contributors) envelope_area = component.envelope.area fill_ratio = component.area / envelope_area if envelope_area else 0.0 # Large connected blocks and irregular chains are administratively # adjacent but not one reliably box-shaped roof target. Preserve their # native instances instead of creating a giant ambiguous detector box. if len(contributors) > 12 or fill_ratio < 0.55: for feature in contributors: retained = json.loads(json.dumps(feature)) retained["properties"]["label_semantics"] = "native_instance_complex_touch_group" merged.append(retained) continue properties = dict(contributors[0]["properties"]) properties.update( { "label_semantics": "visible_touching_roof_instance", "source_feature_ids": source_ids, "source_feature_count": len(source_ids), "source_feature_id": source_ids[0], "roof_group_fill_ratio": fill_ratio, } ) digest = hashlib.sha256(component.normalize().wkb).hexdigest()[:24] merged.append( { "type": "Feature", "id": f"{source_name}:roof-instance:{digest}", "properties": properties, "geometry": mapping(component), } ) return merged def normalize( *, reference_path: Path, raster_path: Path, source_name: str, min_label_px: float, imagery_observed_at: str | None, reference_observed_at: str | None, imagery_valid_to: str | None = None, merge_touching_roofs: bool = False, ) -> tuple[dict[str, Any], dict[str, Any]]: if source_name not in SUPPORTED_SOURCES: raise SystemExit(f"Unsupported governed building source: {source_name}") payload = json.loads(reference_path.read_text(encoding="utf-8-sig")) if payload.get("type") != "FeatureCollection" or not isinstance(payload.get("features"), list): raise SystemExit("Reference must be a GeoJSON FeatureCollection") accepted: list[dict[str, Any]] = [] decisions: list[dict[str, Any]] = [] seen: set[str] = set() counts: Counter[str] = Counter() imagery_dates = [ datetime.fromisoformat(value.replace("Z", "+00:00")) for value in (imagery_observed_at, imagery_valid_to) if value ] # Annual mosaics expose a validity interval but not the flight date for # each pixel. Using the interval end admits buildings created later in the # same year that are visibly absent from the mosaic. The earliest governed # date is therefore the only conservative training-label cutoff. imagery_cutoff = min(imagery_dates) if imagery_dates else None with rasterio.open(raster_path) as raster: if raster.crs is None: raise SystemExit("Raster CRS is required") transformer = Transformer.from_crs("EPSG:4326", raster.crs, always_xy=True) for index, feature in enumerate(payload["features"]): properties = dict(feature.get("properties") or {}) source_feature_id = _feature_id(feature, index) decision = { "source_name": source_name, "source_feature_id": source_feature_id, "source_class": _source_class(properties), "accepted": False, "reason": None, "geometry_repaired": False, } source_creation_at = _source_creation_at(source_name, properties) decision["source_creation_at"] = source_creation_at.isoformat() if source_creation_at else None try: geometry = shape(feature.get("geometry")) except Exception: decision["reason"] = "invalid_geometry_unreadable" decisions.append(decision) counts[decision["reason"]] += 1 continue if not geometry.is_valid: geometry = make_valid(geometry) decision["geometry_repaired"] = True geometry = _polygonal(geometry) if geometry is None or geometry.is_empty or not geometry.is_valid: decision["reason"] = "invalid_geometry_unrepairable" elif (reason := _semantic_exclusion(properties)) is not None: decision["reason"] = reason elif imagery_cutoff and source_creation_at and source_creation_at > imagery_cutoff: decision["reason"] = "created_after_imagery_period" else: metric = shapely_transform(transformer.transform, geometry) min_x, min_y, max_x, max_y = metric.bounds pixel_width = (max_x - min_x) / abs(float(raster.transform.a)) pixel_height = (max_y - min_y) / abs(float(raster.transform.e)) decision["pixel_width"] = pixel_width decision["pixel_height"] = pixel_height if pixel_width < min_label_px or pixel_height < min_label_px: decision["reason"] = "below_resolvable_pixel_size" else: digest = hashlib.sha256(geometry.normalize().wkb).hexdigest() if digest in seen: decision["reason"] = "duplicate_geometry" else: seen.add(digest) decision["accepted"] = True decision["reason"] = "accepted" normalized_properties = dict(properties) normalized_properties.update( { "canonical_class": "building", "source_name": source_name, "reference_layer_name": "buildings", "source_feature_id": source_feature_id, "source_class": decision["source_class"], "label_decision": "accepted", "geometry_repaired": decision["geometry_repaired"], } ) accepted.append( { "type": "Feature", "id": f"{source_name}:{source_feature_id}", "properties": normalized_properties, "geometry": mapping(geometry), } ) decisions.append(decision) counts[str(decision["reason"])] += 1 temporal_mismatch_days = None temporal_alignment_status = "unknown" if imagery_observed_at and reference_observed_at: imagery_date = datetime.fromisoformat(imagery_observed_at.replace("Z", "+00:00")) reference_date = datetime.fromisoformat(reference_observed_at.replace("Z", "+00:00")) temporal_mismatch_days = abs((imagery_date - reference_date).days) temporal_alignment_status = "measured" normalized_features = ( merge_touching_roof_instances(accepted, source_name) if merge_touching_roofs else accepted ) normalized = { "type": "FeatureCollection", "name": f"canonical-building-{source_name}", "features": normalized_features, } audit = { "schema_version": 1, "canonical_class": "building", "source_name": source_name, "reference_path": str(reference_path), "raster_path": str(raster_path), "min_label_px": min_label_px, "imagery_observed_at": imagery_observed_at, "imagery_valid_to": imagery_valid_to, "imagery_feature_creation_cutoff": imagery_cutoff.isoformat() if imagery_cutoff else None, "imagery_feature_creation_cutoff_policy": "earliest_governed_imagery_date", "reference_observed_at": reference_observed_at, "temporal_mismatch_days": temporal_mismatch_days, "temporal_alignment_status": temporal_alignment_status, "input_feature_count": len(payload["features"]), "accepted_source_feature_count": len(accepted), "accepted_feature_count": len(normalized_features), "merge_touching_roofs": merge_touching_roofs, "decision_counts": dict(sorted(counts.items())), "decisions": decisions, } return normalized, audit def main() -> int: parser = argparse.ArgumentParser() parser.add_argument("--reference", type=Path, required=True) parser.add_argument("--raster", type=Path, required=True) parser.add_argument("--source-name", choices=sorted(SUPPORTED_SOURCES), required=True) parser.add_argument("--output-reference", type=Path, required=True) parser.add_argument("--output-audit", type=Path, required=True) parser.add_argument("--min-label-px", type=float, default=3.0) parser.add_argument("--imagery-observed-at") parser.add_argument("--reference-observed-at") parser.add_argument("--imagery-valid-to") parser.add_argument("--merge-touching-roofs", action="store_true") args = parser.parse_args() normalized, audit = normalize( reference_path=args.reference, raster_path=args.raster, source_name=args.source_name, min_label_px=args.min_label_px, imagery_observed_at=args.imagery_observed_at, reference_observed_at=args.reference_observed_at, imagery_valid_to=args.imagery_valid_to, merge_touching_roofs=args.merge_touching_roofs, ) args.output_reference.parent.mkdir(parents=True, exist_ok=True) args.output_audit.parent.mkdir(parents=True, exist_ok=True) args.output_reference.write_text(json.dumps(normalized, ensure_ascii=False), encoding="utf-8") args.output_audit.write_text(json.dumps(audit, ensure_ascii=False, indent=2), encoding="utf-8") print(json.dumps({key: value for key, value in audit.items() if key != "decisions"}, ensure_ascii=False)) return 0 if __name__ == "__main__": raise SystemExit(main())