Blame for app/services/raster.py

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,
    }