Spaces:
Paused
Paused
| from __future__ import annotations | |
| import json | |
| from pathlib import Path | |
| from typing import Any | |
| import numpy as np | |
| import rasterio | |
| from PIL import Image | |
| from pyproj import Transformer | |
| from rasterio.enums import Resampling | |
| from rasterio.features import shapes | |
| from rasterio.transform import Affine | |
| from shapely.geometry import shape | |
| from shapely.ops import transform as shapely_transform | |
| from app.schemas import QualityReport, RasterMetadata | |
| from app.services.spectral import available_indices | |
| def _normalize(array: np.ndarray) -> np.ndarray: | |
| arr = array.astype("float32") | |
| finite = arr[np.isfinite(arr)] | |
| if finite.size == 0: | |
| return np.zeros_like(arr, dtype="float32") | |
| low, high = np.percentile(finite, [2, 98]) | |
| if high <= low: | |
| return np.zeros_like(arr, dtype="float32") | |
| return np.clip((arr - low) / (high - low), 0, 1) | |
| def inspect_raster(path: Path) -> RasterMetadata: | |
| with rasterio.open(path) as src: | |
| sample_h = min(src.height, 512) | |
| sample_w = min(src.width, 512) | |
| mask = src.read_masks(1, out_shape=(sample_h, sample_w), resampling=Resampling.nearest) | |
| nodata_percent = float((mask == 0).mean() * 100) | |
| descriptions = list(src.descriptions or ()) | |
| band_names = [descriptions[i] or f"band_{i + 1}" for i in range(src.count)] | |
| # Sensor modality and format inspection | |
| from app.services.croma_pipeline import inspect_modality | |
| mod_report = inspect_modality(path) | |
| return RasterMetadata( | |
| filename=path.name, | |
| width=src.width, | |
| height=src.height, | |
| bands=src.count, | |
| dtype=str(src.dtypes[0]), | |
| crs=src.crs.to_string() if src.crs else None, | |
| bounds=[float(v) for v in src.bounds], | |
| resolution=[abs(float(src.res[0])), abs(float(src.res[1]))], | |
| nodata=float(src.nodata) if src.nodata is not None else None, | |
| nodata_percent=round(nodata_percent, 3), | |
| band_names=band_names, | |
| available_indices=available_indices(path), | |
| source_format=mod_report.file_format, | |
| file_format=mod_report.file_format, | |
| modality=mod_report.modality, | |
| sensor_verified=mod_report.sensor_verified, | |
| georeferenced=mod_report.georeferenced, | |
| metadata_available=mod_report.metadata_available, | |
| analysis_capabilities=mod_report.analysis_capabilities, | |
| limitations=mod_report.limitations, | |
| ) | |
| def inspect_image_quality(path: Path, max_size: int = 1024) -> dict: | |
| """ | |
| Stage 1 quality diagnostics: compute brightness, dynamic range, | |
| potential cloud/shadow coverage, and saturation statistics. | |
| Returns a dict with quality metrics that feed the arbiter. | |
| """ | |
| try: | |
| with rasterio.open(path) as src: | |
| ratio = min(1.0, max_size / max(src.width, src.height)) | |
| out_h = max(1, round(src.height * ratio)) | |
| out_w = max(1, round(src.width * ratio)) | |
| # Read first 3 bands (or fewer) | |
| n_bands = min(src.count, 3) | |
| bands = [] | |
| for i in range(1, n_bands + 1): | |
| band = src.read(i, out_shape=(out_h, out_w), resampling=Resampling.bilinear).astype("float32") | |
| bands.append(band) | |
| if not bands: | |
| return {"quality_score": 0.5, "diagnostics": "No readable bands"} | |
| stack = np.stack(bands, axis=0) | |
| # Normalize to [0, 1] | |
| max_val = float(stack.max()) | |
| if max_val > 1.5: | |
| stack = stack / (255.0 if max_val <= 255.0 else max_val) | |
| stack = np.clip(stack, 0.0, 1.0) | |
| # Compute luminance (mean across channels) | |
| lum = stack.mean(axis=0) | |
| mean_brightness = float(lum.mean()) | |
| std_brightness = float(lum.std()) | |
| # Dynamic range (2nd to 98th percentile) | |
| p2, p98 = np.percentile(lum, [2, 98]) | |
| dynamic_range = float(p98 - p2) | |
| # Potential cloud pixels (very bright, low contrast) | |
| cloud_suspect = float((lum > 0.90).mean() * 100) | |
| # Potential shadow pixels (very dark) | |
| shadow_suspect = float((lum < 0.05).mean() * 100) | |
| # Saturated pixels (any channel at max) | |
| saturated = float((stack.max(axis=0) > 0.99).mean() * 100) | |
| # Quality score: penalize clouds, shadows, low dynamic range, saturation | |
| quality_score = 1.0 | |
| quality_score -= min(0.3, cloud_suspect / 100 * 0.5) | |
| quality_score -= min(0.2, shadow_suspect / 100 * 0.4) | |
| quality_score -= max(0, 0.3 - dynamic_range) * 0.5 | |
| quality_score -= min(0.1, saturated / 100 * 0.2) | |
| quality_score = round(max(0.1, min(1.0, quality_score)), 3) | |
| return { | |
| "quality_score": quality_score, | |
| "mean_brightness": round(mean_brightness, 4), | |
| "std_brightness": round(std_brightness, 4), | |
| "dynamic_range": round(dynamic_range, 4), | |
| "cloud_suspect_percent": round(cloud_suspect, 2), | |
| "shadow_suspect_percent": round(shadow_suspect, 2), | |
| "saturated_percent": round(saturated, 2), | |
| "image_size": [out_w, out_h], | |
| "n_bands": n_bands, | |
| } | |
| except Exception: | |
| return {"quality_score": 0.5, "diagnostics": "Could not compute quality metrics"} | |
| def validate_inputs(paths: list[Path]) -> tuple[list[RasterMetadata], QualityReport]: | |
| metadata = [inspect_raster(path) for path in paths] | |
| blockers: list[str] = [] | |
| warnings: list[str] = [] | |
| checks: dict[str, Any] = {"file_count": len(paths)} | |
| has_any_crs = any(item.crs is not None for item in metadata) | |
| all_have_crs = all(item.crs is not None for item in metadata) | |
| checks["geospatial"] = all_have_crs | |
| for item in metadata: | |
| if item.crs is None: | |
| warnings.append( | |
| f"{item.filename}: missing CRS — pixel-space analysis only; " | |
| "area/distance values will not be available" | |
| ) | |
| if item.nodata_percent > 40: | |
| warnings.append(f"{item.filename}: {item.nodata_percent:.1f}% NoData") | |
| if item.bands < 1: | |
| blockers.append(f"{item.filename}: no readable bands") | |
| if not item.available_indices: | |
| warnings.append( | |
| f"{item.filename}: no deterministic spectral index is available; name Red/Green/NIR/SWIR bands for NDVI/NDWI/NDBI" | |
| ) | |
| compatible = not blockers | |
| if len(metadata) == 2: | |
| a, b = metadata | |
| same_shape = (a.width, a.height) == (b.width, b.height) | |
| if all_have_crs: | |
| # Full geospatial validation when both images have CRS | |
| same_crs = a.crs == b.crs | |
| same_bounds = bool(np.allclose(a.bounds, b.bounds, rtol=0, atol=max(a.resolution + b.resolution))) | |
| resolution_ratio = max(a.resolution[0], b.resolution[0]) / max(min(a.resolution[0], b.resolution[0]), 1e-9) | |
| intersection_width = max(0.0, min(a.bounds[2], b.bounds[2]) - max(a.bounds[0], b.bounds[0])) | |
| intersection_height = max(0.0, min(a.bounds[3], b.bounds[3]) - max(a.bounds[1], b.bounds[1])) | |
| intersection_area = intersection_width * intersection_height | |
| a_area = max(1e-9, (a.bounds[2] - a.bounds[0]) * (a.bounds[3] - a.bounds[1])) | |
| b_area = max(1e-9, (b.bounds[2] - b.bounds[0]) * (b.bounds[3] - b.bounds[1])) | |
| overlap_ratio = intersection_area / min(a_area, b_area) | |
| alignment_score = 1.0 if same_crs and same_shape and same_bounds else max(0.0, overlap_ratio / resolution_ratio) | |
| checks.update({ | |
| "same_crs": same_crs, | |
| "same_shape": same_shape, | |
| "same_bounds": same_bounds, | |
| "resolution_ratio": round(resolution_ratio, 4), | |
| "geographic_overlap_ratio": round(overlap_ratio, 4), | |
| "alignment_score": round(alignment_score, 4), | |
| }) | |
| if not same_crs: | |
| warnings.append("Paired rasters use different CRS values — will reproject to primary raster CRS") | |
| if not same_shape: | |
| warnings.append(f"Paired rasters have different pixel dimensions ({a.width}×{a.height} vs {b.width}×{b.height}); auto-normalized") | |
| if not same_bounds and overlap_ratio < 0.1: | |
| blockers.append("Paired rasters do not share sufficient geographic overlap (< 10%)") | |
| elif not same_bounds: | |
| warnings.append(f"Paired rasters have partial geographic overlap ({overlap_ratio * 100:.1f}%)") | |
| elif has_any_crs: | |
| # A single georeferenced image cannot establish where the other | |
| # image belongs. Reject paired measurements instead of inventing | |
| # pixel alignment from matching dimensions. | |
| blockers.append("One image has a CRS and the other does not; cannot verify alignment") | |
| alignment_score = 0.0 | |
| checks.update({ | |
| "same_crs": False, | |
| "same_shape": same_shape, | |
| "same_bounds": False, | |
| "alignment_score": alignment_score, | |
| }) | |
| else: | |
| # Neither has CRS — standard pixel-space image analysis (PNG / JPG / etc.) | |
| alignment_score = 1.0 if same_shape else 0.85 | |
| checks.update({ | |
| "same_crs": True, | |
| "same_shape": same_shape, | |
| "same_bounds": True, | |
| "alignment_score": alignment_score, | |
| }) | |
| if not same_shape: | |
| warnings.append(f"Paired images have different dimensions ({a.width}×{a.height} vs {b.width}×{b.height}); auto-normalized") | |
| warnings.append( | |
| "Standard image analysis (no CRS) — coordinates in pixel-space; " | |
| "area measurements will be reported in pixels/percentages" | |
| ) | |
| compatible = not blockers | |
| penalty = 0.20 * len(blockers) + 0.04 * len(warnings) | |
| score = max(0.0, min(1.0, 1.0 - penalty)) | |
| return metadata, QualityReport( | |
| score=round(score, 3), | |
| compatible=compatible, | |
| blockers=blockers, | |
| warnings=warnings, | |
| checks=checks, | |
| ) | |
| def render_preview(path: Path, output_path: Path, max_size: int = 1024) -> Path: | |
| with rasterio.open(path) as src: | |
| ratio = min(1.0, max_size / max(src.width, src.height)) | |
| out_w = max(1, round(src.width * ratio)) | |
| out_h = max(1, round(src.height * ratio)) | |
| indexes = [1, 2, 3] if src.count >= 3 else [1] | |
| data = src.read(indexes, out_shape=(len(indexes), out_h, out_w), resampling=Resampling.bilinear) | |
| if len(indexes) == 1: | |
| gray = (_normalize(data[0]) * 255).astype("uint8") | |
| rgb = np.stack([gray, gray, gray], axis=-1) | |
| else: | |
| rgb = np.stack([(_normalize(data[i]) * 255).astype("uint8") for i in range(3)], axis=-1) | |
| Image.fromarray(rgb, mode="RGB").save(output_path, format="PNG") | |
| return output_path | |
| def _read_gray(path: Path, max_size: int = 1024) -> tuple[np.ndarray, Affine, rasterio.crs.CRS | None]: | |
| with rasterio.open(path) as src: | |
| ratio = min(1.0, max_size / max(src.width, src.height)) | |
| out_w = max(1, round(src.width * ratio)) | |
| out_h = max(1, round(src.height * ratio)) | |
| data = src.read(1, out_shape=(out_h, out_w), resampling=Resampling.bilinear) | |
| transform = src.transform * Affine.scale(src.width / out_w, src.height / out_h) | |
| return _normalize(data), transform, src.crs | |
| def _denoise_binary(mask: np.ndarray) -> np.ndarray: | |
| from scipy.ndimage import uniform_filter | |
| neighbour_count = uniform_filter(mask.astype("float32"), size=3, mode="constant", cval=0.0) * 9 | |
| return mask & (neighbour_count >= 4) | |
| def _area_square_meters(geom: Any, crs: rasterio.crs.CRS | None) -> float | None: | |
| polygon = shape(geom) | |
| if polygon.is_empty or crs is None: | |
| return None | |
| transformer = Transformer.from_crs(crs, "EPSG:6933", always_xy=True) | |
| projected = shapely_transform(transformer.transform, polygon) | |
| return abs(float(projected.area)) | |
| def deterministic_change_detection(before: Path, after: Path, output_dir: Path) -> dict[str, Any]: | |
| """Create a generic radiometric change mask. | |
| This fallback detects changed pixels only. It does not claim what semantic | |
| class changed or whether built-up/water/vegetation increased. | |
| """ | |
| a, transform, crs = _read_gray(before) | |
| b, _, _ = _read_gray(after) | |
| if a.shape != b.shape: | |
| raise ValueError("Rasters must be aligned to the same sampled grid") | |
| delta = np.abs(a - b) | |
| median = float(np.median(delta)) | |
| mad = float(np.median(np.abs(delta - median))) | |
| threshold = max(0.08, median + 3.0 * max(mad, 0.01)) | |
| mask = _denoise_binary(delta > threshold) | |
| changed_pixels = int(mask.sum()) | |
| total_pixels = int(mask.size) | |
| changed_percent = (changed_pixels / total_pixels * 100) if total_pixels else 0.0 | |
| rgba = np.zeros((*mask.shape, 4), dtype="uint8") | |
| rgba[mask] = [239, 68, 68, 190] | |
| mask_path = output_dir / "change_mask.png" | |
| Image.fromarray(rgba, mode="RGBA").save(mask_path, format="PNG") | |
| features = [] | |
| total_area_m2 = 0.0 | |
| area_available = True | |
| for geom, value in shapes(mask.astype("uint8"), mask=mask, transform=transform): | |
| if int(value) != 1: | |
| continue | |
| area = _area_square_meters(geom, crs) | |
| if area is None: | |
| area_available = False | |
| else: | |
| total_area_m2 += area | |
| features.append({ | |
| "type": "Feature", | |
| "geometry": geom, | |
| "properties": {"kind": "generic_change", "area_m2": round(area, 3) if area is not None else None}, | |
| }) | |
| geojson = { | |
| "type": "FeatureCollection", | |
| "features": features, | |
| "properties": {"crs": crs.to_string() if crs else None, "semantic_class": "unknown"}, | |
| } | |
| geojson_path = output_dir / "change_regions.geojson" | |
| geojson_path.write_text(json.dumps(geojson, indent=2), encoding="utf-8") | |
| return { | |
| "changed_pixels": changed_pixels, | |
| "total_pixels": total_pixels, | |
| "changed_percent": round(changed_percent, 3), | |
| "threshold": round(threshold, 4), | |
| "region_count": len(features), | |
| "area_m2": round(total_area_m2, 3) if area_available else None, | |
| "semantic_supported": False, | |
| "mask_path": mask_path, | |
| "geojson_path": geojson_path, | |
| } | |