297 lines
12 KiB
Python
297 lines
12 KiB
Python
#!/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 touching/overlapping footprints 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)
|
|
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],
|
|
}
|
|
)
|
|
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_cutoff = (
|
|
datetime.fromisoformat(imagery_valid_to.replace("Z", "+00:00")) if imagery_valid_to 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,
|
|
"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())
|