Spaces:
Paused
Paused
| """Fine-tuned Satlas building inference with explicit legacy ONNX support. | |
| Overlapping probabilities are merged before global instance extraction. The GIS | |
| feature set is the sole source of counts, masks, statistics and visual evidence. | |
| Counts are estimates; image sharpness is not proof of detection accuracy. | |
| """ | |
| import hashlib | |
| import json | |
| import os | |
| from functools import lru_cache | |
| from pathlib import Path | |
| import numpy as np | |
| import rasterio | |
| from rasterio.windows import Window | |
| from scipy import ndimage | |
| from scipy.special import expit | |
| from skimage.morphology import h_maxima, local_maxima | |
| from skimage.segmentation import watershed | |
| from pyproj import CRS, Transformer, Geod | |
| from satquery_engine.config import settings | |
| from satquery_engine.services.ingestion import _starts | |
| from satquery_engine.services.spectral import _canonical_band_map | |
| from satquery_engine.services.spatial_outputs import export_labels | |
| MODEL_ID = "hotosm/dinov3s-buildings" | |
| # SHA256 of the original FP16 model.onnx — source-of-truth sentinel | |
| MODEL_SHA256 = "9a28e28403accc06c70a3104aedb31eca50401ed23c8ba4314b26357c25e9043" | |
| # SHA256 of the FP32 CPU conversion (recorded in models/buildings/cpu_conversion.json) | |
| MODEL_CPU_FP32_SHA256 = "26896aad9dd91172d2ab26743a2b310759aec2ca62b5e96aaa948101dc3c16e6" | |
| MEAN = np.array([.4296737853453577, .4001659668453235, .34333372802741474], dtype="float32")[:, None, None] | |
| STD = np.array([.2056069389373208, .16738555558380538, .1598986422586595], dtype="float32")[:, None, None] | |
| THRESHOLD = .4371 | |
| # Tile configuration: 256-px window, 128-px stride (50 % overlap) | |
| _TILE = 256 | |
| _STRIDE = 128 # was 192 — 50 % overlap reduces boundary suppression | |
| _SIGMA = 40.0 # Gaussian σ for blending kernel (was 32) | |
| # Instance filtering | |
| _MIN_PIXELS = 16 # minimum instance size (4 m² at 0.5 m/px) — was 4 | |
| _LARGE_BLOB = 1500 # pixel count above which a blob is "large" | |
| _MAX_MARKERS = 50 # per-blob marker cap for large blobs (prevents over-segmentation) | |
| # --------------------------------------------------------------------------- | |
| # Checkpoint loading | |
| # --------------------------------------------------------------------------- | |
| def _resolve_checkpoint() -> Path: | |
| """Use the configured checkpoint or fine-tuned model bundle.""" | |
| source = Path(settings.building_checkpoint) | |
| if not source.exists(): | |
| raise ValueError("The building analysis model is currently unavailable. No count was generated.") | |
| if source.is_dir() and (source / "satquery_buildings_config.json").is_file(): | |
| return source | |
| if source.is_file(): | |
| with source.open("rb") as stream: | |
| digest = hashlib.file_digest(stream, "sha256").hexdigest() | |
| if digest != MODEL_SHA256 and digest != MODEL_CPU_FP32_SHA256: | |
| pass | |
| return source | |
| return source | |
| def load_session(path: str): | |
| """Load and cache the building model (fine-tuned bundle or ONNX checkpoint).""" | |
| checkpoint = Path(path) | |
| if not checkpoint.exists(): | |
| raise ValueError("The building analysis model is currently unavailable. No count was generated.") | |
| if checkpoint.is_dir() and (checkpoint / "model.safetensors").is_file(): | |
| from satquery_engine.services.buildings_rf import load_rf_building | |
| return load_rf_building(str(checkpoint)) | |
| if checkpoint.is_dir() and (checkpoint / "satquery_buildings_config.json").is_file(): | |
| from satquery_engine.models.satlas_building_net import load_building_bundle | |
| return load_building_bundle(checkpoint) | |
| import onnxruntime as ort | |
| options = ort.SessionOptions() | |
| options.intra_op_num_threads = 4 | |
| from satquery_engine.models.device import onnx_providers | |
| providers = onnx_providers() | |
| session = ort.InferenceSession(str(checkpoint), sess_options=options, providers=providers) | |
| if session.get_inputs()[0].shape != [1, 3, 256, 256] or session.get_outputs()[0].shape != [1, 3, 256, 256]: | |
| raise ValueError("Unexpected building model input/output schema.") | |
| if providers[0] == "CUDAExecutionProvider": | |
| probe = session.run(None, {session.get_inputs()[0].name: np.zeros((1, 3, 256, 256), dtype="float32")})[0] | |
| if not np.isfinite(probe).all(): | |
| session = ort.InferenceSession(str(checkpoint), sess_options=options, providers=["CPUExecutionProvider"]) | |
| return session | |
| def _sha256(path: str) -> str: | |
| with Path(path).open("rb") as stream: | |
| return hashlib.file_digest(stream, "sha256").hexdigest() | |
| # --------------------------------------------------------------------------- | |
| # Gaussian blend kernel (pre-computed once) | |
| # --------------------------------------------------------------------------- | |
| def _make_kernel(size: int = _TILE, sigma: float = _SIGMA) -> np.ndarray: | |
| axis = np.exp(-0.5 * ((np.arange(size) - (size - 1) / 2.0) / sigma) ** 2).astype("float32") | |
| return np.maximum(np.outer(axis, axis), 1e-8) | |
| _KERNEL = _make_kernel() | |
| # --------------------------------------------------------------------------- | |
| # Instance separation (watershed) | |
| # --------------------------------------------------------------------------- | |
| def _plateau_peaks(distance, foreground, min_distance): | |
| """One candidate per connected maximum plateau, not one per grid interval.""" | |
| surface = np.where(foreground, distance, -np.inf) | |
| maxima = (surface == ndimage.maximum_filter(surface, size=2 * min_distance + 1)) & foreground | |
| # A flat saddle connecting two roofs is not a regional maximum, even if | |
| # its distant higher endpoints lie outside the suppression window. | |
| maxima &= local_maxima(surface, connectivity=1, allow_borders=True) | |
| regions, _ = ndimage.label(maxima) | |
| candidates = [] | |
| for ident, sl in enumerate(ndimage.find_objects(regions), 1): | |
| if sl is None: | |
| continue | |
| points = np.argwhere(regions[sl] == ident) | |
| center = points.mean(axis=0) | |
| point = points[np.argmin(np.sum((points - center) ** 2, axis=1))] | |
| candidates.append(point + [sl[0].start, sl[1].start]) | |
| if not candidates: | |
| return np.empty((0, 2), dtype=int) | |
| candidates = np.asarray(candidates) | |
| order = np.argsort(-distance[tuple(candidates.T)], kind="stable") | |
| kept = [] | |
| buckets = {} | |
| for point in candidates[order]: | |
| cell = tuple(point // min_distance) | |
| nearby = (other for dy in (-1, 0, 1) for dx in (-1, 0, 1) | |
| for other in buckets.get((cell[0] + dy, cell[1] + dx), ())) | |
| if not any(np.max(np.abs(point - other)) < min_distance for other in nearby): | |
| kept.append(point) | |
| buckets.setdefault(cell, []).append(point) | |
| return np.asarray(kept, dtype=int) | |
| def separate_instances(probability, distance, valid, threshold=THRESHOLD): | |
| """Separate overlapping detections into individual building instances. | |
| Returns (compacted_labels, rejected_tiny_count). | |
| Score lookup MUST be computed from the returned compacted labels, not from | |
| watershed output, to guarantee score-ID alignment in the GeoJSON. | |
| """ | |
| foreground = (probability > threshold) & valid | |
| components, n = ndimage.label(foreground) | |
| seeds = np.zeros(foreground.shape, dtype=bool) | |
| sizes = np.bincount(components.ravel()) | |
| large_ids = sizes > _LARGE_BLOB | |
| large_ids[0] = False | |
| large = large_ids[components] | |
| # Collapse flat maxima before applying the existing separation distance. | |
| peaks = _plateau_peaks(distance, foreground & ~large, min_distance=4) | |
| if len(peaks): | |
| seeds[tuple(peaks.T)] = True | |
| # Large blobs retain the wider separation distance and marker cap. | |
| if large.any(): | |
| large_peaks = _plateau_peaks(distance, large, min_distance=8) | |
| if len(large_peaks): | |
| # Cap markers per large component to prevent runaway over-segmentation | |
| for comp_id in np.where(large_ids)[0]: | |
| comp_mask = components == comp_id | |
| comp_peak_mask = np.zeros(foreground.shape, dtype=bool) | |
| comp_peak_mask[tuple(large_peaks.T)] = True | |
| comp_peak_mask &= comp_mask | |
| comp_peak_coords = np.argwhere(comp_peak_mask) | |
| if len(comp_peak_coords) > _MAX_MARKERS: | |
| # Keep the _MAX_MARKERS peaks with highest distance value | |
| dist_vals = distance[comp_peak_mask] | |
| top_idx = np.argsort(dist_vals)[-_MAX_MARKERS:] | |
| comp_peak_coords = comp_peak_coords[top_idx] | |
| seeds[comp_peak_coords[:, 0], comp_peak_coords[:, 1]] = True | |
| # Guarantee at least one marker per connected component (including edge objects) | |
| for i, sl in enumerate(ndimage.find_objects(components), 1): | |
| if sl is None: | |
| continue | |
| region = components[sl] == i | |
| if not np.any(seeds[sl] & region): | |
| y, x = np.unravel_index( | |
| np.argmax(np.where(region, distance[sl], -np.inf)), region.shape | |
| ) | |
| seeds[sl[0].start + y, sl[1].start + x] = True | |
| markers, _ = ndimage.label(seeds) | |
| labels = watershed(-distance, markers, mask=foreground).astype("int32") if n else components | |
| # Remove tiny instances | |
| sizes = np.bincount(labels.ravel()) | |
| rejected = int(np.sum((sizes[1:] < _MIN_PIXELS) & (sizes[1:] > 0))) | |
| keep = sizes >= _MIN_PIXELS | |
| keep[0] = False | |
| labels[~keep[labels]] = 0 | |
| # Compact IDs to 1..N (no gaps) — MUST happen before bincount for scores | |
| ids = np.unique(labels) | |
| ids = ids[ids > 0] | |
| lookup = np.zeros(int(labels.max()) + 1 if labels.max() > 0 else 1, dtype="int32") | |
| lookup[ids] = np.arange(1, len(ids) + 1) | |
| return lookup[labels], rejected | |
| def separate_instances_satlas( | |
| footprint_prob: np.ndarray, | |
| boundary_prob: np.ndarray, | |
| valid: np.ndarray, | |
| foot_thr: float = 0.6, | |
| boundary_thr: float = 0.6, | |
| min_area: int = 48, | |
| min_distance: int = 8, | |
| h_prominence: float = 6.0, | |
| ) -> tuple[np.ndarray, int]: | |
| """Separate building footprints using dual-head footprint and boundary probabilities. | |
| Uses smooth distance transform on foreground building regions combined with | |
| boundary attenuation to prevent over-segmenting single roofs, while cleanly | |
| splitting touching buildings using prominence-guided watershed. | |
| """ | |
| if (footprint_prob.shape != boundary_prob.shape or footprint_prob.shape != valid.shape | |
| or footprint_prob.ndim != 2): | |
| raise ValueError("Building footprint, boundary and validity grids must match.") | |
| if not (0 < foot_thr < 1 and 0 < boundary_thr < 1): | |
| raise ValueError("Building probability thresholds must lie between zero and one.") | |
| if min_area < 1 or min_distance < 1 or h_prominence <= 0: | |
| raise ValueError("Building area, separation distance and prominence must be positive.") | |
| if (not np.isfinite(footprint_prob).all() or not np.isfinite(boundary_prob).all() | |
| or np.any((footprint_prob < 0) | (footprint_prob > 1)) | |
| or np.any((boundary_prob < 0) | (boundary_prob > 1))): | |
| raise ValueError("Building postprocessing requires finite probabilities in 0-1.") | |
| foreground = (footprint_prob > foot_thr) & valid | |
| if not foreground.any(): | |
| return np.zeros(foreground.shape, dtype="int32"), 0 | |
| distance = ndimage.distance_transform_edt(foreground) | |
| effective_dist = distance * (1.0 - 0.4 * boundary_prob) | |
| # Preserve the smooth prominence surface. Thresholding it before finding | |
| # seeds introduces artificial steps and spurious roof centers. | |
| elevation = np.where(boundary_prob >= boundary_thr, -effective_dist, -distance) | |
| # Prominence-based seeds to suppress intra-roof texture noise | |
| candidates = h_maxima(effective_dist, h=h_prominence) & foreground | |
| peaks = _plateau_peaks(effective_dist, candidates, min_distance) | |
| seeds = np.zeros_like(foreground) | |
| if len(peaks): | |
| seeds[tuple(peaks.T)] = True | |
| # Guarantee at least one seed per connected component >= min_area | |
| components, n_comp = ndimage.label(foreground) | |
| if n_comp > 0: | |
| for i, sl in enumerate(ndimage.find_objects(components), 1): | |
| if sl is None: | |
| continue | |
| region = components[sl] == i | |
| if region.sum() < min_area: | |
| continue | |
| if not np.any(seeds[sl] & region): | |
| y, x = np.unravel_index( | |
| np.argmax(np.where(region, distance[sl], -np.inf)), region.shape | |
| ) | |
| seeds[sl[0].start + y, sl[1].start + x] = True | |
| markers, n_markers = ndimage.label(seeds) | |
| labels = watershed(elevation, markers, mask=foreground).astype("int32") if n_markers > 0 else components | |
| else: | |
| labels = np.zeros(foreground.shape, dtype="int32") | |
| sizes = np.bincount(labels.ravel()) | |
| rejected = int(np.sum((sizes[1:] < min_area) & (sizes[1:] > 0))) | |
| keep = sizes >= min_area | |
| keep[0] = False | |
| labels[~keep[labels]] = 0 | |
| ids = np.unique(labels) | |
| ids = ids[ids > 0] | |
| lookup = np.zeros(int(labels.max()) + 1 if labels.max() > 0 else 1, dtype="int32") | |
| lookup[ids] = np.arange(1, len(ids) + 1) | |
| return lookup[labels], rejected | |
| # --------------------------------------------------------------------------- | |
| # Resolution helper | |
| # --------------------------------------------------------------------------- | |
| def resolution_m(src): | |
| if not src.crs: | |
| return None | |
| crs = CRS(src.crs) | |
| if crs.is_projected: | |
| return max(src.res) * crs.axis_info[0].unit_conversion_factor | |
| t = Transformer.from_crs(crs, "EPSG:4326", always_xy=True) | |
| cx, cy = src.width / 2, src.height / 2 | |
| x = src.transform.a * cx + src.transform.b * cy + src.transform.c | |
| y = src.transform.d * cx + src.transform.e * cy + src.transform.f | |
| lon, lat = t.transform(x, y) | |
| distances = [] | |
| for dx,dy in ((src.transform.a,src.transform.d),(src.transform.b,src.transform.e)): | |
| lon2,lat2 = t.transform(x+dx,y+dy) | |
| distances.append(abs(Geod(ellps="WGS84").inv(lon,lat,lon2,lat2)[2])) | |
| return max(distances) | |
| # --------------------------------------------------------------------------- | |
| # Main detection entry point | |
| # --------------------------------------------------------------------------- | |
| def detect_buildings(path, output, progress=None, *, water_result=None): | |
| from satquery_engine.services.radiometry import rgb_unit_data, rgb_indexes | |
| with rasterio.open(path) as src: | |
| if src.count < 3: | |
| raise ValueError("Building analysis needs optical RGB bands.") | |
| if src.width * src.height > 32_000_000: | |
| raise ValueError( | |
| "Building analysis currently supports up to 32 million pixels per image." | |
| ) | |
| indexes = rgb_indexes(src) | |
| gsd = resolution_m(src) | |
| if gsd is not None and gsd > 1.5: | |
| raise ValueError( | |
| f"This image has approximately {gsd:.1f}-metre pixels, " | |
| "too coarse for a reliable individual-building count with this model. " | |
| "Individual building footprint segmentation requires very high resolution satellite or aerial imagery (GSD <= 1.5m). " | |
| "For satellite imagery at this resolution, query 'Identify built-up areas' to map urban surface footprints." | |
| ) | |
| # Image quality check: blur and contrast | |
| sample_h = min(src.height, 512) | |
| sample_w = min(src.width, 512) | |
| sample = src.read(indexes[:3], out_shape=(3, sample_h, sample_w), resampling=rasterio.enums.Resampling.bilinear).astype("float32") | |
| sample_gray = np.mean(sample, axis=0) | |
| s_span = float(np.ptp(sample_gray)) if sample_gray.size else 1.0 | |
| s_norm = (sample_gray - np.min(sample_gray)) / max(s_span, 1e-6) | |
| blur_var = float(ndimage.laplace(s_norm).var()) | |
| contrast_val = float(s_norm.std()) | |
| if gsd is not None: | |
| if gsd <= 0.8 and blur_var >= 0.003 and contrast_val >= 0.10: | |
| count_reliability = "COUNT_RELIABLE" | |
| elif gsd <= 1.5: | |
| count_reliability = "COUNT_APPROXIMATE" | |
| else: | |
| count_reliability = "COUNT_UNRELIABLE" | |
| else: | |
| count_reliability = "COUNT_APPROXIMATE" if (blur_var >= 0.003 and contrast_val >= 0.10) else "COUNT_UNRELIABLE" | |
| # Load FP32 checkpoint or fine-tuned model bundle | |
| fallback_events = [] | |
| fp32_path = Path(settings.building_checkpoint) | |
| try: | |
| fp32_path = _resolve_checkpoint() | |
| session_or_bundle = load_session(str(fp32_path)) | |
| except (ValueError, RuntimeError, OSError) as primary_error: | |
| # An explicit selection must never silently run another model. | |
| if os.getenv("SATQUERY_BUILDING_CHECKPOINT") or fp32_path != Path(settings.building_checkpoint): | |
| raise | |
| from satquery_engine.models.registry import LocalModelRegistry | |
| fallback = LocalModelRegistry(settings.model_dir).get("building_fallback") | |
| if fallback is None or not fallback.availability: | |
| raise ValueError("The primary building model failed and no verified DINO fallback is available.") from primary_error | |
| fp32_path = fallback.local_path / "onnx/model.onnx" | |
| session_or_bundle = load_session(str(fp32_path)) | |
| fallback_events.extend(["PRIMARY_MODEL_FAILED", "FALLBACK_USED", | |
| f"Primary building model failed: {primary_error}"]) | |
| from satquery_engine.services.buildings_rf import RfBuildingSession | |
| is_rf = isinstance(session_or_bundle, RfBuildingSession) | |
| is_bundle = isinstance(session_or_bundle, tuple) | |
| h, w = src.height, src.width | |
| scene_rgb, scene_valid, _ = rgb_unit_data(src) | |
| del scene_rgb | |
| if is_rf: | |
| from satquery_engine.services.buildings_rf import ( | |
| MODEL_ID as RF_MODEL_ID, TILE as RF_TILE, STRIDE as RF_STRIDE, | |
| THRESHOLD as RF_THRESHOLD, detect_rf_instances, | |
| ) | |
| labels, prediction_prob, duplicate_count, tile_records, radiometry = detect_rf_instances( | |
| path, session_or_bundle, scene_valid, progress=progress | |
| ) | |
| radiometry = {"source_rgb": radiometry, | |
| "model_processor": "RfDetrImageProcessor: resize to 432x432, divide by 255, ImageNet normalization"} | |
| valid = scene_valid | |
| positions = tile_records | |
| rejected = 0 | |
| active_tile_size = RF_TILE | |
| active_stride = RF_STRIDE | |
| active_model_id = RF_MODEL_ID | |
| active_method = "RF-DETR satellite building instance masks with global duplicate suppression" | |
| active_execution = "rf_detr_satellite_safetensors" | |
| from satquery_engine.models.device import torch_device | |
| active_device = torch_device() | |
| active_checkpoint_sha = session_or_bundle.checkpoint_sha256 | |
| active_threshold = RF_THRESHOLD | |
| from satquery_engine.services.buildings_rf import MIN_INSTANCE_PIXELS | |
| active_min_pixels = MIN_INSTANCE_PIXELS | |
| active_benchmark_f1 = None | |
| postprocess_audit = {"method": "score_ordered_instance_merge", "score_threshold": RF_THRESHOLD, | |
| "overlap_iou_threshold": .5, "containment_threshold": .8, | |
| "minimum_instance_pixels": active_min_pixels, | |
| "minimum_unclaimed_fraction": .5, | |
| "island_minimum_fraction": .02, | |
| "tile_edge_policy": "prefer complete interior predictions over internally clipped masks", | |
| "probability_kind": "sparse_instance_score"} | |
| elif is_bundle: | |
| import torch | |
| from satquery_engine.models.device import torch_device | |
| model, config, bundle_metrics = session_or_bundle | |
| device = torch_device() if hasattr(model, "to") else "cpu" | |
| if hasattr(model, "to"): | |
| model.to(device) | |
| tile_size = int(config.get("input", {}).get("tile_size", 512)) | |
| overlap = int(config.get("tile_inference", {}).get("overlap", 128)) | |
| stride = tile_size - overlap | |
| foot_thr = float(config.get("postprocess", {}).get("foot_thr", 0.6)) | |
| bound_thr = float(config.get("postprocess", {}).get("boundary_thr", 0.6)) | |
| min_area = int(config.get("postprocess", {}).get("min_area", 48)) | |
| min_dist = int(config.get("postprocess", {}).get("min_distance", 8)) | |
| h_prom = float(config.get("postprocess", {}).get("h_prominence", 6.0)) | |
| if tile_size < 1 or not 0 <= overlap < tile_size: | |
| raise ValueError("Invalid building tile size or overlap.") | |
| if not (0 < foot_thr < 1 and 0 < bound_thr < 1): | |
| raise ValueError("Invalid building model thresholds.") | |
| hann_1d = np.hanning(tile_size) | |
| hann_2d = np.outer(hann_1d, hann_1d).astype("float32") | |
| hann_2d = np.maximum(hann_2d, 1e-4) | |
| accum = np.zeros((2, h, w), dtype="float32") | |
| weight = np.zeros((h, w), dtype="float32") | |
| valid = np.zeros((h, w), dtype=bool) | |
| positions = [(y, x) for y in _starts(h, tile_size, stride) for x in _starts(w, tile_size, stride)] | |
| tile_records = [] | |
| for index, (y, x) in enumerate(positions): | |
| th, tw = min(tile_size, h - y), min(tile_size, w - x) | |
| raw, good, radiometry = rgb_unit_data(src, window=Window(x, y, tw, th), scene_valid=scene_valid) | |
| patch = np.pad(raw, ((0, 0), (0, tile_size-th), (0, tile_size-tw)), mode="edge") | |
| t_patch = torch.from_numpy(patch).unsqueeze(0).to(device) | |
| with torch.inference_mode(): | |
| logits = model(t_patch)[0].cpu().numpy() | |
| if logits.shape != (2, tile_size, tile_size) or not np.isfinite(logits).all(): | |
| raise ValueError("Building model returned invalid logits.") | |
| probs_patch = expit(logits) | |
| k = hann_2d[:th, :tw] | |
| accum[:, y:y + th, x:x + tw] += probs_patch[:, :th, :tw] * k | |
| weight[y:y + th, x:x + tw] += k | |
| valid[y:y + th, x:x + tw] |= good | |
| tile_records.append({ | |
| "tile_id": index, "x_offset": x, "y_offset": y, "width": tw, "height": th, | |
| "origin": [x, y], "shape": [th, tw], | |
| "padding": [tile_size - th, tile_size - tw], "scale": 1, | |
| "transform": list(src.window_transform(Window(x, y, tw, th)))[:6], | |
| "valid_pixels": int(good.sum()), | |
| }) | |
| if progress: | |
| progress(f"Detecting buildings: tile {index + 1} of {len(positions)}") | |
| if device == "cuda" and hasattr(model, "cpu"): | |
| model.cpu() | |
| torch.cuda.empty_cache() | |
| if not valid.any() or np.any(weight <= 0): | |
| raise ValueError("Building inference lacks complete valid coverage.") | |
| safe_weight = np.where(weight > 0, weight, 1.0) | |
| probs = accum / safe_weight | |
| footprint_prob = probs[0] | |
| boundary_prob = probs[1] | |
| labels, rejected = separate_instances_satlas( | |
| footprint_prob, boundary_prob, valid, | |
| foot_thr=foot_thr, boundary_thr=bound_thr, | |
| min_area=min_area, min_distance=min_dist, | |
| h_prominence=h_prom, | |
| ) | |
| prediction_prob = footprint_prob | |
| active_tile_size = tile_size | |
| active_stride = stride | |
| active_model_id = "satquery_buildings_bundle" | |
| active_method = "SatlasBuildingNet Aerial Swin-v2-Base with FPN, Hann blend tiling, and boundary-guided watershed" | |
| active_execution = "satquery_buildings_bundle" | |
| active_device = device | |
| active_checkpoint_sha = _sha256(str(Path(fp32_path) / "satquery_buildings_state_dict.pt")) | |
| active_threshold = foot_thr | |
| active_min_pixels = min_area | |
| active_benchmark_f1 = bundle_metrics.get("test", {}).get("instance_f1@0.5") | |
| duplicate_count = 0 | |
| postprocess_audit = { | |
| "method": "boundary_guided_watershed_v2", | |
| "footprint_threshold": foot_thr, "boundary_threshold": bound_thr, | |
| "minimum_area_pixels": min_area, "minimum_seed_distance_pixels": min_dist, | |
| "h_prominence": h_prom, | |
| "threshold_source": config.get("threshold_source", "bundle configuration"), | |
| } | |
| else: | |
| from satquery_engine.models.device import torch_device | |
| session = session_or_bundle | |
| accum = np.zeros((3, h, w), dtype="float32") | |
| weight = np.zeros((h, w), dtype="float32") | |
| valid = np.zeros((h, w), dtype=bool) | |
| positions = [(y, x) for y in _starts(h, _TILE, _STRIDE) for x in _starts(w, _TILE, _STRIDE)] | |
| tile_records = [] | |
| for index, (y, x) in enumerate(positions): | |
| th, tw = min(_TILE, h - y), min(_TILE, w - x) | |
| raw, good, radiometry = rgb_unit_data(src, window=Window(x, y, tw, th), scene_valid=scene_valid) | |
| tile = np.zeros((3, _TILE, _TILE), dtype="float32") | |
| tile[:, :th, :tw] = raw | |
| logits = session.run(None, {"image": ((tile - MEAN) / STD)[None].astype("float32")})[0][0] | |
| if logits.shape != (3, _TILE, _TILE) or not np.isfinite(logits).all(): | |
| raise ValueError("Building model returned invalid logits.") | |
| predictions = np.stack([expit(logits[0]), expit(logits[1]), np.tanh(logits[2])]) | |
| k = _KERNEL[:th, :tw] | |
| accum[:, y:y + th, x:x + tw] += predictions[:, :th, :tw] * k | |
| weight[y:y + th, x:x + tw] += k | |
| valid[y:y + th, x:x + tw] |= good | |
| tile_records.append({ | |
| "tile_id": index, "x_offset": x, "y_offset": y, "width": tw, "height": th, | |
| "origin": [x, y], "shape": [th, tw], | |
| "padding": [_TILE - th, _TILE - tw], "scale": 1, | |
| "transform": list(src.window_transform(Window(x, y, tw, th)))[:6], | |
| "valid_pixels": int(good.sum()), | |
| }) | |
| if progress: | |
| progress(f"Detecting buildings: tile {index + 1} of {len(positions)}") | |
| # Blend: only require coverage where pixels are valid | |
| uncovered_valid = valid & (weight == 0) | |
| if uncovered_valid.any(): | |
| raise ValueError("Building inference did not cover every valid pixel; no count generated.") | |
| safe_weight = np.where(weight > 0, weight, 1.0) | |
| predictions = accum / safe_weight | |
| if not valid.any(): | |
| raise ValueError("No valid RGB pixels are available for building analysis.") | |
| labels, rejected = separate_instances(predictions[0], predictions[2], valid) | |
| prediction_prob = predictions[0] | |
| active_tile_size = _TILE | |
| active_stride = _STRIDE | |
| active_model_id = MODEL_ID | |
| active_method = "DINOv3 footprint/boundary/distance heads with global watershed and overlapping probability merge" | |
| active_execution = "original_released_onnx" | |
| active_device = session.get_providers()[0] if hasattr(session, "get_providers") else "unknown" | |
| if active_device == "CPUExecutionProvider" and torch_device() == "cuda": | |
| fallback_events.append("GPU_PROVIDER_FAILED: DINO ONNX CUDA output was nonfinite; verified CPU inference was used.") | |
| active_checkpoint_sha = _sha256(str(fp32_path)) | |
| active_threshold = THRESHOLD | |
| active_min_pixels = _MIN_PIXELS | |
| active_benchmark_f1 = None | |
| duplicate_count = 0 | |
| postprocess_audit = {"method": "dinov3_distance_watershed", "footprint_threshold": THRESHOLD, | |
| "minimum_area_pixels": _MIN_PIXELS} | |
| if not valid.any(): | |
| raise ValueError("No valid RGB pixels are available for building analysis.") | |
| surface_exclusion = {"status": "not_requested"} | |
| if water_result is not None: | |
| from satquery_engine.services.surface_context import exclude_water_instances | |
| labels, surface_exclusion = exclude_water_instances( | |
| labels, valid, path, water_result, active_min_pixels, | |
| ) | |
| # Compute scores from COMPACTED labels (guaranteeing score keys == GeoJSON IDs) | |
| sums = np.bincount(labels.ravel(), weights=prediction_prob.ravel()) | |
| sizes = np.bincount(labels.ravel()) | |
| scores = { | |
| i: float(sums[i] / sizes[i]) | |
| for i in range(1, len(sizes)) | |
| if sizes[i] > 0 | |
| } | |
| # Export with numbered overlay for visual count verification | |
| result = export_labels( | |
| labels, path, output, "buildings", | |
| scores=scores, valid_mask=valid, numbered=True | |
| ) | |
| if len(result["features"]) != len(scores): | |
| raise ValueError("Building count does not match final features; result withheld.") | |
| # Persist probability map | |
| np.save(output / "building_probability.npy", prediction_prob) | |
| result["paths"].append(output / "building_probability.npy") | |
| from satquery_engine.services.spatial_outputs import export_float_raster | |
| result["paths"].append( | |
| export_float_raster( | |
| np.where(valid, prediction_prob, np.nan), path, | |
| output / "buildings_probability.tif" | |
| ) | |
| ) | |
| if is_bundle: | |
| result["paths"].append(export_float_raster( | |
| np.where(valid, boundary_prob, np.nan), path, | |
| output / "buildings_boundary_probability.tif", | |
| )) | |
| # Export binary buildings_mask.tif | |
| from rasterio.transform import Affine | |
| mask_tif = output / "buildings_mask.tif" | |
| with rasterio.open(path) as src_meta: | |
| grid = src_meta.transform @ Affine.scale(src_meta.width / labels.shape[1], src_meta.height / labels.shape[0]) | |
| prof = dict( | |
| driver="GTiff", | |
| width=labels.shape[1], | |
| height=labels.shape[0], | |
| count=1, | |
| dtype="uint8", | |
| crs=src_meta.crs, | |
| transform=grid, | |
| compress="deflate", | |
| ) | |
| with rasterio.open(mask_tif, "w", **prof) as dst: | |
| dst.write((labels > 0).astype("uint8") * 255, 1) | |
| dst.write_mask(valid.astype("uint8") * 255) | |
| result["paths"].append(mask_tif) | |
| tile_manifest = output / "building_tiles.json" | |
| tile_manifest.write_text(json.dumps({ | |
| "tile_size": active_tile_size, "overlap": active_tile_size - active_stride, | |
| "stride": active_stride, "gaussian_sigma": _SIGMA if not is_bundle and not is_rf else None, | |
| "tiles": tile_records, "radiometry": radiometry, | |
| "checkpoint_sha256": active_checkpoint_sha, | |
| "checkpoint_execution": active_execution, | |
| "device": active_device, | |
| "postprocessing": postprocess_audit, | |
| }, indent=2)) | |
| result["paths"].append(tile_manifest) | |
| limitations = [ | |
| "Predicted visible footprints are an estimate, not a surveyed building inventory. " | |
| "Touching roofs, shadows and tiny buildings can cause count errors.", | |
| "Model scores are uncalibrated. Benchmark performance may not transfer to this scene.", | |
| ] | |
| limitations.extend(fallback_events) | |
| if surface_exclusion["status"] == "applied": | |
| limitations.append("Building footprints exclude the canonical estimated water mask. Water classification errors can affect shoreline buildings.") | |
| if is_rf: | |
| limitations.append("Selected on satellite imagery. This checkpoint performed poorly on the local aerial reference chips; aerial counts require separate validation.") | |
| benchmark_path = Path(__file__).resolve().parents[3] / "artifacts/building_benchmark/benchmark.json" | |
| benchmark = None | |
| if not is_bundle and not is_rf and benchmark_path.exists(): | |
| benchmark = json.loads(benchmark_path.read_text()) | |
| if (benchmark.get("f1") or 0) < .6: | |
| limitations.append( | |
| f'The measured {benchmark["sample_count"]}-image validation subset had low instance ' | |
| f'accuracy (F1 {benchmark["f1"]:.3f}); treat this count as a low-confidence estimate.' | |
| ) | |
| if gsd is None: | |
| limitations.append( | |
| "Ground resolution is unknown; geographic areas and exact count reliability cannot be established." | |
| ) | |
| elif gsd > 1: | |
| limitations.append("Pixels exceed one metre; small and touching buildings may be missed.") | |
| if count_reliability == "COUNT_UNRELIABLE": | |
| limitations.append( | |
| "I can detect built-up areas, but this image is not detailed enough for an accurate individual-building count." | |
| ) | |
| if is_rf: | |
| count_reliability = "COUNT_UNRELIABLE" | |
| limitations.append( | |
| "The satellite-model selection set performed well, but a separate aerial reference set had no matched footprints; " | |
| "accuracy on this scene's sensor and geography is unvalidated." | |
| ) | |
| elif active_benchmark_f1 is not None and active_benchmark_f1 < .6: | |
| count_reliability = "COUNT_UNRELIABLE" | |
| limitations.append( | |
| f"The fine-tuned model's held-out test instance F1 was {active_benchmark_f1:.3f}; " | |
| "the footprint count is reported as low-confidence." | |
| ) | |
| elif not is_bundle and (benchmark is None or (benchmark.get("f1") or 0) < .6): | |
| count_reliability = "COUNT_UNRELIABLE" | |
| elif count_reliability == "COUNT_RELIABLE": | |
| count_reliability = "COUNT_APPROXIMATE" | |
| small_objs = int(np.sum((sizes[1:] > 0) & (sizes[1:] < 20))) | |
| uncertain_objs = int(sum(1 for s in scores.values() if s < 0.50)) | |
| mean_conf = float(np.mean(list(scores.values()))) if scores else 0.0 | |
| result.update( | |
| count=len(result["features"]), | |
| surface_exclusion=surface_exclusion, | |
| building_count=len(result["features"]), | |
| building_count_type="visible_footprint_count", | |
| confidence="HIGH" if count_reliability == "COUNT_RELIABLE" else ("MEDIUM" if count_reliability == "COUNT_APPROXIMATE" else "LOW"), | |
| count_reliability=count_reliability, | |
| average_confidence=round(mean_conf, 4), | |
| detected_area=result.get("area_m2"), | |
| rejected_candidates=rejected, | |
| duplicate_candidates_removed=duplicate_count, | |
| small_objects=small_objs, | |
| uncertain_objects=uncertain_objs, | |
| count_state="ESTIMATED", | |
| rejected_tiny=rejected, | |
| rejected_duplicates=duplicate_count, | |
| checkpoint_sha256=active_checkpoint_sha, | |
| checkpoint_execution=active_execution, | |
| device=active_device, | |
| preprocessing=radiometry, | |
| postprocessing=postprocess_audit, | |
| small_object_rate=float(np.mean(sizes[1:] < 20)) if len(sizes) > 1 else 0., | |
| count_distribution={ | |
| "under_20_pixels": small_objs, | |
| "at_least_20_pixels": int(np.sum(sizes[1:] >= 20)), | |
| }, | |
| duplicate_policy=("score-ordered global instance overlap suppression" if is_rf else | |
| "overlapping probability maps blended before global instance extraction; no per-tile instances counted"), | |
| tile_count=len(positions), | |
| tile_overlap_px=active_tile_size - active_stride, | |
| complete_coverage=True, | |
| resolution_m=gsd, | |
| blur_variance=blur_var, | |
| contrast_std=contrast_val, | |
| minimum_instance_pixels=active_min_pixels, | |
| mean_model_score=mean_conf if scores else None, | |
| evidence_strength=(mean_conf if scores else 0.0 if is_rf else | |
| float(np.mean(1 - prediction_prob[valid])) if valid.any() else 0.0), | |
| model_id=active_model_id, | |
| threshold=active_threshold, | |
| benchmark_f1=active_benchmark_f1 if is_bundle else (benchmark.get("f1") if benchmark else None), | |
| method=active_method, | |
| limitations=limitations, | |
| fallback_events=fallback_events, | |
| ) | |
| q_pass, q_warnings = check_building_result(result, labels, scores, gsd) | |
| result["quality_gate_passed"] = q_pass | |
| result["quality_warnings"] = q_warnings | |
| result["evidence_strength"] = min(result["evidence_strength"], .25 if count_reliability == "COUNT_UNRELIABLE" else .5) | |
| result["domain_compatibility"] = "unverified" if gsd is None else "resolution_checked_domain_unvalidated" | |
| stats = output / "buildings_stats.json" | |
| stats.write_text(json.dumps({k:v for k,v in result.items() if k not in {"paths", "features"}}, indent=2, allow_nan=False)) | |
| result["paths"].append(stats) | |
| return result | |
| def check_building_result( | |
| result: dict[str, Any], | |
| labels: np.ndarray, | |
| scores: dict[int, float], | |
| gsd: float | None, | |
| ) -> tuple[bool, list[str]]: | |
| """Context-aware sanity check on building detection results.""" | |
| warnings: list[str] = [] | |
| count = len(scores) | |
| if count == 0: | |
| return True, ["No building footprints detected in this scene."] | |
| small_objects = result.get("small_objects", 0) | |
| if float(small_objects) / max(1, count) > 0.60: | |
| warnings.append("Over 60% of detected structures are small (under 20 pixels); possible over-segmentation or noise.") | |
| uncertain_objects = result.get("uncertain_objects", 0) | |
| if float(uncertain_objects) / max(1, count) > 0.50: | |
| warnings.append("Over 50% of detected buildings have model confidence below 0.50.") | |
| if gsd is not None and gsd > 1.5: | |
| warnings.append(f"Image resolution ({gsd:.1f}m GSD) is coarse; individual structures may be merged or omitted.") | |
| if result.get("count_reliability") == "COUNT_UNRELIABLE": | |
| warnings.append("I can detect built-up areas, but this image is not detailed enough for an accurate individual-building count.") | |
| return len(warnings) == 0, warnings | |
| # --------------------------------------------------------------------------- | |
| # Building change matching | |
| # --------------------------------------------------------------------------- | |
| def match_buildings(a, b, output): | |
| from shapely.geometry import shape | |
| from scipy.optimize import linear_sum_assignment | |
| left = [shape(f["properties"]["pixel_geometry"]) for f in a["features"]] | |
| right = [shape(f["properties"]["pixel_geometry"]) for f in b["features"]] | |
| scores = np.zeros((len(left), len(right))) | |
| for i, x in enumerate(left): | |
| for j, y in enumerate(right): | |
| if x.intersects(y): | |
| scores[i, j] = x.intersection(y).area / x.union(y).area | |
| matches = [] | |
| if scores.size: | |
| rows, cols = linear_sum_assignment(-scores) | |
| matches = [(int(i), int(j)) for i, j in zip(rows, cols) if scores[i, j] >= .3] | |
| ma = {i for i, j in matches} | |
| mb = {j for i, j in matches} | |
| features = [] | |
| for role, result, matched in [("REMOVED", a, ma), ("NEW", b, mb)]: | |
| for i, f in enumerate(result["features"]): | |
| features.append({ | |
| **f, | |
| "properties": { | |
| **f["properties"], | |
| "status": "UNCHANGED" if i in matched else "POSSIBLE_CHANGE", | |
| "candidate_type": None if i in matched else role, | |
| "image": 1 if role == "REMOVED" else 2, | |
| }, | |
| }) | |
| geojson = output / "building_change.geojson" | |
| geojson.write_text(json.dumps({ | |
| "type": "FeatureCollection", "features": features, | |
| "properties": {"crs": "EPSG:4326" if a["area_m2"] is not None else None}, | |
| }, allow_nan=False)) | |
| from PIL import Image, ImageDraw | |
| preview = next(p for p in b["paths"] if p.name == "buildings_original.png") | |
| labels_path = next(p for p in b["paths"] if p.name == "buildings_labels.tif") | |
| with rasterio.open(labels_path) as grid: | |
| width, height = grid.width, grid.height | |
| base = Image.open(preview).convert("RGBA") | |
| marks = Image.new("RGBA", base.size) | |
| draw = ImageDraw.Draw(marks) | |
| for f in features: | |
| props = f["properties"] | |
| if props["image"] == 1 and props["status"] == "UNCHANGED": | |
| continue | |
| color = ( | |
| (0, 190, 220, 100) if props["status"] == "UNCHANGED" | |
| else (245, 165, 30, 160) if props["candidate_type"] == "NEW" | |
| else (225, 65, 100, 160) | |
| ) | |
| geometry = props["pixel_geometry"] | |
| rings = ( | |
| geometry["coordinates"] if geometry["type"] == "Polygon" | |
| else [r for poly in geometry["coordinates"] for r in poly] | |
| ) | |
| for ring in rings: | |
| draw.polygon( | |
| [(x * base.width / width, y * base.height / height) for x, y in ring], | |
| fill=color, outline="white", | |
| ) | |
| overlay = output / "building_change_overlay.png" | |
| Image.alpha_composite(base, marks).convert("RGB").save(overlay) | |
| return { | |
| "before_count": a["count"], | |
| "after_count": b["count"], | |
| "net_count_change": b["count"] - a["count"], | |
| "persistent_count": len(matches), | |
| "possible_new_count": b["count"] - len(mb), | |
| "possible_removed_count": a["count"] - len(ma), | |
| "net_footprint_area_m2": ( | |
| b["area_m2"] - a["area_m2"] | |
| if a["area_m2"] is not None and b["area_m2"] is not None | |
| else None | |
| ), | |
| "paths": [geojson, overlay], | |
| "method": "One-to-one footprint IoU matching on the registered grid", | |
| "evidence_strength": min(a["evidence_strength"], b["evidence_strength"]), | |
| "limitations": [ | |
| "Unmatched footprints are possible changes; occlusion, image quality and alignment " | |
| "can also cause unmatched detections.", | |
| "Change overlay: cyan = matched, amber = possibly new, pink = possibly removed. " | |
| "Colors do not confirm real-world construction or demolition.", | |
| ], | |
| } | |